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, VTYPE, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GUROBI(PROBLEM)
A wrapper function providing a 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 MIQPS_MASTER, 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, VTYPE, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GUROBI(PROBLEM) 0006 % A wrapper function providing a standardized interface for using 0007 % GUROBI 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 % 3 = even more verbose progress output 0039 % skip_prices (0) - flag that specifies whether or not to 0040 % skip the price computation stage, in which the problem 0041 % is re-solved for only the continuous variables, with all 0042 % others being constrained to their solved values 0043 % price_stage_warn_tol (1e-7) - tolerance on the objective fcn 0044 % value and primal variable relative match required to avoid 0045 % mis-match warning message 0046 % grb_opt - options struct for GUROBI, value in verbose 0047 % overrides these options 0048 % PROBLEM : The inputs can alternatively be supplied in a single 0049 % PROBLEM struct with fields corresponding to the input arguments 0050 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0051 % 0052 % Outputs: 0053 % X : solution vector 0054 % F : final objective function value 0055 % EXITFLAG : GUROBI exit flag 0056 % 1 = converged 0057 % 0 or negative values = negative of GUROBI exit flag 0058 % (see GUROBI documentation for details) 0059 % OUTPUT : GUROBI output struct 0060 % (see GUROBI documentation for details) 0061 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0062 % multipliers on the constraints, with fields: 0063 % mu_l - lower (left-hand) limit on linear constraints 0064 % mu_u - upper (right-hand) limit on linear constraints 0065 % lower - lower bound on optimization variables 0066 % upper - upper bound on optimization variables 0067 % 0068 % Note the calling syntax is almost identical to that of QUADPROG 0069 % from MathWorks' Optimization Toolbox. The main difference is that 0070 % the linear constraints are specified with A, L, U instead of 0071 % A, B, Aeq, Beq. 0072 % 0073 % Calling syntax options: 0074 % [x, f, exitflag, output, lambda] = ... 0075 % miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0076 % 0077 % x = miqps_gurobi(H, c, A, l, u) 0078 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax) 0079 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0) 0080 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype) 0081 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0082 % x = miqps_gurobi(problem), where problem is a struct with fields: 0083 % H, c, A, l, u, xmin, xmax, x0, vtype, opt 0084 % all fields except 'c', 'A' and 'l' or 'u' are optional 0085 % x = miqps_gurobi(...) 0086 % [x, f] = miqps_gurobi(...) 0087 % [x, f, exitflag] = miqps_gurobi(...) 0088 % [x, f, exitflag, output] = miqps_gurobi(...) 0089 % [x, f, exitflag, output, lambda] = miqps_gurobi(...) 0090 % 0091 % 0092 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0093 % H = [ 1003.1 4.3 6.3 5.9; 0094 % 4.3 2.2 2.1 3.9; 0095 % 6.3 2.1 3.5 4.8; 0096 % 5.9 3.9 4.8 10 ]; 0097 % c = zeros(4,1); 0098 % A = [ 1 1 1 1; 0099 % 0.17 0.11 0.10 0.18 ]; 0100 % l = [1; 0.10]; 0101 % u = [1; Inf]; 0102 % xmin = zeros(4,1); 0103 % x0 = [1; 0; 0; 1]; 0104 % opt = struct('verbose', 2); 0105 % [x, f, s, out, lambda] = miqps_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt); 0106 % 0107 % See also MIQPS_MASTER, GUROBI. 0108 0109 % MP-Opt-Model 0110 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0111 % by Ray Zimmerman, PSERC Cornell 0112 % 0113 % This file is part of MP-Opt-Model. 0114 % Covered by the 3-clause BSD License (see LICENSE file for details). 0115 % See https://github.com/MATPOWER/mp-opt-model for more info. 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_gurobi: 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 %% set up options struct for Gurobi 0197 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0198 g_opt = gurobi_options(opt.grb_opt); 0199 else 0200 g_opt = gurobi_options; 0201 end 0202 if verbose > 1 0203 g_opt.LogToConsole = 1; 0204 g_opt.OutputFlag = 1; 0205 if verbose > 2 0206 g_opt.DisplayInterval = 1; 0207 else 0208 g_opt.DisplayInterval = 100; 0209 end 0210 else 0211 g_opt.LogToConsole = 0; 0212 g_opt.OutputFlag = 0; 0213 end 0214 0215 if ~issparse(A) 0216 A = sparse(A); 0217 end 0218 if issparse(c); 0219 c = full(c); 0220 end 0221 0222 %% split up linear constraints 0223 ieq = find( abs(u-l) <= eps ); %% equality 0224 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0225 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0226 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0227 0228 %% grab some dimensions 0229 nlt = length(ilt); %% number of upper bounded linear inequalities 0230 ngt = length(igt); %% number of lower bounded linear inequalities 0231 nbx = length(ibx); %% number of doubly bounded linear inequalities 0232 neq = length(ieq); %% number of equalities 0233 niq = nlt+ngt+2*nbx; %% number of inequalities 0234 0235 %% set up model 0236 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0237 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0238 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0239 m.lb = xmin; 0240 m.ub = xmax; 0241 m.obj = c'; 0242 if ~isempty(vtype) 0243 m.vtype = vtype; 0244 end 0245 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ... 0246 vtype == 'S' | vtype == 'N')) 0247 mi = 0; 0248 else 0249 mi = 1; 0250 end 0251 0252 %% call the solver 0253 if isempty(H) || ~any(any(H)) 0254 lpqp = 'LP'; 0255 else 0256 lpqp = 'QP'; 0257 if ~issparse(H) 0258 H = sparse(H); 0259 end 0260 m.Q = 0.5 * H; 0261 end 0262 if mi 0263 lpqp = ['MI' lpqp]; 0264 end 0265 if verbose 0266 alg_names = { 0267 'automatic', 0268 'primal simplex', 0269 'dual simplex', 0270 'interior point', 0271 'concurrent', 0272 'deterministic concurrent' 0273 }; 0274 vn = gurobiver; 0275 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0276 vn, alg_names{g_opt.Method+2}, lpqp); 0277 end 0278 results = gurobi(m, g_opt); 0279 switch results.status 0280 case 'LOADED', %% 1 0281 eflag = -1; 0282 case 'OPTIMAL', %% 2, optimal solution found 0283 eflag = 1; 0284 case 'INFEASIBLE', %% 3 0285 eflag = -3; 0286 case 'INF_OR_UNBD', %% 4 0287 eflag = -4; 0288 case 'UNBOUNDED', %% 5 0289 eflag = -5; 0290 case 'CUTOFF', %% 6 0291 eflag = -6; 0292 case 'ITERATION_LIMIT', %% 7 0293 eflag = -7; 0294 case 'NODE_LIMIT', %% 8 0295 eflag = -8; 0296 case 'TIME_LIMIT', %% 9 0297 eflag = -9; 0298 case 'SOLUTION_LIMIT', %% 10 0299 eflag = -10; 0300 case 'INTERRUPTED', %% 11 0301 eflag = -11; 0302 case 'NUMERIC', %% 12 0303 eflag = -12; 0304 case 'SUBOPTIMAL', %% 13 0305 eflag = -13; 0306 case 'INPROGRESS', %% 14 0307 eflag = -14; 0308 otherwise, 0309 eflag = 0; 0310 end 0311 output = results; 0312 0313 %% check for empty results (in case optimization failed) 0314 if ~isfield(results, 'x') || isempty(results.x) 0315 x = NaN(nx, 1); 0316 lam.lower = NaN(nx, 1); 0317 lam.upper = NaN(nx, 1); 0318 else 0319 x = results.x; 0320 lam.lower = zeros(nx, 1); 0321 lam.upper = zeros(nx, 1); 0322 end 0323 if ~isfield(results, 'objval') || isempty(results.objval) 0324 f = NaN; 0325 else 0326 f = results.objval; 0327 end 0328 if ~isfield(results, 'pi') || isempty(results.pi) 0329 pi = NaN(length(m.rhs), 1); 0330 else 0331 pi = results.pi; 0332 end 0333 if ~isfield(results, 'rc') || isempty(results.rc) 0334 rc = NaN(nx, 1); 0335 else 0336 rc = results.rc; 0337 end 0338 0339 kl = find(rc > 0); %% lower bound binding 0340 ku = find(rc < 0); %% upper bound binding 0341 lam.lower(kl) = rc(kl); 0342 lam.upper(ku) = -rc(ku); 0343 lam.eqlin = pi(1:neq); 0344 lam.ineqlin = pi(neq+(1:niq)); 0345 mu_l = zeros(nA, 1); 0346 mu_u = zeros(nA, 1); 0347 0348 %% repackage lambdas 0349 kl = find(lam.eqlin > 0); %% lower bound binding 0350 ku = find(lam.eqlin < 0); %% upper bound binding 0351 0352 mu_l(ieq(kl)) = lam.eqlin(kl); 0353 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0354 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0355 0356 mu_u(ieq(ku)) = -lam.eqlin(ku); 0357 mu_u(ilt) = -lam.ineqlin(1:nlt); 0358 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0359 0360 lambda = struct( ... 0361 'mu_l', mu_l, ... 0362 'mu_u', mu_u, ... 0363 'lower', lam.lower, ... 0364 'upper', lam.upper ... 0365 ); 0366 0367 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0368 if verbose 0369 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0370 end 0371 if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol) 0372 tol = opt.price_stage_warn_tol; 0373 else 0374 tol = 1e-7; 0375 end 0376 if length(vtype) == 1 0377 if vtype == 'I' || vtype == 'B' || vtype == 'N' 0378 k = (1:nx); 0379 elseif vtype == 'S' 0380 k = find(x == 0); 0381 end 0382 else %% length(vtype) == nx 0383 k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ... 0384 (vtype == 'S' & x' == 0)); 0385 end 0386 0387 x(k) = round(x(k)); 0388 xmin(k) = x(k); 0389 xmax(k) = x(k); 0390 x0 = x; 0391 % opt.grb_opt.Method = 0; %% primal simplex 0392 0393 [x_, f_, eflag_, output_, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt); 0394 if eflag ~= eflag_ 0395 error('miqps_gurobi: EXITFLAG from price computation stage = %d', eflag_); 0396 end 0397 if abs(f - f_)/max(abs(f), 1) > tol 0398 warning('miqps_gurobi: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1)); 0399 end 0400 xn = x; 0401 xn(abs(xn)<1) = 1; 0402 [mx, k] = max(abs(x - x_) ./ xn); 0403 if mx > tol 0404 warning('miqps_gurobi: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0405 end 0406 output.price_stage = output_; 0407 end