QPS_CPLEX Quadratic Program Solver based on CPLEX.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A wrapper function providing a MATPOWER standardized interface for using
CPLEXQP or CPLEXLP to solve 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)
verbose (0) - controls level of progress output displayed
0 = no progress output
1 = some progress output
2 = verbose progress output
cplex_opt - options struct for CPLEX, value in
verbose overrides these options
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 : CPLEXQP/CPLEXLP exit flag
(see CPLEXQP and CPLEXLP documentation for details)
OUTPUT : CPLEXQP/CPLEXLP output struct
(see CPLEXQP and CPLEXLP documentation for details)
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_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_cplex(H, c, A, l, u)
x = qps_cplex(H, c, A, l, u, xmin, xmax)
x = qps_cplex(H, c, A, l, u, xmin, xmax, x0)
x = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_cplex(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_cplex(...)
[x, f] = qps_cplex(...)
[x, f, exitflag] = qps_cplex(...)
[x, f, exitflag, output] = qps_cplex(...)
[x, f, exitflag, output, lambda] = qps_cplex(...)
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_cplex(H, c, A, l, u, xmin, [], x0, opt);
See also CPLEXQP, CPLEXLP, CPLEX_OPTIONS.0001 function [x, f, eflag, output, lambda] = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_CPLEX Quadratic Program Solver based on CPLEX. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A wrapper function providing a MATPOWER standardized interface for using 0006 % CPLEXQP or CPLEXLP to solve the following QP (quadratic programming) 0007 % problem: 0008 % 0009 % min 1/2 X'*H*X + C'*X 0010 % X 0011 % 0012 % subject to 0013 % 0014 % L <= A*X <= U (linear constraints) 0015 % XMIN <= X <= XMAX (variable bounds) 0016 % 0017 % Inputs (all optional except H, C, A and L): 0018 % H : matrix (possibly sparse) of quadratic cost coefficients 0019 % C : vector of linear cost coefficients 0020 % A, L, U : define the optional linear constraints. Default 0021 % values for the elements of L and U are -Inf and Inf, 0022 % respectively. 0023 % XMIN, XMAX : optional lower and upper bounds on the 0024 % X variables, defaults are -Inf and Inf, respectively. 0025 % X0 : optional starting value of optimization vector X 0026 % OPT : optional options structure with the following fields, 0027 % all of which are also optional (default values shown in 0028 % parentheses) 0029 % verbose (0) - controls level of progress output displayed 0030 % 0 = no progress output 0031 % 1 = some progress output 0032 % 2 = verbose progress output 0033 % cplex_opt - options struct for CPLEX, value in 0034 % verbose overrides these options 0035 % PROBLEM : The inputs can alternatively be supplied in a single 0036 % PROBLEM struct with fields corresponding to the input arguments 0037 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0038 % 0039 % Outputs: 0040 % X : solution vector 0041 % F : final objective function value 0042 % EXITFLAG : CPLEXQP/CPLEXLP exit flag 0043 % (see CPLEXQP and CPLEXLP documentation for details) 0044 % OUTPUT : CPLEXQP/CPLEXLP output struct 0045 % (see CPLEXQP and CPLEXLP documentation for details) 0046 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0047 % multipliers on the constraints, with fields: 0048 % mu_l - lower (left-hand) limit on linear constraints 0049 % mu_u - upper (right-hand) limit on linear constraints 0050 % lower - lower bound on optimization variables 0051 % upper - upper bound on optimization variables 0052 % 0053 % Note the calling syntax is almost identical to that of QUADPROG 0054 % from MathWorks' Optimization Toolbox. The main difference is that 0055 % the linear constraints are specified with A, L, U instead of 0056 % A, B, Aeq, Beq. 0057 % 0058 % Calling syntax options: 0059 % [x, f, exitflag, output, lambda] = ... 0060 % qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 0061 % 0062 % x = qps_cplex(H, c, A, l, u) 0063 % x = qps_cplex(H, c, A, l, u, xmin, xmax) 0064 % x = qps_cplex(H, c, A, l, u, xmin, xmax, x0) 0065 % x = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 0066 % x = qps_cplex(problem), where problem is a struct with fields: 0067 % H, c, A, l, u, xmin, xmax, x0, opt 0068 % all fields except 'c', 'A' and 'l' or 'u' are optional 0069 % x = qps_cplex(...) 0070 % [x, f] = qps_cplex(...) 0071 % [x, f, exitflag] = qps_cplex(...) 0072 % [x, f, exitflag, output] = qps_cplex(...) 0073 % [x, f, exitflag, output, lambda] = qps_cplex(...) 0074 % 0075 % 0076 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0077 % H = [ 1003.1 4.3 6.3 5.9; 0078 % 4.3 2.2 2.1 3.9; 0079 % 6.3 2.1 3.5 4.8; 0080 % 5.9 3.9 4.8 10 ]; 0081 % c = zeros(4,1); 0082 % A = [ 1 1 1 1; 0083 % 0.17 0.11 0.10 0.18 ]; 0084 % l = [1; 0.10]; 0085 % u = [1; Inf]; 0086 % xmin = zeros(4,1); 0087 % x0 = [1; 0; 0; 1]; 0088 % opt = struct('verbose', 2); 0089 % [x, f, s, out, lam] = qps_cplex(H, c, A, l, u, xmin, [], x0, opt); 0090 % 0091 % See also CPLEXQP, CPLEXLP, CPLEX_OPTIONS. 0092 0093 % MATPOWER 0094 % $Id: qps_cplex.m,v 1.4 2010/12/16 20:59:57 cvs Exp $ 0095 % by Ray Zimmerman, PSERC Cornell 0096 % Copyright (c) 2010 by Power System Engineering Research Center (PSERC) 0097 % 0098 % This file is part of MATPOWER. 0099 % See http://www.pserc.cornell.edu/matpower/ for more info. 0100 % 0101 % MATPOWER is free software: you can redistribute it and/or modify 0102 % it under the terms of the GNU General Public License as published 0103 % by the Free Software Foundation, either version 3 of the License, 0104 % or (at your option) any later version. 0105 % 0106 % MATPOWER is distributed in the hope that it will be useful, 0107 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0108 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0109 % GNU General Public License for more details. 0110 % 0111 % You should have received a copy of the GNU General Public License 0112 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0113 % 0114 % Additional permission under GNU GPL version 3 section 7 0115 % 0116 % If you modify MATPOWER, or any covered work, to interface with 0117 % other modules (such as MATLAB code and MEX-files) available in a 0118 % MATLAB(R) or comparable environment containing parts covered 0119 % under other licensing terms, the licensors of MATPOWER grant 0120 % you additional permission to convey the resulting work. 0121 0122 %% check for CPLEX 0123 % if ~have_fcn('cplexqp') 0124 % error('qps_cplex: requires the Matlab interface for CPLEX'); 0125 % end 0126 0127 %%----- input argument handling ----- 0128 %% gather inputs 0129 if nargin == 1 && isstruct(H) %% problem struct 0130 p = H; 0131 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0132 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0133 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0134 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0135 if isfield(p, 'u'), u = p.u; else, u = []; end 0136 if isfield(p, 'l'), l = p.l; else, l = []; end 0137 if isfield(p, 'A'), A = p.A; else, A = []; end 0138 if isfield(p, 'c'), c = p.c; else, c = []; end 0139 if isfield(p, 'H'), H = p.H; else, H = []; end 0140 else %% individual args 0141 if nargin < 9 0142 opt = []; 0143 if nargin < 8 0144 x0 = []; 0145 if nargin < 7 0146 xmax = []; 0147 if nargin < 6 0148 xmin = []; 0149 end 0150 end 0151 end 0152 end 0153 end 0154 0155 %% define nx, set default values for missing optional inputs 0156 if isempty(H) || ~any(any(H)) 0157 if isempty(A) && isempty(xmin) && isempty(xmax) 0158 error('qps_cplex: LP problem must include constraints or variable bounds'); 0159 else 0160 if ~isempty(A) 0161 nx = size(A, 2); 0162 elseif ~isempty(xmin) 0163 nx = length(xmin); 0164 else % if ~isempty(xmax) 0165 nx = length(xmax); 0166 end 0167 end 0168 else 0169 nx = size(H, 1); 0170 end 0171 if isempty(c) 0172 c = zeros(nx, 1); 0173 end 0174 if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0175 (isempty(u) || all(u == Inf)) 0176 A = sparse(0,nx); %% no limits => no linear constraints 0177 end 0178 nA = size(A, 1); %% number of original linear constraints 0179 if isempty(u) %% By default, linear inequalities are ... 0180 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0181 end 0182 if isempty(l) 0183 l = -Inf * ones(nA, 1); %% ... unbounded below. 0184 end 0185 if isempty(xmin) %% By default, optimization variables are ... 0186 xmin = -Inf * ones(nx, 1); %% ... unbounded below and ... 0187 end 0188 if isempty(xmax) 0189 xmax = Inf * ones(nx, 1); %% ... unbounded above. 0190 end 0191 if isempty(x0) 0192 x0 = zeros(nx, 1); 0193 end 0194 0195 %% default options 0196 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0197 verbose = opt.verbose; 0198 else 0199 verbose = 0; 0200 end 0201 % if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it) 0202 % max_it = opt.max_it; 0203 % else 0204 % max_it = 0; 0205 % end 0206 0207 %% split up linear constraints 0208 ieq = find( abs(u-l) <= eps ); %% equality 0209 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0210 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0211 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0212 Ae = A(ieq, :); 0213 be = u(ieq); 0214 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0215 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0216 0217 %% grab some dimensions 0218 nlt = length(ilt); %% number of upper bounded linear inequalities 0219 ngt = length(igt); %% number of lower bounded linear inequalities 0220 nbx = length(ibx); %% number of doubly bounded linear inequalities 0221 0222 %% set up options struct for CPLEX 0223 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0224 cplex_opt = cplex_options(opt.cplex_opt); 0225 else 0226 cplex_opt = cplex_options; 0227 end 0228 0229 vrb = max([0 verbose-1]); 0230 cplex_opt.barrier.display = vrb; 0231 cplex_opt.conflict.display = vrb; 0232 cplex_opt.mip.display = vrb; 0233 cplex_opt.sifting.display = vrb; 0234 cplex_opt.simplex.display = vrb; 0235 cplex_opt.tune.display = vrb; 0236 % if max_it 0237 % cplex_opt. %% not sure what to set here 0238 % end 0239 0240 if isempty(Ai) && isempty(Ae) 0241 unconstrained = 1; 0242 Ae = sparse(1, nx); 0243 be = 0; 0244 else 0245 unconstrained = 0; 0246 end 0247 0248 %% call the solver 0249 if verbose 0250 cplex = Cplex('null'); 0251 methods = { 0252 'default', 0253 'primal simplex', 0254 'dual simplex', 0255 'network simplex', 0256 'barrier', 0257 'sifting', 0258 'concurrent' 0259 }; 0260 end 0261 if isempty(H) || ~any(any(H)) 0262 if verbose 0263 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0264 cplex.getVersion, methods{cplex_opt.lpmethod+1}); 0265 end 0266 [x, f, eflag, output, lam] = ... 0267 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0268 else 0269 if verbose 0270 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0271 cplex.getVersion, methods{cplex_opt.qpmethod+1}); 0272 end 0273 [x, f, eflag, output, lam] = ... 0274 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0275 end 0276 0277 %% check for empty results (in case optimization failed) 0278 if isempty(x) 0279 x = NaN(nx, 1); 0280 end 0281 if isempty(f) 0282 f = NaN; 0283 end 0284 if isempty(lam) 0285 lam.ineqlin = NaN(length(bi), 1); 0286 lam.eqlin = NaN(length(be), 1); 0287 lam.lower = NaN(nx, 1); 0288 lam.upper = NaN(nx, 1); 0289 mu_l = NaN(nA, 1); 0290 mu_u = NaN(nA, 1); 0291 else 0292 mu_l = zeros(nA, 1); 0293 mu_u = zeros(nA, 1); 0294 end 0295 if unconstrained 0296 lam.eqlin = []; 0297 end 0298 0299 %% repackage lambdas 0300 kl = find(lam.eqlin > 0); %% lower bound binding 0301 ku = find(lam.eqlin < 0); %% upper bound binding 0302 0303 mu_l(ieq(kl)) = lam.eqlin(kl); 0304 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0305 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0306 0307 mu_u(ieq(ku)) = -lam.eqlin(ku); 0308 mu_u(ilt) = -lam.ineqlin(1:nlt); 0309 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0310 0311 lambda = struct( ... 0312 'mu_l', mu_l, ... 0313 'mu_u', mu_u, ... 0314 'lower', lam.lower, ... 0315 'upper', lam.upper ... 0316 );