t_qps_matpower

T_QPS_MATPOWER  Tests of QPS_MATPOWER QP solvers.

0001 function t_qps_matpower(quiet)
0002 %T_QPS_MATPOWER  Tests of QPS_MATPOWER QP solvers.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   $Id: t_qps_matpower.m 2644 2015-03-11 19:34:22Z ray $
0009 %
0010 %   This file is part of MATPOWER.
0011 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0012 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0013 
0014 if nargin < 1
0015     quiet = 0;
0016 end
0017 
0018 algs = {'BPMPD', 'MIPS', 250, 'IPOPT', 'OT', 'CPLEX', 'MOSEK', 'GUROBI', 'CLP', 'GLPK'};
0019 names = {'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'linprog/quadprog', 'CPLEX', 'MOSEK', 'Gurobi', 'CLP', 'glpk'};
0020 check = {'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi', 'clp', 'glpk'};
0021 does_qp = [1 1 1 1 1 1 1 1 1 0];
0022 
0023 n = 36;
0024 nqp = 28;
0025 t_begin(n*length(algs), quiet);
0026 
0027 for k = 1:length(algs)
0028     if ~isempty(check{k}) && ~have_fcn(check{k})
0029         t_skip(n, sprintf('%s not installed', names{k}));
0030     else
0031         opt = struct('verbose', 0, 'alg', algs{k});
0032         if strcmp(names{k}, 'MIPS') || strcmp(names{k}, 'sc-MIPS')
0033             opt.mips_opt.comptol = 1e-8;
0034         end
0035 %         if strcmp(names{k}, 'linprog/quadprog')
0036 %             opt.verbose = 2;
0037 %             opt.linprog_opt.Algorithm = 'interior-point';
0038 %             opt.linprog_opt.Algorithm = 'active-set';
0039 %             opt.linprog_opt.Algorithm = 'simplex';
0040 %             opt.linprog_opt.Algorithm = 'dual-simplex';
0041 %         end
0042         if strcmp(names{k}, 'CPLEX')
0043 %           alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
0044             alg = 2;        %% use dual simplex
0045             mpopt = mpoption('cplex.lpmethod', alg, 'cplex.qpmethod', min([4 alg]));
0046             opt.cplex_opt = cplex_options([], mpopt);
0047         end
0048         if strcmp(names{k}, 'MOSEK')
0049             mpopt = mpoption;
0050 %             sc = mosek_symbcon;
0051 %             alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
0052 %             mpopt = mpoption(mpopt, 'mosek.lp_alg', alg );
0053             mpopt = mpoption(mpopt, 'mosek.gap_tol', 1e-10);
0054             opt.mosek_opt = mosek_options([], mpopt);
0055         end
0056 
0057         t = sprintf('%s - 3-d LP : ', names{k});
0058         %% based on example from 'doc linprog'
0059         c = [-5; -4; -6];
0060         A = [1 -1  1;
0061              -3  -2  -4;
0062              3  2  0];
0063         l = [-Inf; -42; -Inf];
0064         u = [20; Inf; 30];
0065         xmin = [0; 0; 0];
0066         x0 = [];
0067         [x, f, s, out, lam] = qps_matpower([], c, A, l, u, xmin, [], [], opt);
0068         t_is(s, 1, 12, [t 'success']);
0069         t_is(x, [0; 15; 3], 6, [t 'x']);
0070         t_is(f, -78, 6, [t 'f']);
0071         t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
0072         t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
0073         if strcmp(algs{k}, 'CLP') && ~have_fcn('opti_clp')
0074             t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0075         else
0076             t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0077             t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);
0078         end
0079 
0080         if does_qp(k)
0081             t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0082             %% from http://www.akiti.ca/QuadProgEx0Constr.html
0083             H = [5 -2 -1; -2 4 3; -1 3 5];
0084             c = [2; -35; -47];
0085             x0 = [0; 0; 0];
0086             [x, f, s, out, lam] = qps_matpower(H, c, [], [], [], [], [], [], opt);
0087             t_is(s, 1, 12, [t 'success']);
0088             t_is(x, [3; 5; 7], 8, [t 'x']);
0089             t_is(f, -249, 13, [t 'f']);
0090             t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0091             t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0092             t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0093             t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0094         
0095             t = sprintf('%s - constrained 2-d QP : ', names{k});
0096             %% example from 'doc quadprog'
0097             H = [   1   -1;
0098                     -1  2   ];
0099             c = [-2; -6];
0100             A = [   1   1;
0101                     -1  2;
0102                     2   1   ];
0103             l = [];
0104             u = [2; 2; 3];
0105             xmin = [0; 0];
0106             x0 = [];
0107             [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0108             t_is(s, 1, 12, [t 'success']);
0109             t_is(x, [2; 4]/3, 7, [t 'x']);
0110             t_is(f, -74/9, 6, [t 'f']);
0111             t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0112             t_is(lam.mu_u, [28;4;0]/9, 7, [t 'lam.mu_u']);
0113             if strcmp(algs{k}, 'CLP') && ~have_fcn('opti_clp')
0114                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0115             else
0116                 t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
0117                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0118             end
0119 
0120             t = sprintf('%s - constrained 4-d QP : ', names{k});
0121             %% from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm
0122             H = [   1003.1  4.3     6.3     5.9;
0123                     4.3     2.2     2.1     3.9;
0124                     6.3     2.1     3.5     4.8;
0125                     5.9     3.9     4.8     10  ];
0126             c = zeros(4,1);
0127             A = [   1       1       1       1;
0128                     0.17    0.11    0.10    0.18    ];
0129             l = [1; 0.10];
0130             u = [1; Inf];
0131             xmin = zeros(4,1);
0132             x0 = [1; 0; 0; 1];
0133             [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0134             t_is(s, 1, 12, [t 'success']);
0135             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0136             t_is(f, 3.29/3, 6, [t 'f']);
0137             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0138             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0139             if strcmp(algs{k}, 'CLP') && ~have_fcn('opti_clp')
0140                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0141             else
0142                 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0143                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0144             end
0145 
0146             t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0147             p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0148             [x, f, s, out, lam] = qps_matpower(p);
0149             t_is(s, 1, 12, [t 'success']);
0150             t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0151             t_is(f, 3.29/3, 6, [t 'f']);
0152             t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0153             t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0154             if strcmp(algs{k}, 'CLP') && ~have_fcn('opti_clp')
0155                 t_skip(2, [t 'lam.lower/upper : MEXCLP does not return multipliers on var bounds']);
0156             else
0157                 t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0158                 t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0159             end
0160         else
0161             t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
0162         end
0163 
0164         t = sprintf('%s - infeasible LP : ', names{k});
0165         p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0166         [x, f, s, out, lam] = qps_matpower(p);
0167         t_ok(s <= 0, [t 'no success']);
0168     end
0169 end
0170 
0171 t_end;