0001 function t_qps_matpower(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'};
0017 names = {'DEFAULT', 'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'linprog/quadprog', 'CPLEX', 'MOSEK', 'Gurobi', 'CLP', 'glpk'};
0018 check = {[], 'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi', 'clp', 'glpk'};
0019 does_qp = [1 1 1 1 1 1 1 1 1 1 0];
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 if strcmp(names{k}, 'DEFAULT') || strcmp(names{k}, 'MIPS') || strcmp(names{k}, 'sc-MIPS')
0053 opt.mips_opt.comptol = 1e-8;
0054 end
0055
0056
0057
0058
0059
0060
0061
0062 if strcmp(names{k}, '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 end
0069 if strcmp(names{k}, 'MOSEK')
0070
0071
0072
0073
0074 mpopt.mosek.gap_tol = 1e-10;
0075
0076
0077
0078
0079 vnum = have_feature('mosek', 'vnum');
0080 if vnum >= 8
0081
0082
0083
0084
0085 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0086 end
0087
0088 opt.mosek_opt = mosek_options([], mpopt);
0089 end
0090
0091 t = sprintf('%s - 3-d LP : ', names{k});
0092
0093 c = [-5; -4; -6];
0094 A = [1 -1 1;
0095 -3 -2 -4;
0096 3 2 0];
0097 l = [-Inf; -42; -Inf];
0098 u = [20; Inf; 30];
0099 xmin = [0; 0; 0];
0100 x0 = [];
0101 [x, f, s, out, lam] = qps_matpower([], c, A, l, u, xmin, [], [], opt);
0102 t_is(s, 1, 12, [t 'success']);
0103 t_is(x, [0; 15; 3], 6, [t 'x']);
0104 t_is(f, -78, 6, [t 'f']);
0105 t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0106 t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0107 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0108 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0109 else
0110 t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0111 t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0112 end
0113
0114 if does_qp(k)
0115 t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0116
0117 H = [5 -2 -1; -2 4 3; -1 3 5];
0118 c = [2; -35; -47];
0119 x0 = [0; 0; 0];
0120 [x, f, s, out, lam] = qps_matpower(H, c, [], [], [], [], [], [], opt);
0121 t_is(s, 1, 12, [t 'success']);
0122 t_is(x, [3; 5; 7], 8, [t 'x']);
0123 t_is(f, -249, 13, [t 'f']);
0124 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0125 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0126 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0127 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0128
0129 t = sprintf('%s - constrained 2-d QP : ', names{k});
0130
0131 H = [ 1 -1;
0132 -1 2 ];
0133 c = [-2; -6];
0134 A = [ 1 1;
0135 -1 2;
0136 2 1 ];
0137 l = [];
0138 u = [2; 2; 3];
0139 xmin = [0; 0];
0140 x0 = [];
0141 [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0142 t_is(s, 1, 12, [t 'success']);
0143 t_is(x, [2; 4]/3, 7, [t 'x']);
0144 t_is(f, -74/9, 6, [t 'f']);
0145 t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0146 t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
0147 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0148 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0149 else
0150 t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0151 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0152 end
0153
0154 t = sprintf('%s - constrained 4-d QP : ', names{k});
0155
0156 H = [ 1003.1 4.3 6.3 5.9;
0157 4.3 2.2 2.1 3.9;
0158 6.3 2.1 3.5 4.8;
0159 5.9 3.9 4.8 10 ];
0160 c = zeros(4,1);
0161 A = [ 1 1 1 1;
0162 0.17 0.11 0.10 0.18 ];
0163 l = [1; 0.10];
0164 u = [1; Inf];
0165 xmin = zeros(4,1);
0166 x0 = [1; 0; 0; 1];
0167 [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0168 t_is(s, 1, 12, [t 'success']);
0169 t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0170 t_is(f, 3.29/3, 6, [t 'f']);
0171 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0172 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0173 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0174 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0175 else
0176 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0177 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0178 end
0179
0180 t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0181 p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0182 [x, f, s, out, lam] = qps_matpower(p);
0183 t_is(s, 1, 12, [t 'success']);
0184 t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0185 t_is(f, 3.29/3, 6, [t 'f']);
0186 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0187 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0188 if strcmp(algs{k}, 'CLP') && ~have_feature('opti_clp')
0189 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0190 else
0191 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0192 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0193 end
0194 else
0195 t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0196 end
0197
0198 t = sprintf('%s - infeasible LP : ', names{k});
0199 p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0200 [x, f, s, out, lam] = qps_matpower(p);
0201 t_ok(s <= 0, [t 'no success']);
0202 end
0203 end
0204
0205 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0206 warning(s1.state, diff_alg_warn_id);
0207 end
0208
0209 t_end;