MIQPS_GUROBI Mixed Integer Quadratic Program Solver based on GUROBI.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
MIQPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A wrapper function providing a MATPOWER standardized interface for using
GUROBI 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
3 = even more 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
grb_opt - options struct for GUROBI, 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 : GUROBI exit flag
1 = converged
0 or negative values = negative of GUROBI exit flag
(see GUROBI documentation for details)
OUTPUT : GUROBI output struct
(see GUROBI 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_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_gurobi(H, c, A, l, u)
x = miqps_gurobi(H, c, A, l, u, xmin, xmax)
x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0)
x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype)
x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
x = miqps_gurobi(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_gurobi(...)
[x, f] = miqps_gurobi(...)
[x, f, exitflag] = miqps_gurobi(...)
[x, f, exitflag, output] = miqps_gurobi(...)
[x, f, exitflag, output, lambda] = miqps_gurobi(...)
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_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt);
See also GUROBI.0001 function [x, f, eflag, output, lambda] = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0002 %MIQPS_GUROBI Mixed Integer Quadratic Program Solver based on GUROBI. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % MIQPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A wrapper function providing a MATPOWER standardized interface for using 0006 % GUROBI 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 % 3 = even more 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 % grb_opt - options struct for GUROBI, value in 0046 % verbose 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 : GUROBI exit flag 0055 % 1 = converged 0056 % 0 or negative values = negative of GUROBI exit flag 0057 % (see GUROBI documentation for details) 0058 % OUTPUT : GUROBI output struct 0059 % (see GUROBI documentation for details) 0060 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0061 % multipliers on the constraints, with fields: 0062 % mu_l - lower (left-hand) limit on linear constraints 0063 % mu_u - upper (right-hand) limit on linear constraints 0064 % lower - lower bound on optimization variables 0065 % upper - upper bound on optimization variables 0066 % 0067 % Note the calling syntax is almost identical to that of QUADPROG 0068 % from MathWorks' Optimization Toolbox. The main difference is that 0069 % the linear constraints are specified with A, L, U instead of 0070 % A, B, Aeq, Beq. 0071 % 0072 % Calling syntax options: 0073 % [x, f, exitflag, output, lambda] = ... 0074 % miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0075 % 0076 % x = miqps_gurobi(H, c, A, l, u) 0077 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax) 0078 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0) 0079 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype) 0080 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0081 % x = miqps_gurobi(problem), where problem is a struct with fields: 0082 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0083 % all fields except 'c', 'A' and 'l' or 'u' are optional 0084 % x = miqps_gurobi(...) 0085 % [x, f] = miqps_gurobi(...) 0086 % [x, f, exitflag] = miqps_gurobi(...) 0087 % [x, f, exitflag, output] = miqps_gurobi(...) 0088 % [x, f, exitflag, output, lambda] = miqps_gurobi(...) 0089 % 0090 % 0091 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0092 % H = [ 1003.1 4.3 6.3 5.9; 0093 % 4.3 2.2 2.1 3.9; 0094 % 6.3 2.1 3.5 4.8; 0095 % 5.9 3.9 4.8 10 ]; 0096 % c = zeros(4,1); 0097 % A = [ 1 1 1 1; 0098 % 0.17 0.11 0.10 0.18 ]; 0099 % l = [1; 0.10]; 0100 % u = [1; Inf]; 0101 % xmin = zeros(4,1); 0102 % x0 = [1; 0; 0; 1]; 0103 % opt = struct('verbose', 2); 0104 % [x, f, s, out, lambda] = miqps_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt); 0105 % 0106 % See also GUROBI. 0107 0108 % MATPOWER 0109 % Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC) 0110 % by Ray Zimmerman, PSERC Cornell 0111 % 0112 % This file is part of MATPOWER. 0113 % Covered by the 3-clause BSD License (see LICENSE file for details). 0114 % See https://matpower.org for more info. 0115 0116 %%----- input argument handling ----- 0117 %% gather inputs 0118 if nargin == 1 && isstruct(H) %% problem struct 0119 p = H; 0120 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0121 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0122 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0123 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0124 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0125 if isfield(p, 'u'), u = p.u; else, u = []; end 0126 if isfield(p, 'l'), l = p.l; else, l = []; end 0127 if isfield(p, 'A'), A = p.A; else, A = []; end 0128 if isfield(p, 'c'), c = p.c; else, c = []; end 0129 if isfield(p, 'H'), H = p.H; else, H = []; end 0130 else %% individual args 0131 if nargin < 10 0132 opt = []; 0133 if nargin < 9 0134 vtype = []; 0135 if nargin < 8 0136 x0 = []; 0137 if nargin < 7 0138 xmax = []; 0139 if nargin < 6 0140 xmin = []; 0141 end 0142 end 0143 end 0144 end 0145 end 0146 end 0147 0148 %% define nx, set default values for missing optional inputs 0149 if isempty(H) || ~any(any(H)) 0150 if isempty(A) && isempty(xmin) && isempty(xmax) 0151 error('miqps_gurobi: LP problem must include constraints or variable bounds'); 0152 else 0153 if ~isempty(A) 0154 nx = size(A, 2); 0155 elseif ~isempty(xmin) 0156 nx = length(xmin); 0157 else % if ~isempty(xmax) 0158 nx = length(xmax); 0159 end 0160 end 0161 else 0162 nx = size(H, 1); 0163 end 0164 if isempty(c) 0165 c = zeros(nx, 1); 0166 end 0167 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0168 (isempty(u) || all(u == Inf))) 0169 A = sparse(0,nx); %% no limits => no linear constraints 0170 end 0171 nA = size(A, 1); %% number of original linear constraints 0172 if isempty(u) %% By default, linear inequalities are ... 0173 u = Inf(nA, 1); %% ... unbounded above and ... 0174 end 0175 if isempty(l) 0176 l = -Inf(nA, 1); %% ... unbounded below. 0177 end 0178 if isempty(xmin) %% By default, optimization variables are ... 0179 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0180 end 0181 if isempty(xmax) 0182 xmax = Inf(nx, 1); %% ... unbounded above. 0183 end 0184 if isempty(x0) 0185 x0 = zeros(nx, 1); 0186 end 0187 0188 %% default options 0189 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0190 verbose = opt.verbose; 0191 else 0192 verbose = 0; 0193 end 0194 0195 %% set up options struct for Gurobi 0196 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0197 g_opt = gurobi_options(opt.grb_opt); 0198 else 0199 g_opt = gurobi_options; 0200 end 0201 if verbose > 1 0202 g_opt.LogToConsole = 1; 0203 g_opt.OutputFlag = 1; 0204 if verbose > 2 0205 g_opt.DisplayInterval = 1; 0206 else 0207 g_opt.DisplayInterval = 100; 0208 end 0209 else 0210 g_opt.LogToConsole = 0; 0211 g_opt.OutputFlag = 0; 0212 end 0213 0214 if ~issparse(A) 0215 A = sparse(A); 0216 end 0217 if issparse(c); 0218 c = full(c); 0219 end 0220 0221 %% split up linear constraints 0222 ieq = find( abs(u-l) <= eps ); %% equality 0223 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0224 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0225 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0226 0227 %% grab some dimensions 0228 nlt = length(ilt); %% number of upper bounded linear inequalities 0229 ngt = length(igt); %% number of lower bounded linear inequalities 0230 nbx = length(ibx); %% number of doubly bounded linear inequalities 0231 neq = length(ieq); %% number of equalities 0232 niq = nlt+ngt+2*nbx; %% number of inequalities 0233 0234 %% set up model 0235 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0236 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0237 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0238 m.lb = xmin; 0239 m.ub = xmax; 0240 m.obj = c'; 0241 if ~isempty(vtype) 0242 m.vtype = vtype; 0243 end 0244 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ... 0245 vtype == 'S' | vtype == 'N')) 0246 mi = 0; 0247 else 0248 mi = 1; 0249 end 0250 0251 %% call the solver 0252 if isempty(H) || ~any(any(H)) 0253 lpqp = 'LP'; 0254 else 0255 lpqp = 'QP'; 0256 if ~issparse(H) 0257 H = sparse(H); 0258 end 0259 m.Q = 0.5 * H; 0260 end 0261 if mi 0262 lpqp = ['MI' lpqp]; 0263 end 0264 if verbose 0265 alg_names = { 0266 'automatic', 0267 'primal simplex', 0268 'dual simplex', 0269 'interior point', 0270 'concurrent', 0271 'deterministic concurrent' 0272 }; 0273 vn = gurobiver; 0274 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0275 vn, alg_names{g_opt.Method+2}, lpqp); 0276 end 0277 results = gurobi(m, g_opt); 0278 switch results.status 0279 case 'LOADED', %% 1 0280 eflag = -1; 0281 case 'OPTIMAL', %% 2, optimal solution found 0282 eflag = 1; 0283 case 'INFEASIBLE', %% 3 0284 eflag = -3; 0285 case 'INF_OR_UNBD', %% 4 0286 eflag = -4; 0287 case 'UNBOUNDED', %% 5 0288 eflag = -5; 0289 case 'CUTOFF', %% 6 0290 eflag = -6; 0291 case 'ITERATION_LIMIT', %% 7 0292 eflag = -7; 0293 case 'NODE_LIMIT', %% 8 0294 eflag = -8; 0295 case 'TIME_LIMIT', %% 9 0296 eflag = -9; 0297 case 'SOLUTION_LIMIT', %% 10 0298 eflag = -10; 0299 case 'INTERRUPTED', %% 11 0300 eflag = -11; 0301 case 'NUMERIC', %% 12 0302 eflag = -12; 0303 case 'SUBOPTIMAL', %% 13 0304 eflag = -13; 0305 case 'INPROGRESS', %% 14 0306 eflag = -14; 0307 otherwise, 0308 eflag = 0; 0309 end 0310 output = results; 0311 0312 %% check for empty results (in case optimization failed) 0313 if ~isfield(results, 'x') || isempty(results.x) 0314 x = NaN(nx, 1); 0315 lam.lower = NaN(nx, 1); 0316 lam.upper = NaN(nx, 1); 0317 else 0318 x = results.x; 0319 lam.lower = zeros(nx, 1); 0320 lam.upper = zeros(nx, 1); 0321 end 0322 if ~isfield(results, 'objval') || isempty(results.objval) 0323 f = NaN; 0324 else 0325 f = results.objval; 0326 end 0327 if ~isfield(results, 'pi') || isempty(results.pi) 0328 pi = NaN(length(m.rhs), 1); 0329 else 0330 pi = results.pi; 0331 end 0332 if ~isfield(results, 'rc') || isempty(results.rc) 0333 rc = NaN(nx, 1); 0334 else 0335 rc = results.rc; 0336 end 0337 0338 kl = find(rc > 0); %% lower bound binding 0339 ku = find(rc < 0); %% upper bound binding 0340 lam.lower(kl) = rc(kl); 0341 lam.upper(ku) = -rc(ku); 0342 lam.eqlin = pi(1:neq); 0343 lam.ineqlin = pi(neq+(1:niq)); 0344 mu_l = zeros(nA, 1); 0345 mu_u = zeros(nA, 1); 0346 0347 %% repackage lambdas 0348 kl = find(lam.eqlin > 0); %% lower bound binding 0349 ku = find(lam.eqlin < 0); %% upper bound binding 0350 0351 mu_l(ieq(kl)) = lam.eqlin(kl); 0352 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0353 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0354 0355 mu_u(ieq(ku)) = -lam.eqlin(ku); 0356 mu_u(ilt) = -lam.ineqlin(1:nlt); 0357 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0358 0359 lambda = struct( ... 0360 'mu_l', mu_l, ... 0361 'mu_u', mu_u, ... 0362 'lower', lam.lower, ... 0363 'upper', lam.upper ... 0364 ); 0365 0366 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0367 if verbose 0368 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0369 end 0370 if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol) 0371 tol = opt.price_stage_warn_tol; 0372 else 0373 tol = 1e-7; 0374 end 0375 if length(vtype) == 1 0376 if vtype == 'I' || vtype == 'B' || vtype == 'N' 0377 k = (1:nx); 0378 elseif vtype == 'S' 0379 k = find(x == 0); 0380 end 0381 else %% length(vtype) == nx 0382 k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ... 0383 (vtype == 'S' & x' == 0)); 0384 end 0385 0386 x(k) = round(x(k)); 0387 xmin(k) = x(k); 0388 xmax(k) = x(k); 0389 x0 = x; 0390 % opt.grb_opt.Method = 0; %% primal simplex 0391 0392 [x_, f_, eflag_, output_, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt); 0393 if eflag ~= eflag_ 0394 error('miqps_gurobi: EXITFLAG from price computation stage = %d', eflag_); 0395 end 0396 if abs(f - f_)/max(abs(f), 1) > tol 0397 warning('miqps_gurobi: 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_gurobi: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0404 end 0405 output.price_stage = output_; 0406 end