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 %   $Id: t_qps_matpower.m,v 1.12 2011/07/05 20:34:59 cvs Exp $
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %   Copyright (c) 2010-2011 by Power System Engineering Research Center (PSERC)
0008 %
0009 %   This file is part of MATPOWER.
0010 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0011 %
0012 %   MATPOWER is free software: you can redistribute it and/or modify
0013 %   it under the terms of the GNU General Public License as published
0014 %   by the Free Software Foundation, either version 3 of the License,
0015 %   or (at your option) any later version.
0016 %
0017 %   MATPOWER is distributed in the hope that it will be useful,
0018 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0019 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0020 %   GNU General Public License for more details.
0021 %
0022 %   You should have received a copy of the GNU General Public License
0023 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0024 %
0025 %   Additional permission under GNU GPL version 3 section 7
0026 %
0027 %   If you modify MATPOWER, or any covered work, to interface with
0028 %   other modules (such as MATLAB code and MEX-files) available in a
0029 %   MATLAB(R) or comparable environment containing parts covered
0030 %   under other licensing terms, the licensors of MATPOWER grant
0031 %   you additional permission to convey the resulting work.
0032 
0033 if nargin < 1
0034     quiet = 0;
0035 end
0036 
0037 algs = [100 200 250 400 300 500 600 700];
0038 names = {'BPMPD_MEX', 'MIPS', 'sc-MIPS', 'IPOPT', 'quadprog', 'CPLEX', 'MOSEK', 'Gurobi'};
0039 check = {'bpmpd', [], [], 'ipopt', 'quadprog', 'cplex', 'mosek', 'gurobi'};
0040 
0041 n = 36;
0042 t_begin(n*length(algs), quiet);
0043 
0044 for k = 1:length(algs)
0045     if ~isempty(check{k}) && ~have_fcn(check{k})
0046         t_skip(n, sprintf('%s not installed', names{k}));
0047     else
0048         opt = struct('verbose', 0, 'alg', algs(k));
0049         if strcmp(names{k}, 'MIPS') || strcmp(names{k}, 'sc-MIPS')
0050             opt.mips_opt.comptol = 1e-8;
0051         end
0052 %         if strcmp(names{k}, 'quadprog')
0053 %         end
0054         if strcmp(names{k}, 'CPLEX')
0055 %           alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
0056             alg = 2;        %% use dual simplex
0057             mpopt = mpoption('CPLEX_LPMETHOD', alg, 'CPLEX_QPMETHOD', min([4 alg]));
0058             opt.cplex_opt = cplex_options([], mpopt);
0059         end
0060         if strcmp(names{k}, 'MOSEK')
0061 %             alg = 5;        %% use dual simplex
0062             mpopt = mpoption;
0063 %             mpopt = mpoption(mpopt, 'MOSEK_LP_ALG', alg );
0064             mpopt = mpoption(mpopt, 'MOSEK_GAP_TOL', 1e-9);
0065             opt.mosek_opt = mosek_options([], mpopt);
0066         end
0067 
0068         t = sprintf('%s - 3-d LP : ', names{k});
0069         %% example from 'doc linprog'
0070         c = [-5; -4; -6];
0071         A = [1 -1  1;
0072              3  2  4;
0073              3  2  0];
0074         l = [];
0075         u = [20; 42; 30];
0076         xmin = [0; 0; 0];
0077         x0 = [];
0078         [x, f, s, out, lam] = qps_matpower([], c, A, l, u, xmin, [], [], opt);
0079         t_is(s, 1, 12, [t 'success']);
0080         t_is(x, [0; 15; 3], 6, [t 'x']);
0081         t_is(f, -78, 6, [t 'f']);
0082         t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0083         t_is(lam.mu_u, [0;1.5;0.5], 9, [t 'lam.mu_u']);
0084         t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
0085         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0086 
0087         t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
0088         %% from http://www.akiti.ca/QuadProgEx0Constr.html
0089         H = [5 -2 -1; -2 4 3; -1 3 5];
0090         c = [2; -35; -47];
0091         x0 = [0; 0; 0];
0092         [x, f, s, out, lam] = qps_matpower(H, c, [], [], [], [], [], [], opt);
0093         t_is(s, 1, 12, [t 'success']);
0094         t_is(x, [3; 5; 7], 8, [t 'x']);
0095         t_is(f, -249, 13, [t 'f']);
0096         t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
0097         t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
0098         t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
0099         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0100         
0101         t = sprintf('%s - constrained 2-d QP : ', names{k});
0102         %% example from 'doc quadprog'
0103         H = [   1   -1;
0104                 -1  2   ];
0105         c = [-2; -6];
0106         A = [   1   1;
0107                 -1  2;
0108                 2   1   ];
0109         l = [];
0110         u = [2; 2; 3];
0111         xmin = [0; 0];
0112         x0 = [];
0113         [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0114         t_is(s, 1, 12, [t 'success']);
0115         t_is(x, [2; 4]/3, 7, [t 'x']);
0116         t_is(f, -74/9, 6, [t 'f']);
0117         t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
0118         t_is(lam.mu_u, [28;4;0]/9, 7, [t 'lam.mu_u']);
0119         t_is(lam.lower, zeros(size(x)), 8, [t 'lam.lower']);
0120         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0121 
0122         t = sprintf('%s - constrained 4-d QP : ', names{k});
0123         %% from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm
0124         H = [   1003.1  4.3     6.3     5.9;
0125                 4.3     2.2     2.1     3.9;
0126                 6.3     2.1     3.5     4.8;
0127                 5.9     3.9     4.8     10  ];
0128         c = zeros(4,1);
0129         A = [   1       1       1       1;
0130                 0.17    0.11    0.10    0.18    ];
0131         l = [1; 0.10];
0132         u = [1; Inf];
0133         xmin = zeros(4,1);
0134         x0 = [1; 0; 0; 1];
0135         [x, f, s, out, lam] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);
0136         t_is(s, 1, 12, [t 'success']);
0137         t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0138         t_is(f, 3.29/3, 6, [t 'f']);
0139         t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0140         t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0141         t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0142         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0143 
0144         t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
0145         p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
0146         [x, f, s, out, lam] = qps_matpower(p);
0147         t_is(s, 1, 12, [t 'success']);
0148         t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
0149         t_is(f, 3.29/3, 6, [t 'f']);
0150         t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
0151         t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
0152         t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
0153         t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
0154 
0155         t = sprintf('%s - infeasible LP : ', names{k});
0156         p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
0157         [x, f, s, out, lam] = qps_matpower(p);
0158         t_ok(s <= 0, [t 'no success']);
0159     end
0160 end
0161 
0162 t_end;