0001 function t_om_solve_qps(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
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 = 30;
0022 nqp = 21;
0023 t_begin(27+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
0056
0057
0058
0059
0060
0061
0062 case 'CPLEX'
0063
0064 alg = 2;
0065 mpopt.cplex.lpmethod = alg;
0066 mpopt.cplex.qpmethod = min([4 alg]);
0067 opt.cplex_opt = cplex_options([], mpopt);
0068 case 'MOSEK'
0069
0070
0071
0072
0073 mpopt.mosek.gap_tol = 1e-10;
0074
0075
0076
0077
0078 vnum = have_feature('mosek', 'vnum');
0079 if vnum >= 8
0080
0081
0082
0083
0084 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0085 end
0086
0087 opt.mosek_opt = mosek_options([], mpopt);
0088 case 'OSQP'
0089 opt.osqp_opt.polish = 1;
0090
0091
0092
0093
0094
0095 end
0096
0097 t = sprintf('%s - 3-d LP : ', names{k});
0098
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 om = opt_model;
0108 om.add_var('x', 3, x0, xmin);
0109 om.add_quad_cost('c', [], c);
0110 om.add_lin_constraint('Ax', A, l, u);
0111 [x, f, s, out, lam] = om.solve(opt);
0112 t_is(s, 1, 12, [t 'success']);
0113 t_is(x, [0; 15; 3], 6, [t 'x']);
0114 t_is(f, -78, 6, [t 'f']);
0115 t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0116 t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0117 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0118 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0119 else
0120 t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0121 t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0122 end
0123 t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);
0124
0125 if does_qp(k)
0126 t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0127
0128 H = [5 -2 -1; -2 4 3; -1 3 5];
0129 c = [2; -35; -47];
0130 x0 = [0; 0; 0];
0131 om = opt_model;
0132 om.add_var('x', 3, x0);
0133 om.add_quad_cost('c', H, c);
0134 [x, f, s, out, lam] = om.solve(opt);
0135 t_is(s, 1, 12, [t 'success']);
0136 t_is(x, [3; 5; 7], 7, [t 'x']);
0137 t_is(f, -249, 13, [t 'f']);
0138 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0139 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0140 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0141 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0142
0143 t = sprintf('%s - constrained 2-d QP : ', names{k});
0144
0145 H = [ 1 -1;
0146 -1 2 ];
0147 c = [-2; -6];
0148 A = [ 1 1;
0149 -1 2;
0150 2 1 ];
0151 l = [];
0152 u = [2; 2; 3];
0153 xmin = [0; 0];
0154 x0 = [];
0155 om = opt_model;
0156 om.add_var('x', 2, x0, xmin);
0157 om.add_quad_cost('c', H, c);
0158 om.add_lin_constraint('Ax', A, l, u);
0159 [x, f, s, out, lam] = om.solve(opt);
0160 t_is(s, 1, 12, [t 'success']);
0161 t_is(x, [2; 4]/3, 7, [t 'x']);
0162 t_is(f, -74/9, 6, [t 'f']);
0163 t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0164 t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
0165 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0166 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0167 else
0168 t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0169 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0170 end
0171
0172 t = sprintf('%s - constrained 4-d QP : ', names{k});
0173
0174 H = [ 1003.1 4.3 6.3 5.9;
0175 4.3 2.2 2.1 3.9;
0176 6.3 2.1 3.5 4.8;
0177 5.9 3.9 4.8 10 ];
0178 c = zeros(4,1);
0179 A = [ 1 1 1 1;
0180 0.17 0.11 0.10 0.18 ];
0181 l = [1; 0.10];
0182 u = [1; Inf];
0183 xmin = zeros(4,1);
0184 x0 = [1; 0; 0; 1];
0185 om = opt_model;
0186 om.add_var('x', 4, x0, xmin);
0187 om.add_quad_cost('c', H, c);
0188 om.add_lin_constraint('Ax', A, l, u);
0189 [x, f, s, out, lam] = om.solve(opt);
0190 t_is(s, 1, 12, [t 'success']);
0191 t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0192 t_is(f, 3.29/3, 6, [t 'f']);
0193 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0194 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0195 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0196 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0197 else
0198 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0199 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0200 end
0201 else
0202 t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0203 end
0204
0205 t = sprintf('%s - infeasible LP : ', names{k});
0206 p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0207 [x, f, s, out, lam] = qps_master(p);
0208 t_ok(s <= 0, [t 'no success']);
0209 end
0210 end
0211
0212 t = 'om.soln.';
0213 opt.alg = 'DEFAULT';
0214
0215 H = [ 1003.1 4.3 6.3 5.9;
0216 4.3 2.2 2.1 3.9;
0217 6.3 2.1 3.5 4.8;
0218 5.9 3.9 4.8 10 ];
0219 c = zeros(4,1);
0220 A = [ 1 1 1 1;
0221 0.17 0.11 0.10 0.18 ];
0222 l = [1; 0.10];
0223 u = [1; Inf];
0224 xmin = zeros(4,1);
0225 x0 = [1; 0; 0; 1];
0226 om = opt_model;
0227 om.add_var('x', 4, x0, xmin);
0228 om.add_quad_cost('c', H, c);
0229 om.add_lin_constraint('Ax', A, l, u);
0230 opt.parse_soln = 1;
0231 [x, f, s, out, lam] = om.solve(opt);
0232 t_is(om.soln.x, x, 14, [t 'x']);
0233 t_is(om.soln.f, f, 14, [t 'f']);
0234 t_is(om.soln.eflag, s, 14, [t 'eflag']);
0235 t_ok(strcmp(om.soln.output.alg, out.alg), [t 'output.alg']);
0236 t_is(om.soln.lambda.lower, lam.lower, 14, [t 'om.soln.lambda.lower']);
0237 t_is(om.soln.lambda.upper, lam.upper, 14, [t 'om.soln.lambda.upper']);
0238 t_is(om.soln.lambda.mu_l, lam.mu_l, 14, [t 'om.soln.lambda.mu_l']);
0239 t_is(om.soln.lambda.mu_u, lam.mu_u, 14, [t 'om.soln.lambda.mu_u']);
0240
0241 t = 'om.get_soln(''var'', ''x'') : ';
0242 [x1, mu_l, mu_u] = om.get_soln('var', 'x');
0243 t_is(x1, x, 14, [t 'x']);
0244 t_is(mu_l, lam.lower, 14, [t 'mu_l']);
0245 t_is(mu_u, lam.upper, 14, [t 'mu_u']);
0246
0247 t = 'om.get_soln(''var'', ''mu_l'', ''x'') : ';
0248 t_is(om.get_soln('var', 'mu_l', 'x'), lam.lower, 14, [t 'mu_l']);
0249
0250 t = 'om.get_soln(''lin'', ''Ax'') : ';
0251 [g, mu_l] = om.get_soln('lin', 'Ax');
0252 t_is(g{1}, A*x-u, 14, [t 'A * x - u']);
0253 t_is(g{2}, l-A*x, 14, [t 'l - A * x']);
0254 t_is(mu_l, lam.mu_l, 14, [t 'mu_l']);
0255
0256 t = 'om.get_soln(''lin'', {''mu_u'', ''mu_l'', ''Ax_u''}, ''Ax'') : ';
0257 [mu_u, mu_l, g] = om.get_soln('lin', {'mu_u', 'mu_l', 'Ax_u'}, 'Ax');
0258 t_is(g, A*x-u, 14, [t 'A * x - u']);
0259 t_is(mu_l, lam.mu_l, 14, [t 'mu_l']);
0260 t_is(mu_u, lam.mu_u, 14, [t 'mu_u']);
0261
0262 t = 'om.get_soln(''qdc'', ''c'') : ';
0263 [f1, df, d2f] = om.get_soln('qdc', 'c');
0264 t_is(f1, f, 14, [t 'f']);
0265 t_is(df, H*x+c, 14, [t 'df']);
0266 t_is(d2f, H, 14, [t 'd2f']);
0267
0268 t = 'om.get_soln(''qdc'', ''df'', ''c'') : ';
0269 df = om.get_soln('qdc', 'df', 'c');
0270 t_is(df, H*x+c, 14, [t 'df']);
0271
0272 t = 'parse_soln : ';
0273 t_is(om.soln.var.val.x, om.get_soln('var', 'x'), 14, [t 'var.val.x']);
0274 t_is(om.soln.var.mu_l.x, om.get_soln('var', 'mu_l', 'x'), 14, [t 'var.mu_l.x']);
0275 t_is(om.soln.var.mu_u.x, om.get_soln('var', 'mu_u', 'x'), 14, [t 'var.mu_u.x']);
0276 t_is(om.soln.lin.mu_l.Ax, om.get_soln('lin', 'mu_l', 'Ax'), 14, [t 'lin.mu_l.Ax']);
0277 t_is(om.soln.lin.mu_u.Ax, om.get_soln('lin', 'mu_u', 'Ax'), 14, [t 'lin.mu_u.Ax']);
0278
0279 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0280 warning(s1.state, diff_alg_warn_id);
0281 end
0282
0283 t_end;