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