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