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, lambda] = 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 verbose 0034 % 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, lambda] = 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 2340 2014-06-27 19:15:10Z ray $ 0095 % by Ray Zimmerman, PSERC Cornell 0096 % Copyright (c) 2010-2013 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 0202 %% split up linear constraints 0203 ieq = find( abs(u-l) <= eps ); %% equality 0204 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0205 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0206 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0207 Ae = A(ieq, :); 0208 be = u(ieq); 0209 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0210 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0211 0212 %% grab some dimensions 0213 nlt = length(ilt); %% number of upper bounded linear inequalities 0214 ngt = length(igt); %% number of lower bounded linear inequalities 0215 nbx = length(ibx); %% number of doubly bounded linear inequalities 0216 0217 %% set up options struct for CPLEX 0218 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0219 cplex_opt = cplex_options(opt.cplex_opt); 0220 else 0221 cplex_opt = cplex_options; 0222 end 0223 0224 cplex = Cplex('null'); 0225 vstr = cplex.getVersion; 0226 [s,e,tE,m,t] = regexp(vstr, '(\d+\.\d+)\.'); 0227 vnum = str2num(t{1}{1}); 0228 vrb = max([0 verbose-1]); 0229 cplex_opt.barrier.display = vrb; 0230 cplex_opt.conflict.display = vrb; 0231 cplex_opt.mip.display = vrb; 0232 cplex_opt.sifting.display = vrb; 0233 cplex_opt.simplex.display = vrb; 0234 cplex_opt.tune.display = vrb; 0235 if vrb && vnum > 12.2 0236 cplex_opt.diagnostics = 'on'; 0237 end 0238 if verbose > 2 0239 cplex_opt.Display = 'iter'; 0240 elseif verbose > 1 0241 cplex_opt.Display = 'on'; 0242 elseif verbose > 0 0243 cplex_opt.Display = 'off'; 0244 end 0245 0246 if isempty(Ai) && isempty(Ae) 0247 unconstrained = 1; 0248 Ae = sparse(1, nx); 0249 be = 0; 0250 else 0251 unconstrained = 0; 0252 end 0253 0254 %% call the solver 0255 if verbose 0256 methods = { 0257 'default', 0258 'primal simplex', 0259 'dual simplex', 0260 'network simplex', 0261 'barrier', 0262 'sifting', 0263 'concurrent' 0264 }; 0265 end 0266 if isempty(H) || ~any(any(H)) 0267 if verbose 0268 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0269 vstr, methods{cplex_opt.lpmethod+1}); 0270 end 0271 [x, f, eflag, output, lam] = ... 0272 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0273 else 0274 if verbose 0275 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0276 vstr, methods{cplex_opt.qpmethod+1}); 0277 end 0278 %% ensure H is numerically symmetric 0279 if ~isequal(H, H') 0280 H = (H + H')/2; 0281 end 0282 [x, f, eflag, output, lam] = ... 0283 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0284 end 0285 0286 %% workaround for eflag == 5 (which we have seen return infeasible results) 0287 %% cplexstatus: 6 0288 %% cplexstatusstring: 'non-optimal' 0289 %% message: 'Solution with numerical issues' 0290 if eflag > 1 0291 warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n cplexstatus: %d\n cplexstatusstring: ''%s''\n message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message); 0292 eflag = -100 - eflag; 0293 end 0294 0295 %% check for empty results (in case optimization failed) 0296 if isempty(x) 0297 x = NaN(nx, 1); 0298 end 0299 if isempty(f) 0300 f = NaN; 0301 end 0302 if isempty(lam) 0303 lam.ineqlin = NaN(length(bi), 1); 0304 lam.eqlin = NaN(length(be), 1); 0305 lam.lower = NaN(nx, 1); 0306 lam.upper = NaN(nx, 1); 0307 mu_l = NaN(nA, 1); 0308 mu_u = NaN(nA, 1); 0309 else 0310 mu_l = zeros(nA, 1); 0311 mu_u = zeros(nA, 1); 0312 end 0313 if unconstrained 0314 lam.eqlin = []; 0315 end 0316 0317 %% negate prices depending on version 0318 if vnum < 12.3 0319 lam.eqlin = -lam.eqlin; 0320 lam.ineqlin = -lam.ineqlin; 0321 end 0322 0323 %% repackage lambdas 0324 kl = find(lam.eqlin < 0); %% lower bound binding 0325 ku = find(lam.eqlin > 0); %% upper bound binding 0326 0327 mu_l(ieq(kl)) = -lam.eqlin(kl); 0328 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0329 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0330 0331 mu_u(ieq(ku)) = lam.eqlin(ku); 0332 mu_u(ilt) = lam.ineqlin(1:nlt); 0333 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0334 0335 lambda = struct( ... 0336 'mu_l', mu_l, ... 0337 'mu_u', mu_u, ... 0338 'lower', lam.lower, ... 0339 'upper', lam.upper ... 0340 );