Description of t_mips

0001 function t_mips_pardiso(quiet)
0002 %T_MIPS_PARDISO  Tests of MIPS NLP solver, using PARDISO as linear solver.
0003 
0004 %   MIPS
0005 %   Copyright (c) 2010-2018, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MIPS.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://github.com/MATPOWER/mips for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 num_tests = 60;
0017 
0018 t_begin(num_tests, quiet);
0019 
0020 if have_pardiso()
0021 
0022 opt = struct('linsolver', 'PARDISO');
0023 
0024 t = 'unconstrained banana function : ';
0025 %% from MATLAB Optimization Toolbox's bandem.m
0026 f_fcn = @(x)f2(x);
0027 x0 = [-1.9; 2];
0028 % [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], [], [], struct('verbose', 2));
0029 [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], [], [], opt);
0030 t_is(s, 1, 13, [t 'success']);
0031 t_is(x, [1; 1], 13, [t 'x']);
0032 t_is(f, 0, 13, [t 'f']);
0033 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0034 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0035 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0036 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0037 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0038 
0039 t = 'unconstrained 3-d quadratic : ';
0040 %% from http://www.akiti.ca/QuadProgEx0Constr.html
0041 f_fcn = @(x)f3(x);
0042 x0 = [0; 0; 0];
0043 % [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], [], [], struct('verbose', 2));
0044 [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], [], [], opt);
0045 t_is(s, 1, 13, [t 'success']);
0046 t_is(x, [3; 5; 7], 13, [t 'x']);
0047 t_is(f, -244, 13, [t 'f']);
0048 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0049 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0050 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0051 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0052 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0053 
0054 t = 'constrained 4-d QP : ';
0055 %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0056 f_fcn = @(x)f4(x);
0057 x0 = [1; 0; 0; 1];
0058 A = [   1       1       1       1;
0059         0.17    0.11    0.10    0.18    ];
0060 l = [1; 0.10];
0061 u = [1; Inf];
0062 xmin = zeros(4,1);
0063 % [x, f, s, out, lam] = mips(f_fcn, x0, A, l, u, xmin, [], [], [], struct('verbose', 2));
0064 [x, f, s, out, lam] = mips(f_fcn, x0, A, l, u, xmin, [], [], [], opt);
0065 t_is(s, 1, 13, [t 'success']);
0066 t_is(x, [0; 2.8; 0.2; 0]/3, 6, [t 'x']);
0067 t_is(f, 3.29/3, 6, [t 'f']);
0068 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0069 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0070 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0071 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0072 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0073 
0074 H = [   1003.1  4.3     6.3     5.9;
0075         4.3     2.2     2.1     3.9;
0076         6.3     2.1     3.5     4.8;
0077         5.9     3.9     4.8     10  ];
0078 c = zeros(4,1);
0079 % %% check with quadprog (for dev testing only)
0080 % [x, f, s, out, lam] = quadprog(H,c,-A(2,:), -0.10, A(1,:), 1, xmin);
0081 % t_is(s, 1, 13, [t 'success']);
0082 % t_is(x, [0; 2.8; 0.2; 0]/3, 6, [t 'x']);
0083 % t_is(f, 3.29/3, 6, [t 'f']);
0084 % t_is(lam.eqlin, -6.58/3, 6, [t 'lam.eqlin']);
0085 % t_is(lam.ineqlin, 0, 13, [t 'lam.ineqlin']);
0086 % t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0087 % t_is(lam.upper, [0;0;0;0], 13, [t 'lam.upper']);
0088 
0089 t = 'constrained 2-d nonlinear : ';
0090 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0091 f_fcn = @(x)f5(x);
0092 gh_fcn = @(x)gh5(x);
0093 hess_fcn = @(x, lam, cost_mult)hess5(x, lam, cost_mult);
0094 x0 = [1.1; 0];
0095 xmin = zeros(2, 1);
0096 % xmax = 3 * ones(2, 1);
0097 % [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], xmin, [], gh_fcn, hess_fcn, struct('verbose', 2));
0098 [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], xmin, [], gh_fcn, hess_fcn, opt);
0099 t_is(s, 1, 13, [t 'success']);
0100 t_is(x, [1; 1], 6, [t 'x']);
0101 t_is(f, -2, 6, [t 'f']);
0102 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0103 t_is(lam.ineqnonlin, [0;0.5], 6, [t 'lam.ineqnonlin']);
0104 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0105 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0106 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0107 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0108 % %% check with fmincon (for dev testing only)
0109 % % fmoptions = optimset('Algorithm', 'interior-point');
0110 % % [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], xmin, [], gh_fcn, fmoptions);
0111 % [x, f, s, out, lam] = fmincon(f_fcn, x0, [], [], [], [], [], [], gh_fcn);
0112 % t_is(s, 1, 13, [t 'success']);
0113 % t_is(x, [1; 1], 4, [t 'x']);
0114 % t_is(f, -2, 6, [t 'f']);
0115 % t_is(lam.ineqnonlin, [0;0.5], 6, [t 'lam.ineqnonlin']);
0116 
0117 t = 'constrained 3-d nonlinear : ';
0118 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0119 f_fcn = @(x)f6(x);
0120 gh_fcn = @(x)gh6(x);
0121 hess_fcn = @(x, lam, cost_mult)hess6(x, lam, cost_mult);
0122 x0 = [1; 1; 0];
0123 % [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], gh_fcn, hess_fcn, struct('verbose', 2, 'comptol', 1e-9));
0124 [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], [], [], gh_fcn, hess_fcn, opt);
0125 t_is(s, 1, 13, [t 'success']);
0126 t_is(x, [1.58113883; 2.23606798; 1.58113883], 6, [t 'x']);
0127 t_is(f, -5*sqrt(2), 6, [t 'f']);
0128 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0129 t_is(lam.ineqnonlin, [0;sqrt(2)/2], 7, [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)), 13, [t 'lam.lower']);
0133 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0134 
0135 t = 'constrained 3-d nonlinear (struct) : ';
0136 p = struct('f_fcn', f_fcn, 'x0', x0, 'gh_fcn', gh_fcn, 'hess_fcn', hess_fcn);
0137 [x, f, s, out, lam] = mips(p);
0138 t_is(s, 1, 13, [t 'success']);
0139 t_is(x, [1.58113883; 2.23606798; 1.58113883], 6, [t 'x']);
0140 t_is(f, -5*sqrt(2), 6, [t 'f']);
0141 t_is(out.hist(end).compcond, 0, 6, [t 'compcond']);
0142 t_is(lam.ineqnonlin, [0;sqrt(2)/2], 7, [t 'lam.ineqnonlin']);
0143 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0144 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0145 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0146 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0147 
0148 t = 'constrained 4-d nonlinear : ';
0149 %% Hock & Schittkowski test problem #71
0150 f_fcn = @(x)f7(x);
0151 gh_fcn = @(x)gh7(x);
0152 hess_fcn = @(x, lam, sigma)hess7(x, lam, sigma);
0153 x0 = [1; 5; 5; 1];
0154 xmin = ones(4, 1);
0155 xmax = 5 * xmin;
0156 % [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], xmin, xmax, gh_fcn, hess_fcn, struct('verbose', 2, 'comptol', 1e-9));
0157 [x, f, s, out, lam] = mips(f_fcn, x0, [], [], [], xmin, xmax, gh_fcn, hess_fcn, opt);
0158 t_is(s, 1, 13, [t 'success']);
0159 t_is(x, [1; 4.7429994; 3.8211503; 1.3794082], 6, [t 'x']);
0160 t_is(f, 17.0140173, 6, [t 'f']);
0161 t_is(lam.eqnonlin, 0.1614686, 5, [t 'lam.eqnonlin']);
0162 t_is(lam.ineqnonlin, 0.55229366, 5, [t 'lam.ineqnonlin']);
0163 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0164 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0165 t_is(lam.lower, [1.08787121024; 0; 0; 0], 5, [t 'lam.lower']);
0166 t_is(lam.upper, zeros(size(x)), 7, [t 'lam.upper']);
0167 
0168 else
0169     t_skip(num_tests, 'PARDISO not available.');
0170 end
0171 
0172 t_end;
0173 
0174 
0175 % %%-----  eg99 : linearly constrained fmincon example, mips can't solve  -----
0176 % function [f, df, d2f] = eg99(x)
0177 % f = -x(1)*x(2)*x(3);
0178 % df = -[ x(2)*x(3);
0179 %         x(1)*x(3);
0180 %         x(1)*x(2)   ];
0181 % d2f = -[    0       x(3)    x(2);
0182 %             x(3)    0       x(1);
0183 %             x(2)    x(1)    0   ];
0184 % end
0185 %
0186 % x0 = [10;10;10];
0187 % A = [1 2 2];
0188 % l = 0;
0189 % u = 72;
0190 % fmoptions = optimset('Display', 'testing');
0191 % fmoptions = optimset(fmoptions, 'Algorithm', 'interior-point');
0192 % [x, f, s, out, lam] = fmincon(f_fcn, x0, [-A; A], [-l; u], [], [], [], [], [], fmoptions);
0193 % t_is(x, [24; 12; 12], 13, t);
0194 % t_is(f, -3456, 13, t);
0195 
0196 
0197 %% unconstrained banana function
0198 %% from MATLAB Optimization Toolbox's bandem.m
0199 function [f, df, d2f] = f2(x)
0200     a = 100;
0201     f = a*(x(2)-x(1)^2)^2+(1-x(1))^2;
0202     df = [  4*a*(x(1)^3 - x(1)*x(2)) + 2*x(1)-2;
0203             2*a*(x(2) - x(1)^2)                     ];
0204     d2f = 4*a*[ 3*x(1)^2 - x(2) + 1/(2*a),  -x(1);
0205                 -x(1)                       1/2       ];
0206 
0207 
0208 %% unconstrained 3-d quadratic
0209 %% from http://www.akiti.ca/QuadProgEx0Constr.html
0210 function [f, df, d2f] = f3(x)
0211     H = [5 -2 -1; -2 4 3; -1 3 5];
0212     c = [2; -35; -47];
0213     f = 1/2 * x'*H*x + c'*x + 5;
0214     df = H*x + c;
0215     d2f = H;
0216 
0217 
0218 %% constrained 4-d QP
0219 %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0220 function [f, df, d2f] = f4(x)
0221     H = [   1003.1  4.3     6.3     5.9;
0222             4.3     2.2     2.1     3.9;
0223             6.3     2.1     3.5     4.8;
0224             5.9     3.9     4.8     10  ];
0225     c = zeros(4,1);
0226     f = 1/2 * x'*H*x + c'*x;
0227     df = H*x + c;
0228     d2f = H;
0229 
0230 
0231 %% constrained 2-d nonlinear
0232 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#2-dimensional_example
0233 function [f, df, d2f] = f5(x)
0234     c = -[1; 1];
0235     f = c'*x;
0236     df = c;
0237     d2f = zeros(2,2);
0238 
0239 function [h, g, dh, dg] = gh5(x)
0240     h = [ -1 -1; 1 1] * x.^2 + [1; -2];
0241     dh = 2 * [-x(1) x(1); -x(2) x(2)];
0242     g = []; dg = [];
0243 
0244 function Lxx = hess5(x, lam, cost_mult)
0245     mu = lam.ineqnonlin;
0246     Lxx = 2*[-1 1]*mu*eye(2);
0247 
0248 
0249 %% constrained 3-d nonlinear
0250 %% from https://en.wikipedia.org/wiki/Nonlinear_programming#3-dimensional_example
0251 function [f, df, d2f] = f6(x)
0252     f = -x(1)*x(2) - x(2)*x(3);
0253     df = -[x(2); x(1)+x(3); x(2)];
0254     d2f = -[0 1 0; 1 0 1; 0 1 0];
0255 
0256 function [h, g, dh, dg] = gh6(x)
0257     h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0258     dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
0259     g = []; dg = [];
0260 
0261 function Lxx = hess6(x, lam, cost_mult)
0262     if nargin < 3, cost_mult = 1; end
0263     mu = lam.ineqnonlin;
0264     Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
0265             [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0266 
0267 
0268 %% constrained 4-d nonlinear
0269 %% Hock & Schittkowski test problem #71
0270 function [f, df, d2f] = f7(x)
0271     f = x(1)*x(4)*sum(x(1:3)) + x(3);
0272     df = [ x(1)*x(4) + x(4)*sum(x(1:3));
0273            x(1)*x(4);
0274            x(1)*x(4) + 1;
0275            x(1)*sum(x(1:3)) ];
0276     d2f = sparse([ 2*x(4)        x(4)   x(4)  2*x(1)+x(2)+x(3);
0277               x(4)               0      0     x(1);
0278               x(4)               0      0     x(1);
0279               2*x(1)+x(2)+x(3)  x(1)  x(1)    0
0280         ]);
0281 
0282 function [h, g, dh, dg] = gh7(x)
0283     g = sum(x.^2) - 40;
0284     h = -prod(x) + 25;
0285     dg = 2*x;
0286     dh = -prod(x)./x;
0287 
0288 function Lxx = hess7(x, lam, sigma)
0289     if nargin < 3, sigma = 1; end
0290     lambda = lam.eqnonlin;
0291     mu     = lam.ineqnonlin;
0292     [f, df, d2f] = f7(x);
0293     Lxx = sigma * d2f + lambda*2*speye(4) - ...
0294        mu*sparse([      0     x(3)*x(4) x(2)*x(4) x(2)*x(3);
0295                     x(3)*x(4)     0     x(1)*x(4) x(1)*x(3);
0296                     x(2)*x(4) x(1)*x(4)     0     x(1)*x(2);
0297                     x(2)*x(3) x(1)*x(3) x(1)*x(2)     0  ]);
0298 
0299 
0300 function TorF = have_pardiso()
0301 TorF = have_pardiso_object() || have_pardiso_legacy();
0302 
0303 function TorF = have_pardiso_object()
0304 TorF = exist('pardiso', 'file') == 2;
0305 if TorF
0306     try
0307         id = 1;
0308         A = sparse([1 2; 3 4]);
0309         b = [1;1];
0310         p = pardiso(id, 1, 0);
0311         p.factorize(id, A);
0312         x = p.solve(id, A, b);
0313         p.free(id);
0314         p.clear();
0315         if any(x ~= [-1; 1])
0316             TorF = 0;
0317         end
0318     catch
0319         TorF = 0;
0320     end
0321 end
0322 
0323 function TorF = have_pardiso_legacy()
0324 TorF = exist('pardisoinit', 'file') == 3 && ...
0325         exist('pardisoreorder', 'file') == 3 && ...
0326         exist('pardisofactor', 'file') == 3 && ...
0327         exist('pardisosolve', 'file') == 3 && ...
0328         exist('pardisofree', 'file') == 3;
0329 if TorF
0330     try
0331         A = sparse([1 2; 3 4]);
0332         b = [1;1];
0333         info = pardisoinit(11, 0);
0334         info = pardisoreorder(A, info, false);
0335         info = pardisofactor(A, info, false);
0336         [x, info] = pardisosolve(A, b, info, false);
0337         pardisofree(info);
0338         if any(x ~= [-1; 1])
0339             TorF = 0;
0340         end
0341     catch
0342         info
0343         TorF = 0;
0344     end
0345 end
t_mips_pardiso

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
