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