Home > matpower7.1 > mp-opt-model > lib > t > t_om_solve_miqps.m

t_om_solve_miqps

PURPOSE ^

T_OM_SOLVE_MIQPS Tests of MIQP solvers via OM.SOLVE().

SYNOPSIS ^

function t_om_solve_miqps(quiet)

DESCRIPTION ^

T_OM_SOLVE_MIQPS  Tests of MIQP solvers via OM.SOLVE().

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function t_om_solve_miqps(quiet)
0002 %T_OM_SOLVE_MIQPS  Tests of MIQP solvers via OM.SOLVE().
0003 
0004 %   MP-Opt-Model
0005 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MP-Opt-Model.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://github.com/MATPOWER/mp-opt-model for more info.
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 %           alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
0059             alg = 2;        %% use dual simplex
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 %             sc = mosek_symbcon;
0066 %             alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
0067 %             alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
0068 %             mpopt.mosek.lp_alg = alg;
0069             mpopt.mosek.gap_tol = 1e-10;
0070 %             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-10;
0071 %             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS = 1e-10;
0072 %             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS = 1e-10;
0073 %             mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-10;
0074             vnum = have_feature('mosek', 'vnum');
0075             if vnum >= 8
0076 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1e-10;
0077 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1e-10;
0078 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS = 1e-10;
0079 %                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED = 1e-10;
0080                 mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
0081             end
0082 %             opt.verbose = 3;
0083             opt.mosek_opt = mosek_options([], mpopt);
0084         end
0085 
0086 % opt.verbose = 3;
0087         t = sprintf('%s - 2-d MILP : ', names{k});
0088         %% from MOSEK 6.0 Guided Tour, section  7.13.1, https://docs.mosek.com/6.0/toolbox/node009.html
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             %% from cplexmiqpex.m, CPLEX_Studio_Academic124/cplex/examples/src/matlab/cplexmiqpex.m
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 % opt.verbose = 0;
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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005