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