t_nlps_master

T_NLPS_MASTER  Tests of NLP solvers via NLPS_MASTER().

0001 function t_nlps_master(quiet)
0002 %T_NLPS_MASTER  Tests of NLP solvers via NLPS_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 %%  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 = 54;
0037 
0038 t_begin(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 [x, f, s, out, lam] = nlps_master(f_fcn, x0, [], [], [], [], [], [], [], opt);
0078 t_ok(s > 0, [t 'success']);
0079 t_is(x, [1; 1], 8, [t 'x']);
0080 t_is(f, 0, 13, [t 'f']);
0081 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0082 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0083 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0084 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0085 
0086 t = sprintf('%s - unconstrained 3-d quadratic : ', name);
0087 %% from http://www.akiti.ca/QuadProgEx0Constr.html
0088 f_fcn = @(x)f3(x);
0089 x0 = [0; 0; 0];
0090 [x, f, s, out, lam] = nlps_master(f_fcn, x0, [], [], [], [], [], [], [], opt);
0091 t_ok(s > 0, [t 'success']);
0092 t_is(x, [3; 5; 7], 6, [t 'x']);
0093 t_is(f, -244, 12, [t 'f']);
0094 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0095 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0096 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0097 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0098 
0099 t = sprintf('%s - constrained 4-d QP : ', name);
0100 %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0101 f_fcn = @(x)f4(x);
0102 x0 = [1; 0; 0; 1];
0103 A = [   1       1       1       1;
0104         0.17    0.11    0.10    0.18    ];
0105 l = [1; 0.10];
0106 u = [1; Inf];
0107 xmin = zeros(4,1);
0108 [x, f, s, out, lam] = nlps_master(f_fcn, x0, A, l, u, xmin, [], [], [], opt);
0109 t_ok(s > 0, [t 'success']);
0110 t_is(x, [0; 2.8; 0.2; 0]/3, 6, [t 'x']);
0111 t_is(f, 3.29/3, 6, [t 'f']);
0112 t_is(lam.mu_l, [6.58;0]/3, 5, [t 'lam.mu_l']);
0113 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0114 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0115 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0116 
0117 t = sprintf('%s - constrained 2-d nonlinear : ', name);
0118 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0119 f_fcn = @(x)f5(x);
0120 gh_fcn = @(x)gh5(x);
0121 hess_fcn = @(x, lam, cost_mult)hess5(x, lam, cost_mult);
0122 x0 = [1.1; 0];
0123 xmin = zeros(2, 1);
0124 % xmax = 3 * ones(2, 1);
0125 [x, f, s, out, lam] = nlps_master(f_fcn, x0, [], [], [], xmin, [], gh_fcn, hess_fcn, opt);
0126 t_ok(s > 0, [t 'success']);
0127 t_is(x, [1; 1], 6, [t 'x']);
0128 t_is(f, -2, 6, [t 'f']);
0129 t_is(lam.ineqnonlin, [0;0.5], 6, [t 'lam.ineqnonlin']);
0130 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0131 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0132 t_is(lam.lower, zeros(size(x)), 9, [t 'lam.lower']);
0133 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0134 
0135 t = sprintf('%s - constrained 3-d nonlinear : ', name);
0136 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0137 f_fcn = @(x)f6(x);
0138 gh_fcn = @(x)gh6(x);
0139 hess_fcn = @(x, lam, cost_mult)hess6(x, lam, cost_mult);
0140 x0 = [1; 1; 0];
0141 [x, f, s, out, lam] = nlps_master(f_fcn, x0, [], [], [], [], [], gh_fcn, hess_fcn, opt);
0142 t_ok(s > 0, [t 'success']);
0143 t_is(x, [1.58113883; 2.23606798; 1.58113883], 6, [t 'x']);
0144 t_is(f, -5*sqrt(2), 6, [t 'f']);
0145 t_is(lam.ineqnonlin, [0;sqrt(2)/2], 7, [t 'lam.ineqnonlin']);
0146 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0147 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0148 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0149 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0150 
0151 t = sprintf('%s - constrained 3-d nonlinear (struct) : ', name);
0152 p = struct('f_fcn', f_fcn, 'x0', x0, 'gh_fcn', gh_fcn, 'hess_fcn', hess_fcn, 'opt', opt);
0153 [x, f, s, out, lam] = nlps_master(p);
0154 t_ok(s > 0, [t 'success']);
0155 t_is(x, [1.58113883; 2.23606798; 1.58113883], 6, [t 'x']);
0156 t_is(f, -5*sqrt(2), 6, [t 'f']);
0157 t_is(lam.ineqnonlin, [0;sqrt(2)/2], 7, [t 'lam.ineqnonlin']);
0158 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0159 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0160 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0161 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0162 
0163 t = sprintf('%s - constrained 4-d nonlinear : ', name);
0164 %% Hock & Schittkowski test problem #71
0165 f_fcn = @(x)f7(x);
0166 gh_fcn = @(x)gh7(x);
0167 hess_fcn = @(x, lam, sigma)hess7(x, lam, sigma);
0168 x0 = [1; 5; 5; 1];
0169 xmin = ones(4, 1);
0170 xmax = 5 * xmin;
0171 [x, f, s, out, lam] = nlps_master(f_fcn, x0, [], [], [], xmin, xmax, gh_fcn, hess_fcn, opt);
0172 t_ok(s > 0, [t 'success']);
0173 t_is(x, [1; 4.7429994; 3.8211503; 1.3794082], 6, [t 'x']);
0174 t_is(f, 17.0140173, 6, [t 'f']);
0175 t_is(lam.eqnonlin, 0.1614686, 5, [t 'lam.eqnonlin']);
0176 t_is(lam.ineqnonlin, 0.55229366, 5, [t 'lam.ineqnonlin']);
0177 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0178 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0179 t_is(lam.lower, [1.08787121024; 0; 0; 0], 5, [t 'lam.lower']);
0180 t_is(lam.upper, zeros(size(x)), 7, [t 'lam.upper']);
0181     end
0182 end
0183 
0184 t_end;
0185 
0186 
0187 % %%-----  eg99 : linearly constrained fmincon example, mips can't solve  -----
0188 % function [f, df, d2f] = eg99(x)
0189 % f = -x(1)*x(2)*x(3);
0190 % df = -[ x(2)*x(3);
0191 %         x(1)*x(3);
0192 %         x(1)*x(2)   ];
0193 % d2f = -[    0       x(3)    x(2);
0194 %             x(3)    0       x(1);
0195 %             x(2)    x(1)    0   ];
0196 % end
0197 %
0198 % x0 = [10;10;10];
0199 % A = [1 2 2];
0200 % l = 0;
0201 % u = 72;
0202 % fmoptions = optimset('Display', 'testing');
0203 % fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0204 % [x, f, s, out, lam] = fmincon(f_fcn, x0, [-A; A], [-l; u], [], [], [], [], [], fmoptions);
0205 % t_is(x, [24; 12; 12], 13, t);
0206 % t_is(f, -3456, 13, t);
0207 
0208 
0209 %% unconstrained banana function
0210 %% from MATLAB Optimization Toolbox's bandem.m
0211 function [f, df, d2f] = f2(x)
0212     a = 100;
0213     f = a*(x(2)-x(1)^2)^2+(1-x(1))^2;
0214     df = [  4*a*(x(1)^3 - x(1)*x(2)) + 2*x(1)-2;
0215             2*a*(x(2) - x(1)^2)                     ];
0216     d2f = 4*a*[ 3*x(1)^2 - x(2) + 1/(2*a),  -x(1);
0217                 -x(1)                       1/2       ];
0218 
0219 
0220 %% unconstrained 3-d quadratic
0221 %% from http://www.akiti.ca/QuadProgEx0Constr.html
0222 function [f, df, d2f] = f3(x)
0223     H = [5 -2 -1; -2 4 3; -1 3 5];
0224     c = [2; -35; -47];
0225     f = 1/2 * x'*H*x + c'*x + 5;
0226     df = H*x + c;
0227     d2f = H;
0228 
0229 
0230 %% constrained 4-d QP
0231 %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0232 function [f, df, d2f] = f4(x)
0233     H = [   1003.1  4.3     6.3     5.9;
0234             4.3     2.2     2.1     3.9;
0235             6.3     2.1     3.5     4.8;
0236             5.9     3.9     4.8     10  ];
0237     c = zeros(4,1);
0238     f = 1/2 * x'*H*x + c'*x;
0239     df = H*x + c;
0240     d2f = H;
0241 
0242 
0243 %% constrained 2-d nonlinear
0244 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0245 function [f, df, d2f] = f5(x)
0246     c = -[1; 1];
0247     f = c'*x;
0248     df = c;
0249     d2f = zeros(2,2);
0250 
0251 function [h, g, dh, dg] = gh5(x)
0252     [h, dh] = h5(x);
0253     dh = dh';
0254     g = []; dg = [];
0255 
0256 function [h, dh] = h5(x)
0257     h = [ -1 -1; 1 1] * x.^2 + [1; -2];
0258     dh = 2 * [-x(1) -x(2); x(1) x(2)];
0259 
0260 function Lxx = hess5(x, lam, cost_mult)
0261     Lxx = hess5a(x, lam.ineqnonlin);
0262 
0263 function Lxx = hess5a(x, mu)
0264     Lxx = 2*[-1 1]*mu*eye(2);
0265 
0266 %% constrained 3-d nonlinear
0267 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0268 function [f, df, d2f] = f6(x)
0269     f = -x(1)*x(2) - x(2)*x(3);
0270     df = -[x(2); x(1)+x(3); x(2)];
0271     d2f = -[0 1 0; 1 0 1; 0 1 0];
0272 
0273 function [h, g, dh, dg] = gh6(x)
0274     [h, dh] = h6(x);
0275     dh = dh';
0276     g = []; dg = [];
0277 
0278 function [h, dh] = h6(x)
0279     h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0280     dh = 2 * [x(1) -x(2) x(3); x(1) x(2) x(3)];
0281 
0282 function Lxx = hess6(x, lam, cost_mult)
0283     if nargin < 3, cost_mult = 1; end
0284     Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + hess6a(x, lam.ineqnonlin);
0285 
0286 function d2h = hess6a(x, mu)
0287     d2h = [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0288 
0289 
0290 %% constrained 4-d nonlinear
0291 %% Hock & Schittkowski test problem #71
0292 function [f, df, d2f] = f7(x)
0293     f = x(1)*x(4)*sum(x(1:3)) + x(3);
0294     df = [ x(1)*x(4) + x(4)*sum(x(1:3));
0295            x(1)*x(4);
0296            x(1)*x(4) + 1;
0297            x(1)*sum(x(1:3)) ];
0298     d2f = sparse([ 2*x(4)        x(4)   x(4)  2*x(1)+x(2)+x(3);
0299               x(4)               0      0     x(1);
0300               x(4)               0      0     x(1);
0301               2*x(1)+x(2)+x(3)  x(1)  x(1)    0
0302         ]);
0303 
0304 function [h, g, dh, dg] = gh7(x)
0305     [g, dg] = g7(x);
0306     [h, dh] = h7(x);
0307     dg = dg';
0308     dh = dh';
0309 
0310 function [g, dg] = g7(x)
0311     g = sum(x.^2) - 40;
0312     dg = 2*x';
0313 
0314 function [h, dh] = h7(x)
0315     h = -prod(x) + 25;
0316     dh = -(prod(x)./x)';
0317 
0318 function Lxx = hess7(x, lam, sigma)
0319     if nargin < 3, sigma = 1; end
0320     [f, df, d2f] = f7(x);
0321     Lxx = sigma * d2f + hess7g(x, lam.eqnonlin) + hess7h(x, lam.ineqnonlin);
0322 
0323 function d2g = hess7g(x, lam)
0324     d2g = lam*2*speye(4);
0325 
0326 function d2h = hess7h(x, mu)
0327     d2h = -mu*sparse([      0     x(3)*x(4) x(2)*x(4) x(2)*x(3);
0328                         x(3)*x(4)     0     x(1)*x(4) x(1)*x(3);
0329                         x(2)*x(4) x(1)*x(4)     0     x(1)*x(2);
0330                         x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ]);