0001 function t_om_solve_nlps(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 1
0013 quiet = 0;
0014 end
0015
0016
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]) = [];
0031 end
0032
0033
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
0056
0057
0058 opt.fmincon_opt.tol_f = 1e-9;
0059 case 'IPOPT'
0060 opt.ipopt_opt = opts;
0061
0062
0063
0064
0065 opt.ipopt_opt.compl_inf_tol = 1e-9;
0066
0067
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
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
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
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
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
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
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
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
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
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
0301
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
0311
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
0324
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
0347
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
0371
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 ]);