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 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0095 % by Ray Zimmerman, PSERC Cornell 0096 % 0097 % $Id: qps_cplex.m 2661 2015-03-20 17:02:46Z ray $ 0098 % 0099 % This file is part of MATPOWER. 0100 % Covered by the 3-clause BSD License (see LICENSE file for details). 0101 % See http://www.pserc.cornell.edu/matpower/ for more info. 0102 0103 %% check for CPLEX 0104 % if ~have_fcn('cplexqp') 0105 % error('qps_cplex: requires the Matlab interface for CPLEX'); 0106 % end 0107 0108 %%----- input argument handling ----- 0109 %% gather inputs 0110 if nargin == 1 && isstruct(H) %% problem struct 0111 p = H; 0112 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0113 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0114 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0115 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0116 if isfield(p, 'u'), u = p.u; else, u = []; end 0117 if isfield(p, 'l'), l = p.l; else, l = []; end 0118 if isfield(p, 'A'), A = p.A; else, A = []; end 0119 if isfield(p, 'c'), c = p.c; else, c = []; end 0120 if isfield(p, 'H'), H = p.H; else, H = []; end 0121 else %% individual args 0122 if nargin < 9 0123 opt = []; 0124 if nargin < 8 0125 x0 = []; 0126 if nargin < 7 0127 xmax = []; 0128 if nargin < 6 0129 xmin = []; 0130 end 0131 end 0132 end 0133 end 0134 end 0135 0136 %% define nx, set default values for missing optional inputs 0137 if isempty(H) || ~any(any(H)) 0138 if isempty(A) && isempty(xmin) && isempty(xmax) 0139 error('qps_cplex: LP problem must include constraints or variable bounds'); 0140 else 0141 if ~isempty(A) 0142 nx = size(A, 2); 0143 elseif ~isempty(xmin) 0144 nx = length(xmin); 0145 else % if ~isempty(xmax) 0146 nx = length(xmax); 0147 end 0148 end 0149 else 0150 nx = size(H, 1); 0151 end 0152 if isempty(c) 0153 c = zeros(nx, 1); 0154 end 0155 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0156 (isempty(u) || all(u == Inf))) 0157 A = sparse(0,nx); %% no limits => no linear constraints 0158 end 0159 nA = size(A, 1); %% number of original linear constraints 0160 if isempty(u) %% By default, linear inequalities are ... 0161 u = Inf(nA, 1); %% ... unbounded above and ... 0162 end 0163 if isempty(l) 0164 l = -Inf(nA, 1); %% ... unbounded below. 0165 end 0166 if isempty(xmin) %% By default, optimization variables are ... 0167 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0168 end 0169 if isempty(xmax) 0170 xmax = Inf(nx, 1); %% ... unbounded above. 0171 end 0172 if isempty(x0) 0173 x0 = zeros(nx, 1); 0174 end 0175 0176 %% default options 0177 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0178 verbose = opt.verbose; 0179 else 0180 verbose = 0; 0181 end 0182 0183 %% split up linear constraints 0184 ieq = find( abs(u-l) <= eps ); %% equality 0185 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0186 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0187 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0188 Ae = A(ieq, :); 0189 be = u(ieq); 0190 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0191 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0192 0193 %% grab some dimensions 0194 nlt = length(ilt); %% number of upper bounded linear inequalities 0195 ngt = length(igt); %% number of lower bounded linear inequalities 0196 nbx = length(ibx); %% number of doubly bounded linear inequalities 0197 0198 %% set up options struct for CPLEX 0199 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0200 cplex_opt = cplex_options(opt.cplex_opt); 0201 else 0202 cplex_opt = cplex_options; 0203 end 0204 0205 vstr = have_fcn('cplex', 'vstr'); 0206 vnum = have_fcn('cplex', 'vnum'); 0207 vrb = max([0 verbose-1]); 0208 cplex_opt.barrier.display = vrb; 0209 cplex_opt.conflict.display = vrb; 0210 cplex_opt.mip.display = vrb; 0211 cplex_opt.sifting.display = vrb; 0212 cplex_opt.simplex.display = vrb; 0213 cplex_opt.tune.display = vrb; 0214 if vrb && vnum > 12.002 0215 cplex_opt.diagnostics = 'on'; 0216 end 0217 if verbose > 2 0218 cplex_opt.Display = 'iter'; 0219 elseif verbose > 1 0220 cplex_opt.Display = 'on'; 0221 elseif verbose > 0 0222 cplex_opt.Display = 'off'; 0223 end 0224 0225 if isempty(Ai) && isempty(Ae) 0226 unconstrained = 1; 0227 Ae = sparse(1, nx); 0228 be = 0; 0229 else 0230 unconstrained = 0; 0231 end 0232 0233 %% call the solver 0234 if verbose 0235 methods = { 0236 'default', 0237 'primal simplex', 0238 'dual simplex', 0239 'network simplex', 0240 'barrier', 0241 'sifting', 0242 'concurrent' 0243 }; 0244 end 0245 if isempty(H) || ~any(any(H)) 0246 if verbose 0247 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0248 vstr, methods{cplex_opt.lpmethod+1}); 0249 end 0250 [x, f, eflag, output, lam] = ... 0251 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0252 else 0253 if verbose 0254 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0255 vstr, methods{cplex_opt.qpmethod+1}); 0256 end 0257 %% ensure H is numerically symmetric 0258 if ~isequal(H, H') 0259 H = (H + H')/2; 0260 end 0261 [x, f, eflag, output, lam] = ... 0262 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0263 end 0264 0265 %% workaround for eflag == 5 (which we have seen return infeasible results) 0266 %% cplexstatus: 6 0267 %% cplexstatusstring: 'non-optimal' 0268 %% message: 'Solution with numerical issues' 0269 if eflag > 1 0270 warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n cplexstatus: %d\n cplexstatusstring: ''%s''\n message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message); 0271 eflag = -100 - eflag; 0272 end 0273 0274 %% check for empty results (in case optimization failed) 0275 if isempty(x) 0276 x = NaN(nx, 1); 0277 end 0278 if isempty(f) 0279 f = NaN; 0280 end 0281 if isempty(lam) 0282 lam.ineqlin = NaN(length(bi), 1); 0283 lam.eqlin = NaN(length(be), 1); 0284 lam.lower = NaN(nx, 1); 0285 lam.upper = NaN(nx, 1); 0286 mu_l = NaN(nA, 1); 0287 mu_u = NaN(nA, 1); 0288 else 0289 mu_l = zeros(nA, 1); 0290 mu_u = zeros(nA, 1); 0291 end 0292 if unconstrained 0293 lam.eqlin = []; 0294 end 0295 0296 %% negate prices depending on version 0297 if vnum < 12.003 0298 lam.eqlin = -lam.eqlin; 0299 lam.ineqlin = -lam.ineqlin; 0300 end 0301 0302 %% repackage lambdas 0303 kl = find(lam.eqlin < 0); %% lower bound binding 0304 ku = find(lam.eqlin > 0); %% upper bound binding 0305 0306 mu_l(ieq(kl)) = -lam.eqlin(kl); 0307 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0308 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0309 0310 mu_u(ieq(ku)) = lam.eqlin(ku); 0311 mu_u(ilt) = lam.ineqlin(1:nlt); 0312 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0313 0314 lambda = struct( ... 0315 'mu_l', mu_l, ... 0316 'mu_u', mu_u, ... 0317 'lower', lam.lower, ... 0318 'upper', lam.upper ... 0319 );