QPS_CPLEX Quadratic Program Solver based on CPLEX.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CPLEX(PROBLEM)
A wrapper function providing a 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 https://v8doc.sas.com/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 QPS_MASTER, 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 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_CPLEX(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % CPLEXQP or CPLEXLP to solve the following QP (quadratic programming) 0008 % problem: 0009 % 0010 % min 1/2 X'*H*X + C'*X 0011 % X 0012 % 0013 % subject to 0014 % 0015 % L <= A*X <= U (linear constraints) 0016 % XMIN <= X <= XMAX (variable bounds) 0017 % 0018 % Inputs (all optional except H, C, A and L): 0019 % H : matrix (possibly sparse) of quadratic cost coefficients 0020 % C : vector of linear cost coefficients 0021 % A, L, U : define the optional linear constraints. Default 0022 % values for the elements of L and U are -Inf and Inf, 0023 % respectively. 0024 % XMIN, XMAX : optional lower and upper bounds on the 0025 % X variables, defaults are -Inf and Inf, respectively. 0026 % X0 : optional starting value of optimization vector X 0027 % OPT : optional options structure with the following fields, 0028 % all of which are also optional (default values shown in 0029 % parentheses) 0030 % verbose (0) - controls level of progress output displayed 0031 % 0 = no progress output 0032 % 1 = some progress output 0033 % 2 = verbose progress output 0034 % cplex_opt - options struct for CPLEX, value in verbose 0035 % overrides these options 0036 % PROBLEM : The inputs can alternatively be supplied in a single 0037 % PROBLEM struct with fields corresponding to the input arguments 0038 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0039 % 0040 % Outputs: 0041 % X : solution vector 0042 % F : final objective function value 0043 % EXITFLAG : CPLEXQP/CPLEXLP exit flag 0044 % (see CPLEXQP and CPLEXLP documentation for details) 0045 % OUTPUT : CPLEXQP/CPLEXLP output struct 0046 % (see CPLEXQP and CPLEXLP documentation for details) 0047 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0048 % multipliers on the constraints, with fields: 0049 % mu_l - lower (left-hand) limit on linear constraints 0050 % mu_u - upper (right-hand) limit on linear constraints 0051 % lower - lower bound on optimization variables 0052 % upper - upper bound on optimization variables 0053 % 0054 % Note the calling syntax is almost identical to that of QUADPROG 0055 % from MathWorks' Optimization Toolbox. The main difference is that 0056 % the linear constraints are specified with A, L, U instead of 0057 % A, B, Aeq, Beq. 0058 % 0059 % Calling syntax options: 0060 % [x, f, exitflag, output, lambda] = ... 0061 % qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 0062 % 0063 % x = qps_cplex(H, c, A, l, u) 0064 % x = qps_cplex(H, c, A, l, u, xmin, xmax) 0065 % x = qps_cplex(H, c, A, l, u, xmin, xmax, x0) 0066 % x = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt) 0067 % x = qps_cplex(problem), where problem is a struct with fields: 0068 % H, c, A, l, u, xmin, xmax, x0, opt 0069 % all fields except 'c', 'A' and 'l' or 'u' are optional 0070 % x = qps_cplex(...) 0071 % [x, f] = qps_cplex(...) 0072 % [x, f, exitflag] = qps_cplex(...) 0073 % [x, f, exitflag, output] = qps_cplex(...) 0074 % [x, f, exitflag, output, lambda] = qps_cplex(...) 0075 % 0076 % 0077 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0078 % H = [ 1003.1 4.3 6.3 5.9; 0079 % 4.3 2.2 2.1 3.9; 0080 % 6.3 2.1 3.5 4.8; 0081 % 5.9 3.9 4.8 10 ]; 0082 % c = zeros(4,1); 0083 % A = [ 1 1 1 1; 0084 % 0.17 0.11 0.10 0.18 ]; 0085 % l = [1; 0.10]; 0086 % u = [1; Inf]; 0087 % xmin = zeros(4,1); 0088 % x0 = [1; 0; 0; 1]; 0089 % opt = struct('verbose', 2); 0090 % [x, f, s, out, lambda] = qps_cplex(H, c, A, l, u, xmin, [], x0, opt); 0091 % 0092 % See also QPS_MASTER, CPLEXQP, CPLEXLP, CPLEX_OPTIONS. 0093 0094 % MP-Opt-Model 0095 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0096 % by Ray Zimmerman, PSERC Cornell 0097 % 0098 % This file is part of MP-Opt-Model. 0099 % Covered by the 3-clause BSD License (see LICENSE file for details). 0100 % See https://github.com/MATPOWER/mp-opt-model for more info. 0101 0102 %% check for CPLEX 0103 % if ~have_feature('cplexqp') 0104 % error('qps_cplex: requires the MATLAB interface for CPLEX'); 0105 % end 0106 0107 %%----- input argument handling ----- 0108 %% gather inputs 0109 if nargin == 1 && isstruct(H) %% problem struct 0110 p = H; 0111 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0112 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0113 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0114 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0115 if isfield(p, 'u'), u = p.u; else, u = []; end 0116 if isfield(p, 'l'), l = p.l; else, l = []; end 0117 if isfield(p, 'A'), A = p.A; else, A = []; end 0118 if isfield(p, 'c'), c = p.c; else, c = []; end 0119 if isfield(p, 'H'), H = p.H; else, H = []; end 0120 else %% individual args 0121 if nargin < 9 0122 opt = []; 0123 if nargin < 8 0124 x0 = []; 0125 if nargin < 7 0126 xmax = []; 0127 if nargin < 6 0128 xmin = []; 0129 end 0130 end 0131 end 0132 end 0133 end 0134 0135 %% define nx, set default values for missing optional inputs 0136 if isempty(H) || ~any(any(H)) 0137 if isempty(A) && isempty(xmin) && isempty(xmax) 0138 error('qps_cplex: LP problem must include constraints or variable bounds'); 0139 else 0140 if ~isempty(A) 0141 nx = size(A, 2); 0142 elseif ~isempty(xmin) 0143 nx = length(xmin); 0144 else % if ~isempty(xmax) 0145 nx = length(xmax); 0146 end 0147 end 0148 else 0149 nx = size(H, 1); 0150 end 0151 if isempty(c) 0152 c = zeros(nx, 1); 0153 end 0154 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0155 (isempty(u) || all(u == Inf))) 0156 A = sparse(0,nx); %% no limits => no linear constraints 0157 end 0158 nA = size(A, 1); %% number of original linear constraints 0159 if isempty(u) %% By default, linear inequalities are ... 0160 u = Inf(nA, 1); %% ... unbounded above and ... 0161 end 0162 if isempty(l) 0163 l = -Inf(nA, 1); %% ... unbounded below. 0164 end 0165 if isempty(xmin) %% By default, optimization variables are ... 0166 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0167 end 0168 if isempty(xmax) 0169 xmax = Inf(nx, 1); %% ... unbounded above. 0170 end 0171 if isempty(x0) 0172 x0 = zeros(nx, 1); 0173 end 0174 0175 %% default options 0176 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0177 verbose = opt.verbose; 0178 else 0179 verbose = 0; 0180 end 0181 0182 %% split up linear constraints 0183 ieq = find( abs(u-l) <= eps ); %% equality 0184 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0185 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0186 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0187 Ae = A(ieq, :); 0188 be = u(ieq); 0189 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0190 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0191 0192 %% grab some dimensions 0193 nlt = length(ilt); %% number of upper bounded linear inequalities 0194 ngt = length(igt); %% number of lower bounded linear inequalities 0195 nbx = length(ibx); %% number of doubly bounded linear inequalities 0196 0197 %% set up options struct for CPLEX 0198 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0199 cplex_opt = cplex_options(opt.cplex_opt); 0200 else 0201 cplex_opt = cplex_options; 0202 end 0203 0204 vstr = have_feature('cplex', 'vstr'); 0205 vnum = have_feature('cplex', 'vnum'); 0206 vrb = max([0 verbose-1]); 0207 if vrb && vnum > 12.002 && vnum < 12.007 0208 cplex_opt.diagnostics = 'on'; 0209 end 0210 if verbose > 2 0211 cplex_opt.display = 'iter'; 0212 elseif verbose > 1 0213 cplex_opt.display = 'on'; 0214 elseif verbose > 0 0215 cplex_opt.display = 'off'; 0216 end 0217 0218 if isempty(Ai) && isempty(Ae) 0219 unconstrained = 1; 0220 Ae = sparse(1, nx); 0221 be = 0; 0222 else 0223 unconstrained = 0; 0224 end 0225 0226 %% call the solver 0227 if verbose 0228 alg_names = { 0229 'default', 0230 'primal simplex', 0231 'dual simplex', 0232 'network simplex', 0233 'barrier', 0234 'sifting', 0235 'concurrent' 0236 }; 0237 end 0238 if isempty(H) || ~any(any(H)) 0239 if verbose 0240 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0241 vstr, alg_names{cplex_opt.lpmethod+1}); 0242 end 0243 [x, f, eflag, output, lam] = ... 0244 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0245 else 0246 if verbose 0247 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0248 vstr, alg_names{cplex_opt.qpmethod+1}); 0249 end 0250 %% ensure H is numerically symmetric 0251 if ~isequal(H, H') 0252 H = (H + H')/2; 0253 end 0254 [x, f, eflag, output, lam] = ... 0255 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0256 end 0257 0258 %% workaround for eflag == 5 (which we have seen return infeasible results) 0259 %% cplexstatus: 6 0260 %% cplexstatusstring: 'non-optimal' 0261 %% message: 'Solution with numerical issues' 0262 if eflag > 1 0263 warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n cplexstatus: %d\n cplexstatusstring: ''%s''\n message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message); 0264 eflag = -100 - eflag; 0265 end 0266 0267 %% check for empty results (in case optimization failed) 0268 if isempty(x) 0269 x = NaN(nx, 1); 0270 end 0271 if isempty(f) 0272 f = NaN; 0273 end 0274 if isempty(lam) 0275 lam.ineqlin = NaN(length(bi), 1); 0276 lam.eqlin = NaN(length(be), 1); 0277 lam.lower = NaN(nx, 1); 0278 lam.upper = NaN(nx, 1); 0279 mu_l = NaN(nA, 1); 0280 mu_u = NaN(nA, 1); 0281 else 0282 mu_l = zeros(nA, 1); 0283 mu_u = zeros(nA, 1); 0284 end 0285 if unconstrained 0286 lam.eqlin = []; 0287 end 0288 0289 %% negate prices depending on version 0290 if vnum < 12.003 0291 lam.eqlin = -lam.eqlin; 0292 lam.ineqlin = -lam.ineqlin; 0293 end 0294 0295 %% repackage lambdas 0296 kl = find(lam.eqlin < 0); %% lower bound binding 0297 ku = find(lam.eqlin > 0); %% upper bound binding 0298 0299 mu_l(ieq(kl)) = -lam.eqlin(kl); 0300 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0301 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0302 0303 mu_u(ieq(ku)) = lam.eqlin(ku); 0304 mu_u(ilt) = lam.ineqlin(1:nlt); 0305 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0306 0307 lambda = struct( ... 0308 'mu_l', mu_l, ... 0309 'mu_u', mu_u, ... 0310 'lower', lam.lower, ... 0311 'upper', lam.upper ... 0312 );