t_miqps_matpower

T_MIQPS_MATPOWER  Tests of MIQPS_MATPOWER MIQP solvers.

0001 function t_miqps_matpower(quiet)
0002 %T_MIQPS_MATPOWER  Tests of MIQPS_MATPOWER MIQP solvers.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MATPOWER.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://matpower.org for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 algs = {'CPLEX', 'MOSEK', 'GUROBI', 'GLPK', 'OT'};
0017 names = {'CPLEX', 'MOSEK', 'Gurobi', 'glpk', 'intlin/lin/quadprog'};
0018 check = {'cplex', 'mosek', 'gurobi', 'glpk', 'intlinprog'};
0019 does_qp = [1 1 1 0 0];
0020 
0021 n = 48;
0022 nqp = 28;
0023 nmiqp = 6;
0024 t_begin(n*length(algs), quiet);
0025 
0026 diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
0027 if have_fcn('quadprog') && have_fcn('quadprog', 'vnum') == 7.005
0028     s1 = warning('query', diff_alg_warn_id);
0029     warning('off', diff_alg_warn_id);
0030 end
0031 
0032 for k = 1:length(algs)
0033     if ~isempty(check{k}) && ~have_fcn(check{k})
0034         t_skip(n, sprintf('%s not installed', names{k}));
0035     else
0036         opt = struct('verbose', 0, 'alg', algs{k});
0037         if strcmp(names{k}, 'CPLEX')
0038 %           alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
0039             alg = 2;        %% use dual simplex
0040             mpopt = mpoption('cplex.lpmethod', alg, 'cplex.qpmethod', min([4 alg]));
0041             opt.cplex_opt = cplex_options([], mpopt);
0042         end
0043         if strcmp(names{k}, 'MOSEK')
0044             mpopt = mpoption;
0045 %             sc = mosek_symbcon;
0046 %             alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
0047 %             alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
0048 %             mpopt = mpoption(mpopt, 'mosek.lp_alg', alg );
0049             mpopt = mpoption(mpopt, 'mosek.gap_tol', 1e-10);
0050 %             mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS', 1e-10);
0051 %             mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS', 1e-10);
0052 %             mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS', 1e-10);
0053 %             mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP', 1e-10);
0054             vnum = have_fcn('mosek', 'vnum');
0055             if vnum >= 8
0056 %                 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS', 1e-10);
0057 %                 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS', 1e-10);
0058 %                 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS', 1e-10);
0059 %                 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED', 1e-10);
0060                 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP', 1e-10);
0061             end
0062 %             opt.verbose = 3;
0063             opt.mosek_opt = mosek_options([], mpopt);
0064         end
0065 
0066         t = sprintf('%s - 3-d LP : ', names{k});
0067         %% based on example from 'doc linprog'
0068         c = [-5; -4; -6];
0069         A = [1 -1  1;
0070              -3  -2  -4;
0071              3  2  0];
0072         l = [-Inf; -42; -Inf];
0073         u = [20; Inf; 30];
0074         xmin = [0; 0; 0];
0075         x0 = [];
0076         [x, f, s, out, lam] = miqps_matpower([], c, A, l, u, xmin, [], [], [], opt);
0077         t_is(s, 1, 12, [t 'success']);
0078         t_is(x, [0; 15; 3], 6, [t 'x']);
0079         t_is(f, -78, 6, [t 'f']);
0080         t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0081         t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0082         t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0083         t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0084 
0085         if does_qp(k)
0086             t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0087             %% from http://www.akiti.ca/QuadProgEx0Constr.html
0088             H = [5 -2 -1; -2 4 3; -1 3 5];
0089             c = [2; -35; -47];
0090             x0 = [0; 0; 0];
0091             [x, f, s, out, lam] = miqps_matpower(H, c, [], [], [], [], [], [], [], opt);
0092             t_is(s, 1, 12, [t 'success']);
0093             t_is(x, [3; 5; 7], 8, [t 'x']);
0094             t_is(f, -249, 13, [t 'f']);
0095             t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0096             t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0097             t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0098             t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0099         
0100             t = sprintf('%s - constrained 2-d QP : ', names{k});
0101             %% example from 'doc quadprog'
0102             H = [   1   -1;
0103                     -1  2   ];
0104             c = [-2; -6];
0105             A = [   1   1;
0106                     -1  2;
0107                     2   1   ];
0108             l = [];
0109             u = [2; 2; 3];
0110             xmin = [0; 0];
0111             x0 = [];
0112             [x, f, s, out, lam] = miqps_matpower(H, c, A, l, u, xmin, [], x0, [], opt);
0113             t_is(s, 1, 12, [t 'success']);
0114             t_is(x, [2; 4]/3, 7, [t 'x']);
0115             t_is(f, -74/9, 6, [t 'f']);
0116             t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0117             t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
0118             t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0119             t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0120 
0121             t = sprintf('%s - constrained 4-d QP : ', names{k});
0122             %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0123             H = [   1003.1  4.3     6.3     5.9;
0124                     4.3     2.2     2.1     3.9;
0125                     6.3     2.1     3.5     4.8;
0126                     5.9     3.9     4.8     10  ];
0127             c = zeros(4,1);
0128             A = [   1       1       1       1;
0129                     0.17    0.11    0.10    0.18    ];
0130             l = [1; 0.10];
0131             u = [1; Inf];
0132             xmin = zeros(4,1);
0133             x0 = [1; 0; 0; 1];
0134             [x, f, s, out, lam] = miqps_matpower(H, c, A, l, u, xmin, [], x0, [], opt);
0135             t_is(s, 1, 12, [t 'success']);
0136             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0137             t_is(f, 3.29/3, 6, [t 'f']);
0138             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0139             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0140             t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0141             t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0142 
0143             t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0144             p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0145             [x, f, s, out, lam] = miqps_matpower(p);
0146             t_is(s, 1, 12, [t 'success']);
0147             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0148             t_is(f, 3.29/3, 6, [t 'f']);
0149             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0150             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0151             t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0152             t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0153         else
0154             t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0155         end
0156 
0157         t = sprintf('%s - infeasible LP : ', names{k});
0158         p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0159         [x, f, s, out, lam] = miqps_matpower(p);
0160         t_ok(s <= 0, [t 'no success']);
0161 
0162 % opt.verbose = 3;
0163         t = sprintf('%s - 2-d MILP : ', names{k});
0164         %% from MOSEK 6.0 Guided Tour, section  7.13.1, https://docs.mosek.com/6.0/toolbox/node009.html#283040944
0165         c = [-2; -3];
0166         A = sparse([195 273; 4 40]);
0167         u = [1365; 140];
0168         xmax = [4; Inf];
0169         vtype = 'I';
0170         p = struct('c', c, 'A', A, 'u', u, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
0171         [x, f, s, out, lam] = miqps_matpower(p);
0172         t_is(s, 1, 12, [t 'success']);
0173         t_is(x, [4; 2], 12, [t 'x']);
0174         t_is(lam.mu_l, [0; 0], 12, [t 'lam.mu_l']);
0175         t_is(lam.mu_u, [0; 0], 12, [t 'lam.mu_u']);
0176         t_is(lam.lower, [0; 0], 12, [t 'lam.lower']);
0177         t_is(lam.upper, [2; 3], 12, [t 'lam.upper']);
0178 
0179         if does_qp(k)
0180             t = sprintf('%s - 4-d MIQP : ', names{k});
0181             %% from cplexmiqpex.m, CPLEX_Studio_Academic124/cplex/examples/src/matlab/cplexmiqpex.m
0182             H = sparse([ 33   6    0    0;
0183                           6  22   11.5  0;
0184                           0  11.5 11    0;
0185                           0   0    0    0]);
0186             c = [-1 -2 -3 -1]';
0187             Aineq = [-1  1  1 10;
0188                1 -3  1  0];
0189             bineq = [20  30]';
0190             Aeq   = [0  1  0 -3.5];
0191             beq   =  0;
0192             xmin    = [ 0;   0;   0; 2];
0193             xmax    = [40; Inf; Inf; 3];
0194             A = sparse([Aeq; Aineq]);
0195             l = [beq; -Inf; -Inf];
0196             u = [beq; bineq];
0197             vtype = 'CCCI';
0198             p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u, ...
0199                 'xmin', xmin, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
0200             [x, f, s, out, lam] = miqps_matpower(p);
0201             t_is(s, 1, 12, [t 'success']);
0202             t_is(x, [7; 7; 0; 2], 7, [t 'x']);
0203             t_is(lam.mu_l, [466; 0; 0], 6, [t 'lam.mu_l']);
0204             t_is(lam.mu_u, [0; 272; 0], 6, [t 'lam.mu_u']);
0205             t_is(lam.lower, [0; 0; 349.5; 4350], 5, [t 'lam.lower']);
0206             t_is(lam.upper, [0; 0; 0; 0], 7, [t 'lam.upper']);
0207         else
0208             t_skip(nmiqp, sprintf('%s does not handle MIQP problems', names{k}));
0209         end
0210 % opt.verbose = 0;
0211     end
0212 end
0213 
0214 if have_fcn('quadprog') && have_fcn('quadprog', 'vnum') == 7.005
0215     warning(s1.state, diff_alg_warn_id);
0216 end
0217 
0218 t_end;