0001 function t_miqps_master(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 1
0013 quiet = 0;
0014 end
0015
0016 algs = {'DEFAULT', 'CPLEX', 'MOSEK', 'GUROBI', 'GLPK', 'OT'};
0017 names = {'DEFAULT', 'CPLEX', 'MOSEK', 'Gurobi', 'glpk', 'intlin/lin/quadprog'};
0018 check = {@have_miqp_solver, 'cplex', 'mosek', 'gurobi', 'glpk', 'intlinprog'};
0019 does_qp = [0 1 1 1 0 0];
0020 if have_feature('gurobi') || have_feature('cplex') || have_feature('mosek')
0021 does_qp(1) = 1;
0022 end
0023
0024 n = 48;
0025 nqp = 28;
0026 nmiqp = 6;
0027 t_begin(n*length(algs), quiet);
0028
0029 diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
0030 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0031 s1 = warning('query', diff_alg_warn_id);
0032 warning('off', diff_alg_warn_id);
0033 end
0034
0035 for k = 1:length(algs)
0036 if ~isempty(check{k}) && ...
0037 (ischar(check{k}) && ~have_feature(check{k}) || ...
0038 isa(check{k}, 'function_handle') && ~check{k}())
0039 t_skip(n, sprintf('%s not installed', names{k}));
0040 else
0041 opt = struct('verbose', 0, 'alg', algs{k});
0042 mpopt = struct( ...
0043 'verbose', 0, ...
0044 'opf', struct( ...
0045 'violation', 5e-6 ), ...
0046 'cplex', struct( ...
0047 'lpmethod', 0, ...
0048 'qpmethod', 0, ...
0049 'opt', 0 ), ...
0050 'mosek', struct( ...
0051 'lp_alg', 0, ...
0052 'max_it', 0, ...
0053 'gap_tol', 0, ...
0054 'max_time', 0, ...
0055 'num_threads', 0, ...
0056 'opt', 0 ) ...
0057 );
0058 switch names{k}
0059 case 'CPLEX'
0060
0061 alg = 2;
0062 mpopt.cplex.lpmethod = alg;
0063 mpopt.cplex.qpmethod = min([4 alg]);
0064 opt.cplex_opt = cplex_options([], mpopt);
0065 case 'MOSEK'
0066
0067
0068
0069
0070 mpopt.mosek.gap_tol = 1e-10;
0071
0072
0073
0074
0075 vnum = have_feature('mosek', 'vnum');
0076 if vnum >= 8
0077
0078
0079
0080
0081 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0082 end
0083
0084 opt.mosek_opt = mosek_options([], mpopt);
0085 end
0086
0087 t = sprintf('%s - 3-d LP : ', names{k});
0088
0089 c = [-5; -4; -6];
0090 A = [1 -1 1;
0091 -3 -2 -4;
0092 3 2 0];
0093 l = [-Inf; -42; -Inf];
0094 u = [20; Inf; 30];
0095 xmin = [0; 0; 0];
0096 x0 = [];
0097 [x, f, s, out, lam] = miqps_master([], c, A, l, u, xmin, [], [], [], opt);
0098 t_is(s, 1, 12, [t 'success']);
0099 t_is(x, [0; 15; 3], 6, [t 'x']);
0100 t_is(f, -78, 6, [t 'f']);
0101 t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0102 t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0103 t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0104 t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0105
0106 if does_qp(k)
0107 t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0108
0109 H = [5 -2 -1; -2 4 3; -1 3 5];
0110 c = [2; -35; -47];
0111 x0 = [0; 0; 0];
0112 [x, f, s, out, lam] = miqps_master(H, c, [], [], [], [], [], [], [], opt);
0113 t_is(s, 1, 12, [t 'success']);
0114 t_is(x, [3; 5; 7], 7, [t 'x']);
0115 t_is(f, -249, 13, [t 'f']);
0116 t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0117 t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0118 t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0119 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0120
0121 t = sprintf('%s - constrained 2-d QP : ', names{k});
0122
0123 H = [ 1 -1;
0124 -1 2 ];
0125 c = [-2; -6];
0126 A = [ 1 1;
0127 -1 2;
0128 2 1 ];
0129 l = [];
0130 u = [2; 2; 3];
0131 xmin = [0; 0];
0132 x0 = [];
0133 [x, f, s, out, lam] = miqps_master(H, c, A, l, u, xmin, [], x0, [], opt);
0134 t_is(s, 1, 12, [t 'success']);
0135 t_is(x, [2; 4]/3, 7, [t 'x']);
0136 t_is(f, -74/9, 6, [t 'f']);
0137 t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0138 t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
0139 t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0140 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0141
0142 t = sprintf('%s - constrained 4-d QP : ', names{k});
0143
0144 H = [ 1003.1 4.3 6.3 5.9;
0145 4.3 2.2 2.1 3.9;
0146 6.3 2.1 3.5 4.8;
0147 5.9 3.9 4.8 10 ];
0148 c = zeros(4,1);
0149 A = [ 1 1 1 1;
0150 0.17 0.11 0.10 0.18 ];
0151 l = [1; 0.10];
0152 u = [1; Inf];
0153 xmin = zeros(4,1);
0154 x0 = [1; 0; 0; 1];
0155 [x, f, s, out, lam] = miqps_master(H, c, A, l, u, xmin, [], x0, [], opt);
0156 t_is(s, 1, 12, [t 'success']);
0157 t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0158 t_is(f, 3.29/3, 6, [t 'f']);
0159 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0160 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0161 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0162 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0163
0164 t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0165 p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0166 [x, f, s, out, lam] = miqps_master(p);
0167 t_is(s, 1, 12, [t 'success']);
0168 t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0169 t_is(f, 3.29/3, 6, [t 'f']);
0170 t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0171 t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0172 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0173 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0174 else
0175 t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0176 end
0177
0178 t = sprintf('%s - infeasible LP : ', names{k});
0179 p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0180 [x, f, s, out, lam] = miqps_master(p);
0181 t_ok(s <= 0, [t 'no success']);
0182
0183
0184 t = sprintf('%s - 2-d MILP : ', names{k});
0185
0186 c = [-2; -3];
0187 A = sparse([195 273; 4 40]);
0188 u = [1365; 140];
0189 xmax = [4; Inf];
0190 vtype = 'I';
0191 p = struct('c', c, 'A', A, 'u', u, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
0192 [x, f, s, out, lam] = miqps_master(p);
0193 t_is(s, 1, 12, [t 'success']);
0194 t_is(x, [4; 2], 12, [t 'x']);
0195 t_is(lam.mu_l, [0; 0], 12, [t 'lam.mu_l']);
0196 t_is(lam.mu_u, [0; 0], 12, [t 'lam.mu_u']);
0197 t_is(lam.lower, [0; 0], 12, [t 'lam.lower']);
0198 t_is(lam.upper, [2; 3], 12, [t 'lam.upper']);
0199
0200 if does_qp(k)
0201 t = sprintf('%s - 4-d MIQP : ', names{k});
0202
0203 H = sparse([ 33 6 0 0;
0204 6 22 11.5 0;
0205 0 11.5 11 0;
0206 0 0 0 0]);
0207 c = [-1 -2 -3 -1]';
0208 Aineq = [-1 1 1 10;
0209 1 -3 1 0];
0210 bineq = [20 30]';
0211 Aeq = [0 1 0 -3.5];
0212 beq = 0;
0213 xmin = [ 0; 0; 0; 2];
0214 xmax = [40; Inf; Inf; 3];
0215 A = sparse([Aeq; Aineq]);
0216 l = [beq; -Inf; -Inf];
0217 u = [beq; bineq];
0218 vtype = 'CCCI';
0219 p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u, ...
0220 'xmin', xmin, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
0221 [x, f, s, out, lam] = miqps_master(p);
0222 t_is(s, 1, 12, [t 'success']);
0223 t_is(x, [7; 7; 0; 2], 7, [t 'x']);
0224 t_is(lam.mu_l, [466; 0; 0], 6, [t 'lam.mu_l']);
0225 t_is(lam.mu_u, [0; 272; 0], 6, [t 'lam.mu_u']);
0226 t_is(lam.lower, [0; 0; 349.5; 4350], 5, [t 'lam.lower']);
0227 t_is(lam.upper, [0; 0; 0; 0], 7, [t 'lam.upper']);
0228 else
0229 t_skip(nmiqp, sprintf('%s does not handle MIQP problems', names{k}));
0230 end
0231
0232 end
0233 end
0234
0235 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0236 warning(s1.state, diff_alg_warn_id);
0237 end
0238
0239 t_end;
0240
0241 function TorF = have_miqp_solver()
0242 TorF = have_feature('cplex') || have_feature('glpk') || ...
0243 have_feature('gurobi') || have_feature('intlinprog') || ...
0244 have_feature('mosek');
0245