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