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

t_om_solve_nlps

PURPOSE ^

T_OM_SOLVE_NLPS Tests of NLP solvers via OM.SOLVE().

SYNOPSIS ^

function t_om_solve_nlps(quiet)

DESCRIPTION ^

T_OM_SOLVE_NLPS  Tests of NLP solvers via OM.SOLVE().

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_om_solve_nlps(quiet)
0002 %T_OM_SOLVE_NLPS  Tests of NLP 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 %%  alg         name        check           opts
0017 cfg = {
0018     {'DEFAULT', 'DEFAULT',  [],             []                          },
0019     {'MIPS',    'MIPS',     [],             []                          },
0020     {'MIPS',    'sc-MIPS',  [],             struct('step_control', 1)   },
0021     {'FMINCON', 'fmincon-2','fmincon_ipm',  struct('alg', 2)            },
0022     {'FMINCON', 'fmincon-3','fmincon_ipm',  struct('alg', 3)            },
0023     {'FMINCON', 'fmincon-4','fmincon_ipm',  struct('alg', 4)            },
0024     {'FMINCON', 'fmincon-5','fmincon_ipm',  struct('alg', 5)            },
0025     {'FMINCON', 'fmincon-6','fmincon_ipm',  struct('alg', 6)            },
0026     {'IPOPT',   'IPOPT',    'ipopt',        []                          },
0027     {'KNITRO',  'Knitro',   'knitro',       []                          },
0028 };
0029 if have_feature('matlab') && have_feature('matlab', 'vnum') <= 8.003
0030     cfg([5;8]) = [];    %% Mac MATLAB 7.14-8.2 crash w/ fmincon alg 3, fails w/6
0031 end
0032 % cfg = {
0033 %     {'KNITRO',  'Knitro',   'knitro',       []                          },
0034 % };
0035 
0036 n = 47;
0037 
0038 t_begin(34+n*length(cfg), quiet);
0039 
0040 for k = 1:length(cfg)
0041     alg   = cfg{k}{1};
0042     name  = cfg{k}{2};
0043     check = cfg{k}{3};
0044     opts  = cfg{k}{4};
0045     if ~isempty(check) && ~have_feature(check)
0046         t_skip(n, sprintf('%s not installed', name));
0047     else
0048         opt = struct('verbose', 0, 'alg', alg);
0049         switch alg
0050             case {'DEFAULT', 'MIPS'}
0051                 opt.mips_opt = opts;
0052                 opt.mips_opt.comptol = 1e-8;
0053             case 'FMINCON'
0054                 opt.fmincon_opt = opts;
0055 %                 opt.fmincon_opt.opts.OptimalityTolerance = 1e-10;
0056 %                 opt.fmincon_opt.opts.TolCon = 1e-8;
0057 %                 opt.fmincon_opt.tol_x = 1e-8;
0058                  opt.fmincon_opt.tol_f = 1e-9;
0059             case 'IPOPT'
0060                 opt.ipopt_opt = opts;
0061 %                 opt.ipopt_opt.acceptable_tol = 1e-11;
0062 %                 opt.ipopt_opt.acceptable_compl_inf_tol = 1e-11;
0063 %                 opt.ipopt_opt.acceptable_constr_viol_tol = 1e-11;
0064 %                 opt.ipopt_opt.dual_inf_tol = 1e-9;
0065                 opt.ipopt_opt.compl_inf_tol   = 1e-9;
0066 %                 opt.ipopt_opt.acceptable_obj_change_tol   = 1e-10;
0067 %                 opt.ipopt_opt.linear_solver   = 'ma57';
0068             case 'KNITRO'
0069                 opt.knitro_opt = opts;
0070                 opt.knitro_opt.tol_f = 1e-8;
0071         end
0072 
0073         t = sprintf('%s - unconstrained banana function : ', name);
0074         %% from MATLAB Optimization Toolbox's bandem.m
0075         f_fcn = @(x)f2(x);
0076         x0 = [-1.9; 2];
0077         om = opt_model;
0078         om.add_var('x', 2, x0);
0079         om.add_nln_cost('f', 1, f_fcn);
0080         [x, f, s, out, lam] = om.solve(opt);
0081         t_ok(s > 0, [t 'success']);
0082         t_is(x, [1; 1], 8, [t 'x']);
0083         t_is(f, 0, 13, [t 'f']);
0084         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0085         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0086         t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0087         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0088 
0089         t = sprintf('%s - unconstrained 3-d quadratic : ', name);
0090         %% from http://www.akiti.ca/QuadProgEx0Constr.html
0091         f_fcn = @(x)f3(x);
0092         x0 = [0; 0; 0];
0093         om = opt_model;
0094         om.add_var('x', 3, x0);
0095         om.add_nln_cost('f', 1, f_fcn);
0096         [x, f, s, out, lam] = om.solve(opt);
0097         t_ok(s > 0, [t 'success']);
0098         t_is(x, [3; 5; 7], 6, [t 'x']);
0099         t_is(f, -244, 12, [t 'f']);
0100         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0101         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0102         t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0103         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0104 
0105         t = sprintf('%s - constrained 4-d QP : ', name);
0106         %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0107         f_fcn = @(x)f4(x);
0108         x0 = [1; 0; 0; 1];
0109         A = [   1       1       1       1;
0110                 0.17    0.11    0.10    0.18    ];
0111         l = [1; 0.10];
0112         u = [1; Inf];
0113         xmin = zeros(4,1);
0114         om = opt_model;
0115         om.add_var('x', 4, x0, xmin);
0116         om.add_nln_cost('f', 1, f_fcn);
0117         om.add_lin_constraint('Ax_b', A, l, u);
0118         [x, f, s, out, lam] = om.solve(opt);
0119         t_ok(s > 0, [t 'success']);
0120         t_is(x, [0; 2.8; 0.2; 0]/3, 6, [t 'x']);
0121         t_is(f, 3.29/3, 6, [t 'f']);
0122         t_is(lam.mu_l, [6.58;0]/3, 5, [t 'lam.mu_l']);
0123         t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0124         t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0125         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0126 
0127         t = sprintf('%s - constrained 2-d nonlinear : ', name);
0128         %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0129         f_fcn = @(x)f5(x);
0130         h_fcn = @(x)h5(x);
0131         d2h_fcn = @(x, lam)hess5a(x, lam);
0132         x0 = [1.1; 0];
0133         xmin = zeros(2, 1);
0134         % xmax = 3 * ones(2, 1);
0135         om = opt_model;
0136         om.add_var('x', 2, x0, xmin);
0137         om.add_nln_cost('f', 1, f_fcn);
0138         om.add_nln_constraint('h', 2, 0, h_fcn, d2h_fcn);
0139         [x, f, s, out, lam] = om.solve(opt);
0140         t_ok(s > 0, [t 'success']);
0141         t_is(x, [1; 1], 6, [t 'x']);
0142         t_is(f, -2, 6, [t 'f']);
0143         t_is(lam.ineqnonlin, [0;0.5], 6, [t 'lam.ineqnonlin']);
0144         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0145         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0146         t_is(lam.lower, zeros(size(x)), 9, [t 'lam.lower']);
0147         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0148 
0149         t = sprintf('%s - constrained 3-d nonlinear : ', name);
0150         %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0151         f_fcn = @(x)f6(x);
0152         h_fcn = @(x)h6(x);
0153         d2h_fcn = @(x, lam)hess6a(x, lam);
0154         x0 = [1; 1; 0];
0155         om = opt_model;
0156         om.add_var('x', 3, x0);
0157         om.add_nln_cost('f', 1, f_fcn);
0158         om.add_nln_constraint('h', 2, 0, h_fcn, d2h_fcn);
0159         [x, f, s, out, lam] = om.solve(opt);
0160         t_ok(s > 0, [t 'success']);
0161         t_is(x, [1.58113883; 2.23606798; 1.58113883], 6, [t 'x']);
0162         t_is(f, -5*sqrt(2), 6, [t 'f']);
0163         t_is(lam.ineqnonlin, [0;sqrt(2)/2], 7, [t 'lam.ineqnonlin']);
0164         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0165         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0166         t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0167         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0168         t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);
0169 
0170         t = sprintf('%s - constrained 4-d nonlinear : ', name);
0171         opt.parse_soln = 1;
0172         %% Hock & Schittkowski test problem #71
0173         f_fcn = @(x)f7(x);
0174         g_fcn = @(x)g7(x);
0175         h_fcn = @(x)h7(x);
0176         d2g_fcn = @(x, lam)hess7g(x, lam);
0177         d2h_fcn = @(x, lam)hess7h(x, lam);
0178         x0 = [1; 5; 5; 1];
0179         xmin = ones(4, 1);
0180         xmax = 5 * xmin;
0181         om = opt_model;
0182         om.add_var('x', 4, x0, xmin, xmax);
0183         om.add_nln_cost('f', 1, f_fcn);
0184         om.add_nln_constraint('g', 1, 1, g_fcn, d2g_fcn);
0185         om.add_nln_constraint('h', 1, 0, h_fcn, d2h_fcn);
0186         [x, f, s, out, lam] = om.solve(opt);
0187         t_ok(s > 0, [t 'success']);
0188         t_is(x, [1; 4.7429994; 3.8211503; 1.3794082], 6, [t 'x']);
0189         t_is(f, 17.0140173, 6, [t 'f']);
0190         t_is(lam.eqnonlin, 0.1614686, 5, [t 'lam.eqnonlin']);
0191         t_is(lam.ineqnonlin, 0.55229366, 5, [t 'lam.ineqnonlin']);
0192         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0193         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0194         t_is(lam.lower, [1.08787121024; 0; 0; 0], 5, [t 'lam.lower']);
0195         t_is(lam.upper, zeros(size(x)), 7, [t 'lam.upper']);
0196         opt.parse_soln = 0;
0197     end
0198 end
0199 
0200 t = 'om.soln.';
0201 t_is(om.soln.x, x, 14, [t 'x']);
0202 t_is(om.soln.f, f, 14, [t 'f']);
0203 t_is(om.soln.eflag, s, 14, [t 'eflag']);
0204 t_ok(strcmp(om.soln.output.alg, out.alg), [t 'output.alg']);
0205 t_is(om.soln.lambda.lower, lam.lower, 14, [t 'om.soln.lambda.lower']);
0206 t_is(om.soln.lambda.upper, lam.upper, 14, [t 'om.soln.lambda.upper']);
0207 t_is(om.soln.lambda.mu_l, lam.mu_l, 14, [t 'om.soln.lambda.mu_l']);
0208 t_is(om.soln.lambda.mu_u, lam.mu_u, 14, [t 'om.soln.lambda.mu_u']);
0209 t_is(om.soln.lambda.eqnonlin, lam.eqnonlin, 14, [t 'om.soln.lambda.eqnonlin']);
0210 t_is(om.soln.lambda.ineqnonlin, lam.ineqnonlin, 14, [t 'om.soln.lambda.ineqnonlin']);
0211 
0212 t = 'om.get_soln(''var'', ''x'') : ';
0213 [x1, mu_l, mu_u] = om.get_soln('var', 'x');
0214 t_is(x1, x, 14, [t 'x']);
0215 t_is(mu_l, lam.lower, 14, [t 'mu_l']);
0216 t_is(mu_u, lam.upper, 14, [t 'mu_u']);
0217 
0218 t = 'om.get_soln(''var'', ''mu_l'', ''x'') : ';
0219 t_is(om.get_soln('var', 'mu_l', 'x'), lam.lower, 14, [t 'mu_l']);
0220 
0221 t = 'om.get_soln(''nle'', ''g'') : ';
0222 [g, lm, dg] = om.get_soln('nle', 'g');
0223 [eg, edg] = g_fcn(x);
0224 t_is(g, eg, 14, [t 'g']);
0225 t_is(dg, edg, 14, [t 'dg']);
0226 t_is(lm, lam.eqnonlin, 14, [t 'lam']);
0227 
0228 t = 'om.get_soln(''nle'', {''lam'', ''g''}, ''g'') : ';
0229 [lm, g] = om.get_soln('nle', {'lam', 'g'}, 'g');
0230 t_is(g, eg, 14, [t 'g']);
0231 t_is(dg, edg, 14, [t 'dg']);
0232 t_is(lm, lam.eqnonlin, 14, [t 'lam']);
0233 
0234 t = 'om.get_soln(''nli'', ''h'') : ';
0235 [h, mu, dh] = om.get_soln('nli', 'h');
0236 [eh, edh] = h_fcn(x);
0237 t_is(h, eh, 14, [t 'h']);
0238 t_is(dh, edh, 14, [t 'dh']);
0239 t_is(mu, lam.ineqnonlin, 14, [t 'mu']);
0240 
0241 t = 'om.get_soln(''nli'', {''dh'', ''mu''}, ''h'') : ';
0242 [dh, mu] = om.get_soln('nli', {'dh', 'mu'}, 'h');
0243 t_is(dh, edh, 14, [t 'dh']);
0244 t_is(mu, lam.ineqnonlin, 14, [t 'mu']);
0245 
0246 t = 'om.get_soln(''nlc'', ''f'') : ';
0247 [f1, df, d2f] = om.get_soln('nlc', 'f');
0248 [ef, edf, ed2f] = f_fcn(x);
0249 t_is(f1, f, 14, [t 'f']);
0250 t_is(df, edf, 14, [t 'df']);
0251 t_is(d2f, ed2f, 14, [t 'd2f']);
0252 
0253 t = 'om.get_soln(''nlc'',  ''df'', ''f'') : ';
0254 df = om.get_soln('nlc', 'df', 'f');
0255 t_is(df, edf, 14, [t 'df']);
0256 
0257 t = 'parse_soln : ';
0258 t_is(om.soln.var.val.x, om.get_soln('var', 'x'), 14, [t 'var.val.x']);
0259 t_is(om.soln.var.mu_l.x, om.get_soln('var', 'mu_l', 'x'), 14, [t 'var.mu_l.x']);
0260 t_is(om.soln.var.mu_u.x, om.get_soln('var', 'mu_u', 'x'), 14, [t 'var.mu_u.x']);
0261 t_is(om.soln.nle.lam.g, om.get_soln('nle', 'lam', 'g'), 14, [t 'nle.lam.g']);
0262 t_is(om.soln.nli.mu.h, om.get_soln('nli', 'mu', 'h'), 14, [t 'nli.mu.h']);
0263 
0264 t_end;
0265 
0266 
0267 % %%-----  eg99 : linearly constrained fmincon example, mips can't solve  -----
0268 % function [f, df, d2f] = eg99(x)
0269 % f = -x(1)*x(2)*x(3);
0270 % df = -[ x(2)*x(3);
0271 %         x(1)*x(3);
0272 %         x(1)*x(2)   ];
0273 % d2f = -[    0       x(3)    x(2);
0274 %             x(3)    0       x(1);
0275 %             x(2)    x(1)    0   ];
0276 % end
0277 %
0278 % x0 = [10;10;10];
0279 % A = [1 2 2];
0280 % l = 0;
0281 % u = 72;
0282 % fmoptions = optimset('Display', 'testing');
0283 % fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0284 % [x, f, s, out, lam] = fmincon(f_fcn, x0, [-A; A], [-l; u], [], [], [], [], [], fmoptions);
0285 % t_is(x, [24; 12; 12], 13, t);
0286 % t_is(f, -3456, 13, t);
0287 
0288 
0289 %% unconstrained banana function
0290 %% from MATLAB Optimization Toolbox's bandem.m
0291 function [f, df, d2f] = f2(x)
0292     a = 100;
0293     f = a*(x(2)-x(1)^2)^2+(1-x(1))^2;
0294     df = [  4*a*(x(1)^3 - x(1)*x(2)) + 2*x(1)-2;
0295             2*a*(x(2) - x(1)^2)                     ];
0296     d2f = 4*a*[ 3*x(1)^2 - x(2) + 1/(2*a),  -x(1);
0297                 -x(1)                       1/2       ];
0298 
0299 
0300 %% unconstrained 3-d quadratic
0301 %% from http://www.akiti.ca/QuadProgEx0Constr.html
0302 function [f, df, d2f] = f3(x)
0303     H = [5 -2 -1; -2 4 3; -1 3 5];
0304     c = [2; -35; -47];
0305     f = 1/2 * x'*H*x + c'*x + 5;
0306     df = H*x + c;
0307     d2f = H;
0308 
0309 
0310 %% constrained 4-d QP
0311 %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0312 function [f, df, d2f] = f4(x)
0313     H = [   1003.1  4.3     6.3     5.9;
0314             4.3     2.2     2.1     3.9;
0315             6.3     2.1     3.5     4.8;
0316             5.9     3.9     4.8     10  ];
0317     c = zeros(4,1);
0318     f = 1/2 * x'*H*x + c'*x;
0319     df = H*x + c;
0320     d2f = H;
0321 
0322 
0323 %% constrained 2-d nonlinear
0324 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0325 function [f, df, d2f] = f5(x)
0326     c = -[1; 1];
0327     f = c'*x;
0328     df = c;
0329     d2f = zeros(2,2);
0330 
0331 function [h, g, dh, dg] = gh5(x)
0332     [h, dh] = h5(x);
0333     dh = dh';
0334     g = []; dg = [];
0335 
0336 function [h, dh] = h5(x)
0337     h = [ -1 -1; 1 1] * x.^2 + [1; -2];
0338     dh = 2 * [-x(1) -x(2); x(1) x(2)];
0339 
0340 function Lxx = hess5(x, lam, cost_mult)
0341     Lxx = hess5a(x, lam.ineqnonlin);
0342 
0343 function Lxx = hess5a(x, mu)
0344     Lxx = 2*[-1 1]*mu*eye(2);
0345 
0346 %% constrained 3-d nonlinear
0347 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0348 function [f, df, d2f] = f6(x)
0349     f = -x(1)*x(2) - x(2)*x(3);
0350     df = -[x(2); x(1)+x(3); x(2)];
0351     d2f = -[0 1 0; 1 0 1; 0 1 0];
0352 
0353 function [h, g, dh, dg] = gh6(x)
0354     [h, dh] = h6(x);
0355     dh = dh';
0356     g = []; dg = [];
0357 
0358 function [h, dh] = h6(x)
0359     h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0360     dh = 2 * [x(1) -x(2) x(3); x(1) x(2) x(3)];
0361 
0362 function Lxx = hess6(x, lam, cost_mult)
0363     if nargin < 3, cost_mult = 1; end
0364     Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + hess6a(x, lam.ineqnonlin);
0365 
0366 function d2h = hess6a(x, mu)
0367     d2h = [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0368 
0369 
0370 %% constrained 4-d nonlinear
0371 %% Hock & Schittkowski test problem #71
0372 function [f, df, d2f] = f7(x)
0373     f = x(1)*x(4)*sum(x(1:3)) + x(3);
0374     df = [ x(1)*x(4) + x(4)*sum(x(1:3));
0375            x(1)*x(4);
0376            x(1)*x(4) + 1;
0377            x(1)*sum(x(1:3)) ];
0378     d2f = sparse([ 2*x(4)        x(4)   x(4)  2*x(1)+x(2)+x(3);
0379               x(4)               0      0     x(1);
0380               x(4)               0      0     x(1);
0381               2*x(1)+x(2)+x(3)  x(1)  x(1)    0
0382         ]);
0383 
0384 function [h, g, dh, dg] = gh7(x)
0385     [g, dg] = g7(x);
0386     [h, dh] = h7(x);
0387     dg = dg';
0388     dh = dh';
0389 
0390 function [g, dg] = g7(x)
0391     g = sum(x.^2) - 40;
0392     dg = 2*x';
0393 
0394 function [h, dh] = h7(x)
0395     h = -prod(x) + 25;
0396     dh = -(prod(x)./x)';
0397 
0398 function Lxx = hess7(x, lam, sigma)
0399     if nargin < 3, sigma = 1; end
0400     [f, df, d2f] = f7(x);
0401     Lxx = sigma * d2f + hess7g(x, lam.eqnonlin) + hess7h(x, lam.ineqnonlin);
0402 
0403 function d2g = hess7g(x, lam)
0404     d2g = lam*2*speye(4);
0405 
0406 function d2h = hess7h(x, mu)
0407     d2h = -mu*sparse([      0     x(3)*x(4) x(2)*x(4) x(2)*x(3);
0408                         x(3)*x(4)     0     x(1)*x(4) x(1)*x(3);
0409                         x(2)*x(4) x(1)*x(4)     0     x(1)*x(2);
0410                         x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ]);

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