t_qps_master

T_QPS_MASTER  Tests of LP/QP solvers via QPS_MASTER().

0001 function t_qps_master(quiet)
0002 %T_QPS_MASTER  Tests of LP/QP solvers via QPS_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 algs = {'DEFAULT', 'BPMPD', 'MIPS', 250, 'IPOPT', 'OT', 'CPLEX', 'MOSEK', 'GUROBI', 'CLP', 'GLPK', 'OSQP'};
0017 names = {'DEFAULT', 'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'linprog/quadprog', 'CPLEX', 'MOSEK', 'Gurobi', 'CLP', 'glpk', 'OSQP'};
0018 check = {[], 'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi', 'clp', 'glpk', 'osqp'};
0019 does_qp = [1 1 1 1 1 1 1 1 1 1 0 1];
0020 
0021 n = 36;
0022 nqp = 28;
0023 t_begin(n*length(algs), quiet);
0024 
0025 diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
0026 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0027     s1 = warning('query', diff_alg_warn_id);
0028     warning('off', diff_alg_warn_id);
0029 end
0030 
0031 for k = 1:length(algs)
0032     if ~isempty(check{k}) && ~have_feature(check{k})
0033         t_skip(n, sprintf('%s not installed', names{k}));
0034     else
0035         opt = struct('verbose', 0, 'alg', algs{k});
0036         mpopt = struct( ...
0037             'verbose', 0, ...
0038             'opf', struct( ...
0039                 'violation', 5e-6 ), ...
0040             'cplex', struct( ...
0041                 'lpmethod', 0, ...
0042                 'qpmethod', 0, ...
0043                 'opt', 0 ), ...
0044             'mosek', struct( ...
0045                 'lp_alg', 0, ...
0046                 'max_it', 0, ...
0047                 'gap_tol', 0, ...
0048                 'max_time', 0, ...
0049                 'num_threads', 0, ...
0050                 'opt', 0 ) ...
0051         );
0052         switch names{k}
0053             case {'DEFAULT', 'MIPS', 'sc-MIPS'}
0054                 opt.mips_opt.comptol = 1e-8;
0055 %             case 'linprog/quadprog'
0056 %                 opt.verbose = 2;
0057 %                 opt.linprog_opt.Algorithm = 'interior-point';
0058 %                 opt.linprog_opt.Algorithm = 'active-set';
0059 %                 opt.linprog_opt.Algorithm = 'simplex';
0060 %                 opt.linprog_opt.Algorithm = 'dual-simplex';
0061 %             end
0062             case 'CPLEX'
0063 %               alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
0064                 alg = 2;        %% use dual simplex
0065                 mpopt.cplex.lpmethod = alg;
0066                 mpopt.cplex.qpmethod = min([4 alg]);
0067                 opt.cplex_opt = cplex_options([], mpopt);
0068             case 'MOSEK'
0069 %                 sc = mosek_symbcon;
0070 %                 alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
0071 %                 alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
0072 %                 mpopt.mosek.lp_alg = alg;
0073                 mpopt.mosek.gap_tol = 1e-10;
0074 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-10;
0075 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS = 1e-10;
0076 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS = 1e-10;
0077 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-10;
0078                 vnum = have_feature('mosek', 'vnum');
0079                 if vnum >= 8
0080 %                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1e-10;
0081 %                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1e-10;
0082 %                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS = 1e-10;
0083 %                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED = 1e-10;
0084                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0085                 end
0086 %                 opt.verbose = 3;
0087                 opt.mosek_opt = mosek_options([], mpopt);
0088             case 'OSQP'
0089                 opt.osqp_opt.polish = 1;
0090 %                 opt.osqp_opt.alpha = 1;
0091 %                 opt.osqp_opt.eps_abs = 1e-8;
0092 %                 opt.osqp_opt.eps_rel = 1e-10;
0093 %                 opt.osqp_opt.eps_prim_inf = 1e-8;
0094 %                 opt.osqp_opt.eps_dual_inf = 1e-8;
0095         end
0096 
0097         t = sprintf('%s - 3-d LP : ', names{k});
0098         %% based on example from 'doc linprog'
0099         c = [-5; -4; -6];
0100         A = [1 -1  1;
0101              -3  -2  -4;
0102              3  2  0];
0103         l = [-Inf; -42; -Inf];
0104         u = [20; Inf; 30];
0105         xmin = [0; 0; 0];
0106         x0 = [];
0107         [x, f, s, out, lam] = qps_master([], c, A, l, u, xmin, [], [], opt);
0108         t_is(s, 1, 12, [t 'success']);
0109         t_is(x, [0; 15; 3], 6, [t 'x']);
0110         t_is(f, -78, 6, [t 'f']);
0111         t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0112         t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0113         if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0114             t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0115         else
0116             t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0117             t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0118         end
0119 
0120         if does_qp(k)
0121             t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0122             %% from http://www.akiti.ca/QuadProgEx0Constr.html
0123             H = [5 -2 -1; -2 4 3; -1 3 5];
0124             c = [2; -35; -47];
0125             x0 = [0; 0; 0];
0126             [x, f, s, out, lam] = qps_master(H, c, [], [], [], [], [], [], opt);
0127             t_is(s, 1, 12, [t 'success']);
0128             t_is(x, [3; 5; 7], 7, [t 'x']);
0129             t_is(f, -249, 13, [t 'f']);
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 = sprintf('%s - constrained 2-d QP : ', names{k});
0136             %% example from 'doc quadprog'
0137             H = [   1   -1;
0138                     -1  2   ];
0139             c = [-2; -6];
0140             A = [   1   1;
0141                     -1  2;
0142                     2   1   ];
0143             l = [];
0144             u = [2; 2; 3];
0145             xmin = [0; 0];
0146             x0 = [];
0147             [x, f, s, out, lam] = qps_master(H, c, A, l, u, xmin, [], x0, opt);
0148             t_is(s, 1, 12, [t 'success']);
0149             t_is(x, [2; 4]/3, 7, [t 'x']);
0150             t_is(f, -74/9, 6, [t 'f']);
0151             t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0152             t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
0153             if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0154                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0155             else
0156                 t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0157                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0158             end
0159 
0160             t = sprintf('%s - constrained 4-d QP : ', names{k});
0161             %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
0162             H = [   1003.1  4.3     6.3     5.9;
0163                     4.3     2.2     2.1     3.9;
0164                     6.3     2.1     3.5     4.8;
0165                     5.9     3.9     4.8     10  ];
0166             c = zeros(4,1);
0167             A = [   1       1       1       1;
0168                     0.17    0.11    0.10    0.18    ];
0169             l = [1; 0.10];
0170             u = [1; Inf];
0171             xmin = zeros(4,1);
0172             x0 = [1; 0; 0; 1];
0173             [x, f, s, out, lam] = qps_master(H, c, A, l, u, xmin, [], x0, opt);
0174             t_is(s, 1, 12, [t 'success']);
0175             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0176             t_is(f, 3.29/3, 6, [t 'f']);
0177             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0178             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0179             if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0180                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0181             else
0182                 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0183                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0184             end
0185 
0186             t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0187             p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0188             [x, f, s, out, lam] = qps_master(p);
0189             t_is(s, 1, 12, [t 'success']);
0190             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0191             t_is(f, 3.29/3, 6, [t 'f']);
0192             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0193             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0194             if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0195                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0196             else
0197                 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0198                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0199             end
0200         else
0201             t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0202         end
0203 
0204         t = sprintf('%s - infeasible LP : ', names{k});
0205         p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0206         [x, f, s, out, lam] = qps_master(p);
0207         t_ok(s <= 0, [t 'no success']);
0208     end
0209 end
0210 
0211 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0212     warning(s1.state, diff_alg_warn_id);
0213 end
0214 
0215 t_end;