0001 function t_miqps_matpower(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
0039 alg = 2;
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
0046
0047
0048
0049 mpopt = mpoption(mpopt, 'mosek.gap_tol', 1e-10);
0050
0051
0052
0053
0054 vnum = have_fcn('mosek', 'vnum');
0055 if vnum >= 8
0056
0057
0058
0059
0060 mpopt = mpoption(mpopt, 'mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP', 1e-10);
0061 end
0062
0063 opt.mosek_opt = mosek_options([], mpopt);
0064 end
0065
0066 t = sprintf('%s - 3-d LP : ', names{k});
0067
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
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
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
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
0163 t = sprintf('%s - 2-d MILP : ', names{k});
0164
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
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
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;