t_nleqs_master

T_NLEQS_MASTER  Tests of NLEQ solvers via NLEQS_MASTER().

0001 function t_nleqs_master(quiet)
0002 %T_NLEQS_MASTER  Tests of NLEQ solvers via NLEQS_MASTER().
0003 
0004 %   MP-Opt-Model
0005 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MP-Opt-Model.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 core_sp_newton = struct( ...
0017     'alg',              'NEWTON', ...
0018     'name',             'Newton''s', ...
0019     'default_max_it',   10, ...
0020     'need_jac',         1, ...
0021     'update_fcn',       @(x, f, J)newton_update_fcn(x, f, J, '')  );
0022 core_sp_gs = struct( ...
0023     'alg',              'GS', ...
0024     'name',             'Gauss-Seidel', ...
0025     'default_max_it',   1000, ...
0026     'need_jac',         0, ...
0027     'update_fcn',       @(x, f)x_update_fcn2(x, f)  );
0028 
0029 if have_feature('matlab')
0030     %%  alg         name        check       opts
0031     cfg = {
0032         {'DEFAULT', 'default',  []          []  },
0033         {'NEWTON',  'Newton',   [],         []  },
0034         {'FD',      'fast-decoupled Newton',[],[]  },
0035         {'FSOLVE',  'fsolve-1', 'fsolve',   []  },
0036         {'FSOLVE',  'fsolve-2', 'fsolve',   struct('Algorithm', 'trust-region-dogleg')  },
0037         {'FSOLVE',  'fsolve-3', 'fsolve',   struct('Algorithm', 'trust-region-reflective')  },
0038         {'FSOLVE',  'fsolve-4', 'fsolve',   struct('Algorithm', 'levenberg-marquardt', 'TolX', 1e-11) },
0039         {'GS',      'Gauss-Seidel',[],      []  },
0040         {'CORE-N',  'Newton-CORE',   [],    core_sp_newton  },
0041         {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs  },
0042     };
0043     if have_feature('matlab', 'vnum') <= 7.010
0044         cfg([6]) = [];  %% MATLAB 7.10 does not work w/ fsolve alg 3
0045     end
0046 else    %% octave
0047     %%  alg         name        check       opts
0048     cfg = {
0049         {'DEFAULT', 'default',  []          []  },
0050         {'NEWTON',  'Newton',   [],         []  },
0051         {'FD',      'fast-decoupled Newton',[],[]  },
0052         {'FSOLVE',  'fsolve', 'fsolve',     struct('TolX', 1e-11)  },
0053         {'GS',      'Gauss-Seidel',[],      []  },
0054         {'CORE-N',  'Newton-CORE',   [],    core_sp_newton  },
0055         {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs  },
0056     };
0057 end
0058 
0059 n = 15;
0060 
0061 t_begin(n*length(cfg), quiet);
0062 
0063 for k = 1:length(cfg)
0064     alg   = cfg{k}{1};
0065     name  = cfg{k}{2};
0066     check = cfg{k}{3};
0067     opts  = cfg{k}{4};
0068     if ~isempty(check) && ~have_feature(check)
0069         t_skip(n, sprintf('%s not installed', name));
0070     else
0071         opt = struct('verbose', 0, 'alg', alg, 'tol', 1e-11);
0072         switch alg
0073             case {'DEFAULT', 'NEWTON'}
0074             case {'FD'}
0075                 opt.fd_opt.jac_approx_fcn = @jac_approx_fcn2;
0076             case 'FSOLVE'
0077                 opt.fsolve_opt = opts;
0078             case 'GS'
0079                 opt.gs_opt.x_update_fcn = @(x, f)x_update_fcn2(x, f);
0080             case {'CORE-N', 'CORE-GS'}
0081                 opt.core_sp = opts;
0082                 opt.alg = 'CORE';
0083         end
0084 
0085         switch alg
0086             case {'DEFAULT', 'NEWTON', 'FSOLVE', 'CORE-N'}
0087                 t = sprintf('%s - 2-d function : ', name);
0088                 x0 = [-1;0];
0089                 [x, f, e, out, jac] = nleqs_master(@f1, x0, opt);
0090                 t_is(e, 1, 12, [t 'success']);
0091                 t_is(x, [-3; 4], 8, [t 'x']);
0092                 t_is(f, 0, 10, [t 'f']);
0093                 switch alg
0094                     case {'DEFAULT', 'CORE-N'}
0095                         out_alg = 'NEWTON';
0096                     otherwise
0097                         out_alg = alg;
0098                 end
0099                 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0100                 eJ = [1 1; 6 1];
0101                 t_is(jac, eJ, 5.8, [t 'jac']);
0102 
0103                 t = sprintf('%s - 2-d function (struct) : ', name);
0104                 p = struct('fcn', @f1, 'x0', [1;0], 'opt', opt);
0105                 [x, f, e] = nleqs_master(p);
0106                 t_is(e, 1, 12, [t 'success']);
0107                 t_is(x, [2; -1], 8, [t 'x']);
0108                 t_is(f, 0, 10, [t 'f']);
0109 
0110                 p.opt.max_it = 3;
0111                 t = sprintf('%s - 2-d function (max_it) : ', name);
0112                 [x, f, e, out] = nleqs_master(p);
0113                 t_is(e, 0, 12, [t 'no success']);
0114                 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0115             otherwise
0116                 t_skip(10, sprintf('not implemented for solver ''%s''', alg));
0117         end
0118         
0119         t = sprintf('%s - 2-d function2 (struct) : ', name);
0120         p = struct('fcn', @f2, 'x0', [1;2], 'opt', opt);
0121         [x, f, e] = nleqs_master(p);
0122         t_is(e, 1, 12, [t 'success']);
0123         t_is(x, [2; 3], 8, [t 'x']);
0124         t_is(f, 0, 10, [t 'f']);
0125 
0126         p.opt.max_it = 3;
0127         t = sprintf('%s - 2-d function2 (max_it) : ', name);
0128         [x, f, e, out] = nleqs_master(p);
0129         t_is(e, 0, 12, [t 'no success']);
0130         t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0131     end
0132 end
0133 
0134 t_end;
0135 
0136 
0137 %% 2-d problem with 2 solutions
0138 %% from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/
0139 function [f, J] = f1(x)
0140 f = [  x(1)   + x(2) - 1;
0141       -x(1)^2 + x(2) + 5    ];
0142 if nargout > 1
0143     J = sparse([1 1; -2*x(1) 1]);
0144 end
0145 
0146 %% another 2-d problem
0147 %% from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg
0148 function [f, J] = f2(x)
0149 f = [  x(1)^2 +   x(1)*x(2)   - 10;
0150        x(2)   + 3*x(1)*x(2)^2 - 57  ];
0151 if nargout > 1
0152     J = sparse([    2*x(1)+x(2)    x(1);
0153                     3*x(2)^2       6*x(1)*x(2)+1    ]);
0154 end
0155 
0156 function JJ = jac_approx_fcn2()
0157 %% for use with fast-decoupled Newton's method
0158 J = [7 2; 27 37];
0159 JJ = {J(1,1), J(2,2)};
0160 
0161 function x = x_update_fcn2(x, f)
0162 %% for use with Gauss-Seidel method
0163 x(1) = sqrt(10 - x(1)*x(2));
0164 x(2) = sqrt((57-x(2))/3/x(1));
0165 
0166 function x = newton_update_fcn(x, f, J, lin_solver)
0167 dx = mplinsolve(J, -f, lin_solver);     %% compute update step
0168 x = x + dx;                             %% update x