Home > matpower7.1 > mp-opt-model > lib > t > t_om_solve_nleqs.m

t_om_solve_nleqs

PURPOSE ^

T_OM_SOLVE_NLEQS Tests of NLEQ solvers via OM.SOLVE().

SYNOPSIS ^

function t_om_solve_nleqs(quiet)

DESCRIPTION ^

T_OM_SOLVE_NLEQS  Tests of NLEQ solvers via OM.SOLVE().

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_om_solve_nleqs(quiet)
0002 %T_OM_SOLVE_NLEQS  Tests of NLEQ solvers via OM.SOLVE().
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 = 18;
0060 
0061 t_begin(14+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                 om = opt_model;
0090                 om.add_var('x', 2, x0);
0091                 om.add_nln_constraint('f', 2, 1, @f1, []);
0092                 [x, f, e, out, jac] = om.solve(opt);
0093                 t_is(e, 1, 12, [t 'success']);
0094                 t_is(x, [-3; 4], 8, [t 'x']);
0095                 t_is(f, 0, 10, [t 'f']);
0096                 switch alg
0097                     case {'DEFAULT', 'CORE-N'}
0098                         out_alg = 'NEWTON';
0099                     otherwise
0100                         out_alg = alg;
0101                 end
0102                 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0103                 eJ = [1 1; 6 1];
0104                 t_is(jac, eJ, 5.8, [t 'jac']);
0105 
0106                 t = sprintf('%s - 2-d function (max_it) : ', name);
0107                 opt.max_it = 3;
0108                 [x, f, e, out] = om.solve(opt);
0109                 t_is(e, 0, 12, [t 'no success']);
0110                 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0111                 opt = rmfield(opt, 'max_it');
0112 
0113                 t = sprintf('%s - 5-d lin/nonlin function : ', name);
0114                 x0_1 = [-1;0];
0115                 x0_2 = [0;0;0];
0116                 A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
0117                 b2 = [-5; 1; -7];
0118                 x2 = [-2; 1; 3];
0119                 om = opt_model;
0120                 om.add_var('x1', 2, x0_1);
0121                 om.add_var('x2', 3, x0_2);
0122                 om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
0123                 om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
0124                 [x, f, e, out, jac] = om.solve(opt);
0125                 t_is(e, 1, 12, [t 'success']);
0126                 t_is(x, [-3; 4; -2; 1; 3], 8, [t 'x']);
0127                 t_is(f, 0, 10, [t 'f']);
0128                 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0129                 eJ = [[1 1; 6 1] zeros(2, 3);
0130                       zeros(3, 2) A2 ];
0131                 t_is(jac, eJ, 5.8, [t 'jac']);
0132             otherwise
0133                 t_skip(12, sprintf('not implemented for solver ''%s''', alg));
0134         end
0135         
0136         t = sprintf('%s - 2-d function2 (struct) : ', name);
0137         x0 = [1;2];
0138         om = opt_model;
0139         om.add_var('x', 2, x0);
0140         om.add_nln_constraint('f', 2, 1, @f2, []);
0141         [x, f, e] = om.solve(opt);
0142         t_is(e, 1, 12, [t 'success']);
0143         t_is(x, [2; 3], 8, [t 'x']);
0144         t_is(f, 0, 10, [t 'f']);
0145         t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);
0146 
0147         opt.max_it = 3;
0148         t = sprintf('%s - 2-d function2 (max_it) : ', name);
0149         [x, f, e, out] = om.solve(opt);
0150         t_is(e, 0, 12, [t 'no success']);
0151         t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0152         opt = rmfield(opt, 'max_it');
0153     end
0154 end
0155 
0156 t = 'om.soln.';
0157 opt.alg = 'DEFAULT';
0158 x0_1 = [-1;0];
0159 x0_2 = [0;0;0];
0160 A2 = sparse([2 -1 0; -3 1 -2; 0 5 -4]);
0161 b2 = [-5; 1; -7];
0162 x2 = [-2; 1; 3];
0163 om = opt_model;
0164 om.add_var('x1', 2, x0_1);
0165 om.add_var('x2', 3, x0_2);
0166 om.add_nln_constraint('f', 2, 1, @f1, [], {'x1'});
0167 om.add_lin_constraint('Ax_b', A2, b2, b2, {'x2'});
0168 opt.parse_soln = 1;
0169 [x, f, e, out, jac] = om.solve(opt);
0170 t_is(om.soln.x, x, 14, [t 'x']);
0171 t_is(om.soln.f, f, 14, [t 'f']);
0172 t_is(om.soln.eflag, e, 14, [t 'eflag']);
0173 t_ok(strcmp(om.soln.output.alg, out.alg), [t 'output.alg']);
0174 t_is(om.soln.jac, jac, 14, [t 'jac']);
0175 
0176 t = 'om.get_soln(''var'', ''x1'') : ';
0177 t_is(om.get_soln('var', 'x1'), x(1:2), 14, [t 'x1']);
0178 
0179 t = 'om.get_soln(''var'', ''x'', ''x2'') : ';
0180 t_is(om.get_soln('var', 'x', 'x2'), x(3:5), 14, [t 'x2']);
0181 
0182 t = 'om.get_soln(''lin'', ''g'', ''Ax_b'') : ';
0183 g = om.get_soln('lin', 'g', 'Ax_b');
0184 t_is(g{1}, f(3:5), 14, [t 'A * x - u']);
0185 t_is(g{2}, f(3:5), 14, [t 'l - A * x']);
0186 
0187 t = 'om.get_soln(''nle'', ''f'') : ';
0188 g = om.get_soln('nle', 'f');
0189 t_is(g, f(1:2), 14, [t 'f']);
0190 
0191 t = 'om.get_soln(''nle'', {''g'', ''dg''}, ''f'') : ';
0192 [g, dg] = om.get_soln('nle', {'g', 'dg'}, 'f');
0193 t_is(g, f(1:2), 14, [t 'f']);
0194 t_is(dg, jac(1:2, 1:2), 14, [t 'jac']);
0195 
0196 t = 'parse_soln : ';
0197 t_is(om.soln.var.val.x1, om.get_soln('var', 'x1'), 14, [t 'var.val.x1']);
0198 t_is(om.soln.var.val.x2, om.get_soln('var', 'x2'), 14, [t 'var.val.x2']);
0199 
0200 t_end;
0201 
0202 
0203 %% 2-d problem with 2 solutions
0204 %% from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/
0205 function [f, J] = f1(x)
0206 if iscell(x)    %% in case it was part of a varset
0207     x = x{1};
0208 end
0209 f = [  x(1)   + x(2) - 1;
0210       -x(1)^2 + x(2) + 5    ];
0211 if nargout > 1
0212     J = sparse([1 1; -2*x(1) 1]);
0213 end
0214 
0215 %% another 2-d problem
0216 %% from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg
0217 function [f, J] = f2(x)
0218 f = [  x(1)^2 +   x(1)*x(2)   - 10;
0219        x(2)   + 3*x(1)*x(2)^2 - 57  ];
0220 if nargout > 1
0221     J = sparse([    2*x(1)+x(2)    x(1);
0222                     3*x(2)^2       6*x(1)*x(2)+1    ]);
0223 end
0224 
0225 function JJ = jac_approx_fcn2()
0226 %% for use with fast-decoupled Newton's method
0227 J = [7 2; 27 37];
0228 JJ = {J(1,1), J(2,2)};
0229 
0230 function x = x_update_fcn2(x, f)
0231 %% for use with Gauss-Seidel method
0232 x(1) = sqrt(10 - x(1)*x(2));
0233 x(2) = sqrt((57-x(2))/3/x(1));
0234 
0235 function x = newton_update_fcn(x, f, J, lin_solver)
0236 dx = mplinsolve(J, -f, lin_solver);     %% compute update step
0237 x = x + dx;                             %% update x

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005