0001 function t_om_solve_miqps(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 = 15;
0025 nmiqp = 7;
0026 t_begin(28+n*length(algs), quiet);
0027
0028 diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
0029 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0030 s1 = warning('query', diff_alg_warn_id);
0031 warning('off', diff_alg_warn_id);
0032 end
0033
0034 for k = 1:length(algs)
0035 if ~isempty(check{k}) && ...
0036 (ischar(check{k}) && ~have_feature(check{k}) || ...
0037 isa(check{k}, 'function_handle') && ~check{k}())
0038 t_skip(n, sprintf('%s not installed', names{k}));
0039 else
0040 opt = struct('verbose', 0, 'alg', algs{k});
0041 mpopt = struct( ...
0042 'verbose', 0, ...
0043 'opf', struct( ...
0044 'violation', 5e-6 ), ...
0045 'cplex', struct( ...
0046 'lpmethod', 0, ...
0047 'qpmethod', 0, ...
0048 'opt', 0 ), ...
0049 'mosek', struct( ...
0050 'lp_alg', 0, ...
0051 'max_it', 0, ...
0052 'gap_tol', 0, ...
0053 'max_time', 0, ...
0054 'num_threads', 0, ...
0055 'opt', 0 ) ...
0056 );
0057 if strcmp(names{k}, 'CPLEX')
0058
0059 alg = 2;
0060 mpopt.cplex.lpmethod = alg;
0061 mpopt.cplex.qpmethod = min([4 alg]);
0062 opt.cplex_opt = cplex_options([], mpopt);
0063 end
0064 if strcmp(names{k}, 'MOSEK')
0065
0066
0067
0068
0069 mpopt.mosek.gap_tol = 1e-10;
0070
0071
0072
0073
0074 vnum = have_feature('mosek', 'vnum');
0075 if vnum >= 8
0076
0077
0078
0079
0080 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0081 end
0082
0083 opt.mosek_opt = mosek_options([], mpopt);
0084 end
0085
0086
0087 t = sprintf('%s - 2-d MILP : ', names{k});
0088
0089 c = [-2; -3];
0090 A = sparse([195 273; 4 40]);
0091 u = [1365; 140];
0092 xmax = [4; Inf];
0093 vtype = 'I';
0094 om = opt_model;
0095 om.add_var('x', 2, [], [], xmax, vtype);
0096 om.add_quad_cost('c', [], c);
0097 om.add_lin_constraint('Ax', A, [], u);
0098 [x, f, s, out, lam] = om.solve(opt);
0099 t_is(s, 1, 12, [t 'success']);
0100 t_is(x, [4; 2], 12, [t 'x']);
0101 t_is(f, -14, 12, [t 'f']);
0102 t_is(lam.mu_l, [0; 0], 12, [t 'lam.mu_l']);
0103 t_is(lam.mu_u, [0; 0], 12, [t 'lam.mu_u']);
0104 t_is(lam.lower, [0; 0], 12, [t 'lam.lower']);
0105 t_is(lam.upper, [2; 3], 12, [t 'lam.upper']);
0106 t_ok(~isfield(om.soln, 'var'), [t 'no parse_soln() outputs']);
0107
0108 if does_qp(k)
0109 t = sprintf('%s - 4-d MIQP : ', names{k});
0110
0111 H = sparse([ 33 6 0 0;
0112 6 22 11.5 0;
0113 0 11.5 11 0;
0114 0 0 0 0]);
0115 c = [-1 -2 -3 -1]';
0116 Aineq = [-1 1 1 10;
0117 1 -3 1 0];
0118 bineq = [20 30]';
0119 Aeq = [0 1 0 -3.5];
0120 beq = 0;
0121 xmin = [ 0; 0; 0; 2];
0122 xmax = [40; Inf; Inf; 3];
0123 A = sparse([Aeq; Aineq]);
0124 l = [beq; -Inf; -Inf];
0125 u = [beq; bineq];
0126 vtype = 'CCCI';
0127 om = opt_model;
0128 om.add_var('x', 4, [], xmin, xmax, vtype);
0129 om.add_quad_cost('c', H, c);
0130 om.add_lin_constraint('Ax', A, l, u);
0131 [x, f, s, out, lam] = om.solve(opt);
0132 t_is(s, 1, 12, [t 'success']);
0133 t_is(x, [7; 7; 0; 2], 7, [t 'x']);
0134 t_is(f, 1618.5, 5, [t 'f']);
0135 t_is(lam.mu_l, [466; 0; 0], 6, [t 'lam.mu_l']);
0136 t_is(lam.mu_u, [0; 272; 0], 6, [t 'lam.mu_u']);
0137 t_is(lam.lower, [0; 0; 349.5; 4350], 5, [t 'lam.lower']);
0138 t_is(lam.upper, [0; 0; 0; 0], 7, [t 'lam.upper']);
0139 else
0140 t_skip(nmiqp, sprintf('%s does not handle MIQP problems', names{k}));
0141 end
0142
0143 end
0144 end
0145
0146 if have_miqp_solver()
0147 t = 'om.soln.';
0148 c = [-2; -3];
0149 A = sparse([195 273; 4 40]);
0150 u = [1365; 140];
0151 xmax = [4; Inf];
0152 vtype = 'I';
0153 om = opt_model;
0154 om.add_var('x', 2, [], [], xmax, vtype);
0155 om.add_quad_cost('c', [], c);
0156 om.add_lin_constraint('Ax', A, [], u);
0157 opt.parse_soln = 1;
0158 [x, f, s, out, lam] = om.solve(opt);
0159 t_is(om.soln.x, x, 14, [t 'x']);
0160 t_is(om.soln.f, f, 14, [t 'f']);
0161 t_is(om.soln.eflag, s, 14, [t 'eflag']);
0162 t_ok(strcmp(om.soln.output.alg, out.alg), [t 'output.alg']);
0163 t_is(om.soln.lambda.lower, lam.lower, 14, [t 'om.soln.lambda.lower']);
0164 t_is(om.soln.lambda.upper, lam.upper, 14, [t 'om.soln.lambda.upper']);
0165 t_is(om.soln.lambda.mu_l, lam.mu_l, 14, [t 'om.soln.lambda.mu_l']);
0166 t_is(om.soln.lambda.mu_u, lam.mu_u, 14, [t 'om.soln.lambda.mu_u']);
0167
0168 t = 'om.get_soln(''var'', ''x'') : ';
0169 [x1, mu_l, mu_u] = om.get_soln('var', 'x');
0170 t_is(x1, x, 14, [t 'x']);
0171 t_is(mu_l, lam.lower, 14, [t 'mu_l']);
0172 t_is(mu_u, lam.upper, 14, [t 'mu_u']);
0173
0174 t = 'om.get_soln(''var'', ''mu_u'', ''x'') : ';
0175 t_is(om.get_soln('var', 'mu_u', 'x'), lam.upper, 14, [t 'mu_l']);
0176
0177 t = 'om.get_soln(''lin'', ''Ax'') : ';
0178 [g, mu_u] = om.get_soln('lin', 'Ax');
0179 t_is(g{1}, A*x-u, 14, [t 'A * x - u']);
0180 t_ok(all(isinf(g{2}) & g{2} < 0), [t 'l - A * x']);
0181 t_is(mu_u, lam.mu_u, 14, [t 'mu_u']);
0182
0183 t = 'om.get_soln(''lin'', {''mu_u'', ''mu_l'', ''g''}, ''Ax'') : ';
0184 [mu_u, mu_l, g] = om.get_soln('lin', {'mu_u', 'mu_l', 'g'}, 'Ax');
0185 t_is(g{1}, A*x-u, 14, [t 'A * x - u']);
0186 t_ok(all(isinf(g{2}) & g{2} < 0), [t 'l - A * x']);
0187 t_is(mu_l, lam.mu_l, 14, [t 'mu_l']);
0188 t_is(mu_u, lam.mu_u, 14, [t 'mu_u']);
0189
0190 t = 'om.get_soln(''qdc'', ''c'') : ';
0191 [f1, df, d2f] = om.get_soln('qdc', 'c');
0192 t_is(sum(f1), f, 14, [t 'f']);
0193 t_is(df, c, 14, [t 'df']);
0194 t_is(d2f, zeros(2,1), 14, [t 'd2f']);
0195
0196 t = 'om.get_soln(''qdc'', ''df'', ''c'') : ';
0197 df = om.get_soln('qdc', 'df', 'c');
0198 t_is(df, c, 14, [t 'df']);
0199
0200 t = 'parse_soln : ';
0201 t_is(om.soln.var.val.x, om.get_soln('var', 'x'), 14, [t 'var.val.x']);
0202 t_is(om.soln.var.mu_l.x, om.get_soln('var', 'mu_l', 'x'), 14, [t 'var.mu_l.x']);
0203 t_is(om.soln.var.mu_u.x, om.get_soln('var', 'mu_u', 'x'), 14, [t 'var.mu_u.x']);
0204 t_is(om.soln.lin.mu_l.Ax, om.get_soln('lin', 'mu_l', 'Ax'), 14, [t 'lin.mu_l.Ax']);
0205 t_is(om.soln.lin.mu_u.Ax, om.get_soln('lin', 'mu_u', 'Ax'), 14, [t 'lin.mu_u.Ax']);
0206 else
0207 t_skip(28, 'no MILP/MIQP solver installed');
0208 end
0209
0210 if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
0211 warning(s1.state, diff_alg_warn_id);
0212 end
0213
0214 t_end;
0215
0216 function TorF = have_miqp_solver()
0217 TorF = have_feature('cplex') || have_feature('glpk') || ...
0218 have_feature('gurobi') || have_feature('intlinprog') || ...
0219 have_feature('mosek');
0220