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