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
'CLP' : CLP
'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 % 'CLP' : CLP 0039 % 'CPLEX' : (or 500) CPLEX 0040 % 'GLPK' : GLPK, (LP problems only, i.e. empty H matrix) 0041 % 'GUROBI' : (or 700) Gurobi 0042 % 'IPOPT' : (or 400) IPOPT 0043 % 'MOSEK' : (or 600) MOSEK 0044 % 'OT' : (or 300) Optimization Toolbox, QUADPROG or LINPROG 0045 % verbose (0) - controls level of progress output displayed 0046 % 0 = no progress output 0047 % 1 = some progress output 0048 % 2 = verbose progress output 0049 % bp_opt - options vector for BP 0050 % cplex_opt - options struct for CPLEX 0051 % glpk_opt - options struct for GLPK 0052 % grb_opt - options struct for GBUROBI_MEX 0053 % ipopt_opt - options struct for IPOPT 0054 % linprog_opt - options struct for LINPROG 0055 % mips_opt - options struct for QPS_MIPS 0056 % mosek_opt - options struct for MOSEK 0057 % quadprog_opt - options struct for QUADPROG 0058 % PROBLEM : The inputs can alternatively be supplied in a single 0059 % PROBLEM struct with fields corresponding to the input arguments 0060 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0061 % 0062 % Outputs: 0063 % X : solution vector 0064 % F : final objective function value 0065 % EXITFLAG : exit flag 0066 % 1 = converged 0067 % 0 or negative values = algorithm specific failure codes 0068 % OUTPUT : output struct with the following fields: 0069 % alg - algorithm code of solver used 0070 % (others) - algorithm specific fields 0071 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0072 % multipliers on the constraints, with fields: 0073 % mu_l - lower (left-hand) limit on linear constraints 0074 % mu_u - upper (right-hand) limit on linear constraints 0075 % lower - lower bound on optimization variables 0076 % upper - upper bound on optimization variables 0077 % 0078 % Note the calling syntax is almost identical to that of QUADPROG 0079 % from MathWorks' Optimization Toolbox. The main difference is that 0080 % the linear constraints are specified with A, L, U instead of 0081 % A, B, Aeq, Beq. 0082 % 0083 % Calling syntax options: 0084 % [x, f, exitflag, output, lambda] = ... 0085 % qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt) 0086 % 0087 % x = qps_matpower(H, c, A, l, u) 0088 % x = qps_matpower(H, c, A, l, u, xmin, xmax) 0089 % x = qps_matpower(H, c, A, l, u, xmin, xmax, x0) 0090 % x = qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt) 0091 % x = qps_matpower(problem), where problem is a struct with fields: 0092 % H, c, A, l, u, xmin, xmax, x0, opt 0093 % all fields except 'c', 'A' and 'l' or 'u' are optional 0094 % x = qps_matpower(...) 0095 % [x, f] = qps_matpower(...) 0096 % [x, f, exitflag] = qps_matpower(...) 0097 % [x, f, exitflag, output] = qps_matpower(...) 0098 % [x, f, exitflag, output, lambda] = qps_matpower(...) 0099 % 0100 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0101 % H = [ 1003.1 4.3 6.3 5.9; 0102 % 4.3 2.2 2.1 3.9; 0103 % 6.3 2.1 3.5 4.8; 0104 % 5.9 3.9 4.8 10 ]; 0105 % c = zeros(4,1); 0106 % A = [ 1 1 1 1; 0107 % 0.17 0.11 0.10 0.18 ]; 0108 % l = [1; 0.10]; 0109 % u = [1; Inf]; 0110 % xmin = zeros(4,1); 0111 % x0 = [1; 0; 0; 1]; 0112 % opt = struct('verbose', 2); 0113 % [x, f, s, out, lambda] = qps_matpower(H, c, A, l, u, xmin, [], x0, opt); 0114 0115 % MATPOWER 0116 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0117 % by Ray Zimmerman, PSERC Cornell 0118 % 0119 % $Id: qps_matpower.m 2644 2015-03-11 19:34:22Z ray $ 0120 % 0121 % This file is part of MATPOWER. 0122 % Covered by the 3-clause BSD License (see LICENSE file for details). 0123 % See http://www.pserc.cornell.edu/matpower/ for more info. 0124 0125 %%----- input argument handling ----- 0126 %% gather inputs 0127 if nargin == 1 && isstruct(H) %% problem struct 0128 p = H; 0129 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0130 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0131 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0132 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0133 if isfield(p, 'u'), u = p.u; else, u = []; end 0134 if isfield(p, 'l'), l = p.l; else, l = []; end 0135 if isfield(p, 'A'), A = p.A; else, A = []; end 0136 if isfield(p, 'c'), c = p.c; else, c = []; end 0137 if isfield(p, 'H'), H = p.H; else, H = []; end 0138 else %% individual args 0139 if nargin < 9 0140 opt = []; 0141 if nargin < 8 0142 x0 = []; 0143 if nargin < 7 0144 xmax = []; 0145 if nargin < 6 0146 xmin = []; 0147 end 0148 end 0149 end 0150 end 0151 end 0152 0153 %% default options 0154 if ~isempty(opt) && isfield(opt, 'alg') && ~isempty(opt.alg) 0155 alg = opt.alg; 0156 %% convert integer codes to string values 0157 if ~ischar(alg) 0158 switch alg 0159 case 0 0160 alg = 'DEFAULT'; 0161 case 100 0162 alg = 'BPMPD'; 0163 case 200 0164 alg = 'MIPS'; 0165 opt.mips_opt.step_control = 0; 0166 case 250 0167 alg = 'MIPS'; 0168 opt.mips_opt.step_control = 1; 0169 case 300 0170 alg = 'OT'; 0171 case 400 0172 alg = 'IPOPT'; 0173 case 500 0174 alg = 'CPLEX'; 0175 case 600 0176 alg = 'MOSEK'; 0177 case 700 0178 alg = 'GUROBI'; 0179 otherwise 0180 error('qps_matpower: %d is not a valid algorithm code', alg); 0181 end 0182 end 0183 else 0184 alg = 'DEFAULT'; 0185 end 0186 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0187 verbose = opt.verbose; 0188 else 0189 verbose = 0; 0190 end 0191 if strcmp(alg, 'DEFAULT') 0192 if have_fcn('cplex') %% use CPLEX by default, if available 0193 alg = 'CPLEX'; 0194 elseif have_fcn('gurobi') %% if not, then Gurobi, if available 0195 alg = 'GUROBI'; 0196 elseif have_fcn('mosek') %% if not, then MOSEK, if available 0197 alg = 'MOSEK'; 0198 elseif have_fcn('bpmpd') %% if not, then BPMPD_MEX, if available 0199 alg = 'BPMPD'; 0200 elseif have_fcn('quadprog') %% if not, then Optimization Tbx, if available 0201 alg = 'OT'; 0202 elseif (isempty(H) || ~any(any(H))) && have_fcn('glpk') %% if not, and 0203 alg = 'GLPK'; %% prob is LP (not QP), then GLPK, if available 0204 else %% otherwise MIPS 0205 alg = 'MIPS'; 0206 end 0207 end 0208 0209 %%----- call the appropriate solver ----- 0210 switch alg 0211 case 'BPMPD' %% use BPMPD_MEX 0212 [x, f, eflag, output, lambda] = ... 0213 qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt); 0214 0215 if eflag == -99 0216 if verbose 0217 fprintf(' Retrying with QPS_MIPS solver ...\n\n'); 0218 end 0219 %% save (incorrect) solution from BPMPD 0220 bpmpd = struct('x', x, 'f', f, 'eflag', eflag, ... 0221 'output', output, 'lambda', lambda); 0222 opt.alg = 'MIPS'; 0223 [x, f, eflag, output, lambda] = ... 0224 qps_matpower(H, c, A, l, u, xmin, xmax, x0, opt); 0225 output.bpmpd = bpmpd; 0226 end 0227 case 'CLP' 0228 [x, f, eflag, output, lambda] = ... 0229 qps_clp(H, c, A, l, u, xmin, xmax, x0, opt); 0230 case 'CPLEX' 0231 [x, f, eflag, output, lambda] = ... 0232 qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt); 0233 case 'GLPK' 0234 [x, f, eflag, output, lambda] = ... 0235 qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt); 0236 case 'GUROBI' 0237 [x, f, eflag, output, lambda] = ... 0238 qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt); 0239 case 'IPOPT' 0240 [x, f, eflag, output, lambda] = ... 0241 qps_ipopt(H, c, A, l, u, xmin, xmax, x0, opt); 0242 case 'MIPS' 0243 %% set up options 0244 if ~isempty(opt) && isfield(opt, 'mips_opt') && ~isempty(opt.mips_opt) 0245 mips_opt = opt.mips_opt; 0246 else 0247 mips_opt = []; 0248 end 0249 mips_opt.verbose = verbose; 0250 0251 %% call solver 0252 [x, f, eflag, output, lambda] = ... 0253 qps_mips(H, c, A, l, u, xmin, xmax, x0, mips_opt); 0254 case 'MOSEK' 0255 [x, f, eflag, output, lambda] = ... 0256 qps_mosek(H, c, A, l, u, xmin, xmax, x0, 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 otherwise 0261 error('qps_matpower: %d is not a valid algorithm code', alg); 0262 end 0263 if ~isfield(output, 'alg') || isempty(output.alg) 0264 output.alg = alg; 0265 end