MIQPS_CPLEX Mixed Integer Quadratic Program Solver based on CPLEX.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_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
VTYPE : character string of length NX (number of elements in X),
or 1 (value applies to all variables in x),
allowed values are 'C' (continuous), 'B' (binary),
'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
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
skip_prices (0) - flag that specifies whether or not to
skip the price computation stage, in which the problem
is re-solved for only the continuous variables, with all
others being constrained to their solved values
price_stage_warn_tol (1e-7) - tolerance on the objective fcn
value and primal variable relative match required to avoid
mis-match warning message
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, vtype, 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] = ...
miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_cplex(H, c, A, l, u)
x = miqps_cplex(H, c, A, l, u, xmin, xmax)
x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0)
x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_cplex(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, vtype, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = miqps_cplex(...)
[x, f] = miqps_cplex(...)
[x, f, exitflag] = miqps_cplex(...)
[x, f, exitflag, output] = miqps_cplex(...)
[x, f, exitflag, output, lambda] = miqps_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] = miqps_cplex(H, c, A, l, u, xmin, [], x0, vtype, opt);
See also MIQPS_MASTER, CPLEXMIQP, CPLEXMILP, CPLEXQP, CPLEXLP,
CPLEX_OPTIONS.0001 function [x, f, eflag, output, lambda] = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_CPLEX Mixed Integer Quadratic Program Solver based on CPLEX. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_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 % VTYPE : character string of length NX (number of elements in X), 0028 % or 1 (value applies to all variables in x), 0029 % allowed values are 'C' (continuous), 'B' (binary), 0030 % 'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer). 0031 % OPT : optional options structure with the following fields, 0032 % all of which are also optional (default values shown in 0033 % parentheses) 0034 % verbose (0) - controls level of progress output displayed 0035 % 0 = no progress output 0036 % 1 = some progress output 0037 % 2 = verbose progress output 0038 % skip_prices (0) - flag that specifies whether or not to 0039 % skip the price computation stage, in which the problem 0040 % is re-solved for only the continuous variables, with all 0041 % others being constrained to their solved values 0042 % price_stage_warn_tol (1e-7) - tolerance on the objective fcn 0043 % value and primal variable relative match required to avoid 0044 % mis-match warning message 0045 % cplex_opt - options struct for CPLEX, value in verbose 0046 % overrides these options 0047 % PROBLEM : The inputs can alternatively be supplied in a single 0048 % PROBLEM struct with fields corresponding to the input arguments 0049 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0050 % 0051 % Outputs: 0052 % X : solution vector 0053 % F : final objective function value 0054 % EXITFLAG : CPLEXQP/CPLEXLP exit flag 0055 % (see CPLEXQP and CPLEXLP documentation for details) 0056 % OUTPUT : CPLEXQP/CPLEXLP output struct 0057 % (see CPLEXQP and CPLEXLP documentation for details) 0058 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0059 % multipliers on the constraints, with fields: 0060 % mu_l - lower (left-hand) limit on linear constraints 0061 % mu_u - upper (right-hand) limit on linear constraints 0062 % lower - lower bound on optimization variables 0063 % upper - upper bound on optimization variables 0064 % 0065 % Note the calling syntax is almost identical to that of QUADPROG 0066 % from MathWorks' Optimization Toolbox. The main difference is that 0067 % the linear constraints are specified with A, L, U instead of 0068 % A, B, Aeq, Beq. 0069 % 0070 % Calling syntax options: 0071 % [x, f, exitflag, output, lambda] = ... 0072 % miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0073 % 0074 % x = miqps_cplex(H, c, A, l, u) 0075 % x = miqps_cplex(H, c, A, l, u, xmin, xmax) 0076 % x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0) 0077 % x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype) 0078 % x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0079 % x = miqps_cplex(problem), where problem is a struct with fields: 0080 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0081 % all fields except 'c', 'A' and 'l' or 'u' are optional 0082 % x = miqps_cplex(...) 0083 % [x, f] = miqps_cplex(...) 0084 % [x, f, exitflag] = miqps_cplex(...) 0085 % [x, f, exitflag, output] = miqps_cplex(...) 0086 % [x, f, exitflag, output, lambda] = miqps_cplex(...) 0087 % 0088 % 0089 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0090 % H = [ 1003.1 4.3 6.3 5.9; 0091 % 4.3 2.2 2.1 3.9; 0092 % 6.3 2.1 3.5 4.8; 0093 % 5.9 3.9 4.8 10 ]; 0094 % c = zeros(4,1); 0095 % A = [ 1 1 1 1; 0096 % 0.17 0.11 0.10 0.18 ]; 0097 % l = [1; 0.10]; 0098 % u = [1; Inf]; 0099 % xmin = zeros(4,1); 0100 % x0 = [1; 0; 0; 1]; 0101 % opt = struct('verbose', 2); 0102 % [x, f, s, out, lambda] = miqps_cplex(H, c, A, l, u, xmin, [], x0, vtype, opt); 0103 % 0104 % See also MIQPS_MASTER, CPLEXMIQP, CPLEXMILP, CPLEXQP, CPLEXLP, 0105 % CPLEX_OPTIONS. 0106 0107 % MP-Opt-Model 0108 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0109 % by Ray Zimmerman, PSERC Cornell 0110 % 0111 % This file is part of MP-Opt-Model. 0112 % Covered by the 3-clause BSD License (see LICENSE file for details). 0113 % See https://github.com/MATPOWER/mp-opt-model for more info. 0114 0115 %% check for CPLEX 0116 % if ~have_feature('cplexqp') 0117 % error('miqps_cplex: requires the MATLAB interface for CPLEX'); 0118 % end 0119 0120 %%----- input argument handling ----- 0121 %% gather inputs 0122 if nargin == 1 && isstruct(H) %% problem struct 0123 p = H; 0124 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0125 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0126 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0127 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0128 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0129 if isfield(p, 'u'), u = p.u; else, u = []; end 0130 if isfield(p, 'l'), l = p.l; else, l = []; end 0131 if isfield(p, 'A'), A = p.A; else, A = []; end 0132 if isfield(p, 'c'), c = p.c; else, c = []; end 0133 if isfield(p, 'H'), H = p.H; else, H = []; end 0134 else %% individual args 0135 if nargin < 10 0136 opt = []; 0137 if nargin < 9 0138 vtype = []; 0139 if nargin < 8 0140 x0 = []; 0141 if nargin < 7 0142 xmax = []; 0143 if nargin < 6 0144 xmin = []; 0145 end 0146 end 0147 end 0148 end 0149 end 0150 end 0151 0152 %% define nx, set default values for missing optional inputs 0153 if isempty(H) || ~any(any(H)) 0154 if isempty(A) && isempty(xmin) && isempty(xmax) 0155 error('miqps_cplex: LP problem must include constraints or variable bounds'); 0156 else 0157 if ~isempty(A) 0158 nx = size(A, 2); 0159 elseif ~isempty(xmin) 0160 nx = length(xmin); 0161 else % if ~isempty(xmax) 0162 nx = length(xmax); 0163 end 0164 end 0165 else 0166 nx = size(H, 1); 0167 end 0168 if isempty(c) 0169 c = zeros(nx, 1); 0170 end 0171 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0172 (isempty(u) || all(u == Inf))) 0173 A = sparse(0,nx); %% no limits => no linear constraints 0174 end 0175 nA = size(A, 1); %% number of original linear constraints 0176 if isempty(u) %% By default, linear inequalities are ... 0177 u = Inf(nA, 1); %% ... unbounded above and ... 0178 end 0179 if isempty(l) 0180 l = -Inf(nA, 1); %% ... unbounded below. 0181 end 0182 if isempty(xmin) %% By default, optimization variables are ... 0183 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0184 end 0185 if isempty(xmax) 0186 xmax = Inf(nx, 1); %% ... unbounded above. 0187 end 0188 if isempty(x0) 0189 x0 = zeros(nx, 1); 0190 end 0191 0192 %% default options 0193 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0194 verbose = opt.verbose; 0195 else 0196 verbose = 0; 0197 end 0198 0199 %% split up linear constraints 0200 ieq = find( abs(u-l) <= eps ); %% equality 0201 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0202 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0203 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0204 Ae = A(ieq, :); 0205 be = u(ieq); 0206 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0207 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0208 0209 %% grab some dimensions 0210 nlt = length(ilt); %% number of upper bounded linear inequalities 0211 ngt = length(igt); %% number of lower bounded linear inequalities 0212 nbx = length(ibx); %% number of doubly bounded linear inequalities 0213 0214 %% set up options struct for CPLEX 0215 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0216 cplex_opt = cplex_options(opt.cplex_opt); 0217 else 0218 cplex_opt = cplex_options; 0219 end 0220 0221 vstr = have_feature('cplex', 'vstr'); 0222 vnum = have_feature('cplex', 'vnum'); 0223 vrb = max([0 verbose-1]); 0224 if vrb && vnum > 12.002 && vnum < 12.007 0225 cplex_opt.diagnostics = 'on'; 0226 end 0227 if verbose > 2 0228 cplex_opt.display = 'iter'; 0229 elseif verbose > 1 0230 cplex_opt.display = 'on'; 0231 elseif verbose > 0 0232 cplex_opt.display = 'off'; 0233 end 0234 0235 if isempty(Ai) && isempty(Ae) 0236 unconstrained = 1; 0237 Ae = sparse(1, nx); 0238 be = 0; 0239 else 0240 unconstrained = 0; 0241 end 0242 0243 %% call the solver 0244 if verbose 0245 alg_names = { 0246 'default', 0247 'primal simplex', 0248 'dual simplex', 0249 'network simplex', 0250 'barrier', 0251 'sifting', 0252 'concurrent' 0253 }; 0254 end 0255 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ... 0256 vtype == 'S' | vtype == 'N')) 0257 mi = 0; 0258 else 0259 mi = 1; 0260 %% expand vtype to nx elements if necessary 0261 if length(vtype) == 1 && nx > 1 0262 vtype = char(vtype * ones(1,nx)); 0263 end 0264 end 0265 0266 if mi 0267 if isempty(H) || ~any(any(H)) 0268 if verbose 0269 fprintf('CPLEX Version %s -- %s MILP solver\n', ... 0270 vstr, alg_names{cplex_opt.lpmethod+1}); 0271 end 0272 [x, f, eflag, output] = ... 0273 cplexmilp(c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt); 0274 lam = []; 0275 else 0276 if verbose 0277 fprintf('CPLEX Version %s -- %s MIQP solver\n', ... 0278 vstr, alg_names{cplex_opt.qpmethod+1}); 0279 end 0280 %% ensure H is numerically symmetric 0281 if ~isequal(H, H') 0282 H = (H + H')/2; 0283 end 0284 [x, f, eflag, output] = ... 0285 cplexmiqp(H, c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt); 0286 end 0287 lam = []; 0288 else 0289 if isempty(H) || ~any(any(H)) 0290 if verbose 0291 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0292 vstr, alg_names{cplex_opt.lpmethod+1}); 0293 end 0294 [x, f, eflag, output, lam] = ... 0295 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0296 else 0297 if verbose 0298 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0299 vstr, alg_names{cplex_opt.qpmethod+1}); 0300 end 0301 %% ensure H is numerically symmetric 0302 if ~isequal(H, H') 0303 H = (H + H')/2; 0304 end 0305 [x, f, eflag, output, lam] = ... 0306 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0307 end 0308 end 0309 0310 %% workaround for eflag == 5 (which we have seen return infeasible results) 0311 %% cplexstatus: 6 0312 %% cplexstatusstring: 'non-optimal' 0313 %% message: 'Solution with numerical issues' 0314 if eflag > 1 0315 warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n cplexstatus: %d\n cplexstatusstring: ''%s''\n message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message); 0316 if eflag == 5 && mi 0317 eflag = 1; %% give it a try for the MI phase 0318 else 0319 eflag = -100 - eflag; 0320 end 0321 end 0322 0323 %% check for empty results (in case optimization failed) 0324 if isempty(x) 0325 x = NaN(nx, 1); 0326 end 0327 if isempty(f) 0328 f = NaN; 0329 end 0330 if isempty(lam) 0331 lam.ineqlin = NaN(length(bi), 1); 0332 lam.eqlin = NaN(length(be), 1); 0333 lam.lower = NaN(nx, 1); 0334 lam.upper = NaN(nx, 1); 0335 mu_l = NaN(nA, 1); 0336 mu_u = NaN(nA, 1); 0337 else 0338 mu_l = zeros(nA, 1); 0339 mu_u = zeros(nA, 1); 0340 end 0341 if unconstrained 0342 lam.eqlin = []; 0343 end 0344 0345 %% negate prices depending on version 0346 if vnum < 12.003 0347 lam.eqlin = -lam.eqlin; 0348 lam.ineqlin = -lam.ineqlin; 0349 end 0350 0351 %% repackage lambdas 0352 kl = find(lam.eqlin < 0); %% lower bound binding 0353 ku = find(lam.eqlin > 0); %% upper bound binding 0354 0355 mu_l(ieq(kl)) = -lam.eqlin(kl); 0356 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0357 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0358 0359 mu_u(ieq(ku)) = lam.eqlin(ku); 0360 mu_u(ilt) = lam.ineqlin(1:nlt); 0361 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0362 0363 lambda = struct( ... 0364 'mu_l', mu_l, ... 0365 'mu_u', mu_u, ... 0366 'lower', lam.lower, ... 0367 'upper', lam.upper ... 0368 ); 0369 0370 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0371 if verbose 0372 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0373 end 0374 if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol) 0375 tol = opt.price_stage_warn_tol; 0376 else 0377 tol = 1e-7; 0378 end 0379 k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ... 0380 (vtype == 'S' & x' == 0)); 0381 x(k) = round(x(k)); 0382 xmin(k) = x(k); 0383 xmax(k) = x(k); 0384 x0 = x; 0385 opt.cplex_opt.lpmethod = 1; %% primal simplex 0386 opt.cplex_opt.qpmethod = 1; %% primal simplex 0387 0388 [x_, f_, eflag_, output_, lambda] = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt); 0389 if eflag ~= eflag_ 0390 error('miqps_cplex: EXITFLAG from price computation stage = %d', eflag_); 0391 end 0392 if abs(f - f_)/max(abs(f), 1) > tol 0393 warning('miqps_cplex: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1)); 0394 end 0395 xn = x; 0396 xn(abs(xn)<1) = 1; 0397 [mx, k] = max(abs(x - x_) ./ xn); 0398 if mx > tol 0399 warning('miqps_cplex: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0400 end 0401 output.price_stage = output_; 0402 end