QPS_MATPOWER Quadratic Program Solver for MATPOWER.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A common wrapper function for various QP solvers.
Solves the following QP (quadratic programming) problem:
min 1/2 X'*H*X + C'*X
X
subject to
L <= A*X <= U (linear constraints)
XMIN <= X <= XMAX (variable bounds)
Inputs (all optional except H, C, A and L):
H : matrix (possibly sparse) of quadratic cost coefficients
C : vector of linear cost coefficients
A, L, U : define the optional linear constraints. Default
values for the elements of L and U are -Inf and Inf,
respectively.
XMIN, XMAX : optional lower and upper bounds on the
X variables, defaults are -Inf and Inf, respectively.
X0 : optional starting value of optimization vector X
OPT : optional options structure with the following fields,
all of which are also optional (default values shown in
parentheses)
alg ('DEFAULT') : determines which solver to use, can be either
a (new-style) string or an (old-style) numerical alg code
'DEFAULT' : (or 0) automatic, first available of CPLEX,
Gurobi, MOSEK, BPMPD, Opt Tbx, GLPK (LPs only), MIPS
'MIPS' : (or 200) MIPS, MATLAB Interior Point Solver
pure MATLAB implementation of a primal-dual
interior point method, if mips_opt.step_control = 1
(or alg=250) it uses MIPS-sc, a step controlled
variant of MIPS
'BPMPD' : (or 100) BPMPD_MEX
'CPLEX' : (or 500) CPLEX
'GLPK' : GLPK, (LP problems only, i.e. empty H matrix)
'GUROBI' : (or 700) Gurobi
'IPOPT' : (or 400) IPOPT
'MOSEK' : (or 600) MOSEK
'OT' : (or 300) Optimization Toolbox, QUADPROG or LINPROG
verbose (0) - controls level of progress output displayed
0 = no progress output
1 = some progress output
2 = verbose progress output
bp_opt - options vector for BP
cplex_opt - options struct for CPLEX
glpk_opt - options struct for GLPK
grb_opt - options struct for GBUROBI_MEX
ipopt_opt - options struct for IPOPT
linprog_opt - options struct for LINPROG
mips_opt - options struct for QPS_MIPS
mosek_opt - options struct for MOSEK
quadprog_opt - options struct for QUADPROG
PROBLEM : The inputs can alternatively be supplied in a single
PROBLEM struct with fields corresponding to the input arguments
described above: H, c, A, l, u, xmin, xmax, x0, opt
Outputs:
X : solution vector
F : final objective function value
EXITFLAG : exit flag
1 = converged
0 or negative values = algorithm specific failure codes
OUTPUT : output struct with the following fields:
alg - algorithm code of solver used
(others) - algorithm specific fields
LAMBDA : struct containing the Langrange and Kuhn-Tucker
multipliers on the constraints, with fields:
mu_l - lower (left-hand) limit on linear constraints
mu_u - upper (right-hand) limit on linear constraints
lower - lower bound on optimization variables
upper - upper bound on optimization variables
Note the calling syntax is almost identical to that of QUADPROG
from MathWorks' Optimization Toolbox. The main difference is that
the linear constraints are specified with A, L, U instead of
A, B, Aeq, Beq.
Calling syntax options:
[x, f, exitflag, output, lambda] = ...
qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_matpower(H, c, A, l, u)
x = qps_matpower(H, c, A, l, u, xmin, xmax)
x = qps_matpower(H, c, A, l, u, xmin, xmax, x0)
x = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_matpower(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = qps_matpower(...)
[x, f] = qps_matpower(...)
[x, f, exitflag] = qps_matpower(...)
[x, f, exitflag, output] = qps_matpower(...)
[x, f, exitflag, output, lambda] = qps_matpower(...)
Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
H = [ 1003.1 4.3 6.3 5.9;
4.3 2.2 2.1 3.9;
6.3 2.1 3.5 4.8;
5.9 3.9 4.8 10 ];
c = zeros(4,1);
A = [ 1 1 1 1;
0.17 0.11 0.10 0.18 ];
l = [1; 0.10];
u = [1; Inf];
xmin = zeros(4,1);
x0 = [1; 0; 0; 1];
opt = struct('verbose', 2);
[x, f, s, out, lambda] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt);0001 function [x, f, eflag, output, lambda] = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_MATPOWER Quadratic Program Solver for MATPOWER. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_MATPOWER(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A common wrapper function for various QP solvers. 0006 % Solves the following QP (quadratic programming) problem: 0007 % 0008 % min 1/2 X'*H*X + C'*X 0009 % X 0010 % 0011 % subject to 0012 % 0013 % L <= A*X <= U (linear constraints) 0014 % XMIN <= X <= XMAX (variable bounds) 0015 % 0016 % Inputs (all optional except H, C, A and L): 0017 % H : matrix (possibly sparse) of quadratic cost coefficients 0018 % C : vector of linear cost coefficients 0019 % A, L, U : define the optional linear constraints. Default 0020 % values for the elements of L and U are -Inf and Inf, 0021 % respectively. 0022 % XMIN, XMAX : optional lower and upper bounds on the 0023 % X variables, defaults are -Inf and Inf, respectively. 0024 % X0 : optional starting value of optimization vector X 0025 % OPT : optional options structure with the following fields, 0026 % all of which are also optional (default values shown in 0027 % parentheses) 0028 % alg ('DEFAULT') : determines which solver to use, can be either 0029 % a (new-style) string or an (old-style) numerical alg code 0030 % 'DEFAULT' : (or 0) automatic, first available of CPLEX, 0031 % Gurobi, MOSEK, BPMPD, Opt Tbx, GLPK (LPs only), MIPS 0032 % 'MIPS' : (or 200) MIPS, MATLAB Interior Point Solver 0033 % pure MATLAB implementation of a primal-dual 0034 % interior point method, if mips_opt.step_control = 1 0035 % (or alg=250) it uses MIPS-sc, a step controlled 0036 % variant of MIPS 0037 % 'BPMPD' : (or 100) BPMPD_MEX 0038 % 'CPLEX' : (or 500) CPLEX 0039 % 'GLPK' : GLPK, (LP problems only, i.e. empty H matrix) 0040 % 'GUROBI' : (or 700) Gurobi 0041 % 'IPOPT' : (or 400) IPOPT 0042 % 'MOSEK' : (or 600) MOSEK 0043 % 'OT' : (or 300) Optimization Toolbox, QUADPROG or LINPROG 0044 % verbose (0) - controls level of progress output displayed 0045 % 0 = no progress output 0046 % 1 = some progress output 0047 % 2 = verbose progress output 0048 % bp_opt - options vector for BP 0049 % cplex_opt - options struct for CPLEX 0050 % glpk_opt - options struct for GLPK 0051 % grb_opt - options struct for GBUROBI_MEX 0052 % ipopt_opt - options struct for IPOPT 0053 % linprog_opt - options struct for LINPROG 0054 % mips_opt - options struct for QPS_MIPS 0055 % mosek_opt - options struct for MOSEK 0056 % quadprog_opt - options struct for QUADPROG 0057 % PROBLEM : The inputs can alternatively be supplied in a single 0058 % PROBLEM struct with fields corresponding to the input arguments 0059 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0060 % 0061 % Outputs: 0062 % X : solution vector 0063 % F : final objective function value 0064 % EXITFLAG : exit flag 0065 % 1 = converged 0066 % 0 or negative values = algorithm specific failure codes 0067 % OUTPUT : output struct with the following fields: 0068 % alg - algorithm code of solver used 0069 % (others) - algorithm specific fields 0070 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0071 % multipliers on the constraints, with fields: 0072 % mu_l - lower (left-hand) limit on linear constraints 0073 % mu_u - upper (right-hand) limit on linear constraints 0074 % lower - lower bound on optimization variables 0075 % upper - upper bound on optimization variables 0076 % 0077 % Note the calling syntax is almost identical to that of QUADPROG 0078 % from MathWorks' Optimization Toolbox. The main difference is that 0079 % the linear constraints are specified with A, L, U instead of 0080 % A, B, Aeq, Beq. 0081 % 0082 % Calling syntax options: 0083 % [x, f, exitflag, output, lambda] = ... 0084 % qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt) 0085 % 0086 % x = qps_matpower(H, c, A, l, u) 0087 % x = qps_matpower(H, c, A, l, u, xmin, xmax) 0088 % x = qps_matpower(H, c, A, l, u, xmin, xmax, x0) 0089 % x = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt) 0090 % x = qps_matpower(problem), where problem is a struct with fields: 0091 % H, c, A, l, u, xmin, xmax, x0, opt 0092 % all fields except 'c', 'A' and 'l' or 'u' are optional 0093 % x = qps_matpower(...) 0094 % [x, f] = qps_matpower(...) 0095 % [x, f, exitflag] = qps_matpower(...) 0096 % [x, f, exitflag, output] = qps_matpower(...) 0097 % [x, f, exitflag, output, lambda] = qps_matpower(...) 0098 % 0099 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0100 % H = [ 1003.1 4.3 6.3 5.9; 0101 % 4.3 2.2 2.1 3.9; 0102 % 6.3 2.1 3.5 4.8; 0103 % 5.9 3.9 4.8 10 ]; 0104 % c = zeros(4,1); 0105 % A = [ 1 1 1 1; 0106 % 0.17 0.11 0.10 0.18 ]; 0107 % l = [1; 0.10]; 0108 % u = [1; Inf]; 0109 % xmin = zeros(4,1); 0110 % x0 = [1; 0; 0; 1]; 0111 % opt = struct('verbose', 2); 0112 % [x, f, s, out, lambda] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt); 0113 0114 % MATPOWER 0115 % $Id: qps_matpower.m 2385 2014-10-11 15:33:45Z ray $ 0116 % by Ray Zimmerman, PSERC Cornell 0117 % Copyright (c) 2010-2014 by Power System Engineering Research Center (PSERC) 0118 % 0119 % This file is part of MATPOWER. 0120 % See http://www.pserc.cornell.edu/matpower/ for more info. 0121 % 0122 % MATPOWER is free software: you can redistribute it and/or modify 0123 % it under the terms of the GNU General Public License as published 0124 % by the Free Software Foundation, either version 3 of the License, 0125 % or (at your option) any later version. 0126 % 0127 % MATPOWER is distributed in the hope that it will be useful, 0128 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0129 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0130 % GNU General Public License for more details. 0131 % 0132 % You should have received a copy of the GNU General Public License 0133 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0134 % 0135 % Additional permission under GNU GPL version 3 section 7 0136 % 0137 % If you modify MATPOWER, or any covered work, to interface with 0138 % other modules (such as MATLAB code and MEX-files) available in a 0139 % MATLAB(R) or comparable environment containing parts covered 0140 % under other licensing terms, the licensors of MATPOWER grant 0141 % you additional permission to convey the resulting work. 0142 0143 %%----- input argument handling ----- 0144 %% gather inputs 0145 if nargin == 1 && isstruct(H) %% problem struct 0146 p = H; 0147 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0148 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0149 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0150 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0151 if isfield(p, 'u'), u = p.u; else, u = []; end 0152 if isfield(p, 'l'), l = p.l; else, l = []; end 0153 if isfield(p, 'A'), A = p.A; else, A = []; end 0154 if isfield(p, 'c'), c = p.c; else, c = []; end 0155 if isfield(p, 'H'), H = p.H; else, H = []; end 0156 else %% individual args 0157 if nargin < 9 0158 opt = []; 0159 if nargin < 8 0160 x0 = []; 0161 if nargin < 7 0162 xmax = []; 0163 if nargin < 6 0164 xmin = []; 0165 end 0166 end 0167 end 0168 end 0169 end 0170 0171 %% default options 0172 if ~isempty(opt) && isfield(opt, 'alg') && ~isempty(opt.alg) 0173 alg = opt.alg; 0174 %% convert integer codes to string values 0175 if ~ischar(alg) 0176 switch alg 0177 case 0 0178 alg = 'DEFAULT'; 0179 case 100 0180 alg = 'BPMPD'; 0181 case 200 0182 alg = 'MIPS'; 0183 opt.mips_opt.step_control = 0; 0184 case 250 0185 alg = 'MIPS'; 0186 opt.mips_opt.step_control = 1; 0187 case 300 0188 alg = 'OT'; 0189 case 400 0190 alg = 'IPOPT'; 0191 case 500 0192 alg = 'CPLEX'; 0193 case 600 0194 alg = 'MOSEK'; 0195 case 700 0196 alg = 'GUROBI'; 0197 otherwise 0198 error('qps_matpower: %d is not a valid algorithm code', alg); 0199 end 0200 end 0201 else 0202 alg = 'DEFAULT'; 0203 end 0204 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0205 verbose = opt.verbose; 0206 else 0207 verbose = 0; 0208 end 0209 if strcmp(alg, 'DEFAULT') 0210 if have_fcn('cplex') %% use CPLEX by default, if available 0211 alg = 'CPLEX'; 0212 elseif have_fcn('gurobi') %% if not, then Gurobi, if available 0213 alg = 'GUROBI'; 0214 elseif have_fcn('mosek') %% if not, then MOSEK, if available 0215 alg = 'MOSEK'; 0216 elseif have_fcn('bpmpd') %% if not, then BPMPD_MEX, if available 0217 alg = 'BPMPD'; 0218 elseif have_fcn('quadprog') %% if not, then Optimization Tbx, if available 0219 alg = 'OT'; 0220 elseif (isempty(H) || ~any(any(H))) && have_fcn('glpk') %% if not, and 0221 alg = 'GLPK'; %% prob is LP (not QP), then GLPK, if available 0222 else %% otherwise MIPS 0223 alg = 'MIPS'; 0224 end 0225 end 0226 0227 %%----- call the appropriate solver ----- 0228 switch alg 0229 case 'BPMPD' %% use BPMPD_MEX 0230 [x, f, eflag, output, lambda] = ... 0231 qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt); 0232 0233 if eflag == -99 0234 if verbose 0235 fprintf(' Retrying with QPS_MIPS solver ...\n\n'); 0236 end 0237 %% save (incorrect) solution from BPMPD 0238 bpmpd = struct('x', x, 'f', f, 'eflag', eflag, ... 0239 'output', output, 'lambda', lambda); 0240 opt.alg = 'MIPS'; 0241 [x, f, eflag, output, lambda] = ... 0242 qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt); 0243 output.bpmpd = bpmpd; 0244 end 0245 case 'MIPS' 0246 %% set up options 0247 if ~isempty(opt) && isfield(opt, 'mips_opt') && ~isempty(opt.mips_opt) 0248 mips_opt = opt.mips_opt; 0249 else 0250 mips_opt = []; 0251 end 0252 mips_opt.verbose = verbose; 0253 0254 %% call solver 0255 [x, f, eflag, output, lambda] = ... 0256 qps_mips(H, c, A, l, u, xmin, xmax, x0, mips_opt); 0257 case 'OT' %% use QUADPROG or LINPROG from Opt Tbx ver 2.x+ 0258 [x, f, eflag, output, lambda] = ... 0259 qps_ot(H, c, A, l, u, xmin, xmax, x0, opt); 0260 case 'IPOPT' 0261 [x, f, eflag, output, lambda] = ... 0262 qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt); 0263 case 'CPLEX' 0264 [x, f, eflag, output, lambda] = ... 0265 qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt); 0266 case 'MOSEK' 0267 [x, f, eflag, output, lambda] = ... 0268 qps_mosek(H, c, A, l, u, xmin, xmax, x0, opt); 0269 case 'GUROBI' 0270 [x, f, eflag, output, lambda] = ... 0271 qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt); 0272 case 'CLP' 0273 [x, f, eflag, output, lambda] = ... 0274 qps_clp(H, c, A, l, u, xmin, xmax, x0, opt); 0275 case 'GLPK' 0276 [x, f, eflag, output, lambda] = ... 0277 qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt); 0278 otherwise 0279 error('qps_matpower: %d is not a valid algorithm code', alg); 0280 end 0281 if ~isfield(output, 'alg') || isempty(output.alg) 0282 output.alg = alg; 0283 end