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
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 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_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 % grb_opt - options struct for GUROBI, value in 0043 % verbose overrides these options 0044 % PROBLEM : The inputs can alternatively be supplied in a single 0045 % PROBLEM struct with fields corresponding to the input arguments 0046 % described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt 0047 % 0048 % Outputs: 0049 % X : solution vector 0050 % F : final objective function value 0051 % EXITFLAG : GUROBI exit flag 0052 % 1 = converged 0053 % 0 or negative values = negative of GUROBI exit flag 0054 % (see GUROBI documentation for details) 0055 % OUTPUT : GUROBI output struct 0056 % (see GUROBI 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_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0072 % 0073 % x = miqps_gurobi(H, c, A, l, u) 0074 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax) 0075 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0) 0076 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype) 0077 % x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt) 0078 % x = miqps_gurobi(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_gurobi(...) 0082 % [x, f] = miqps_gurobi(...) 0083 % [x, f, exitflag] = miqps_gurobi(...) 0084 % [x, f, exitflag, output] = miqps_gurobi(...) 0085 % [x, f, exitflag, output, lambda] = miqps_gurobi(...) 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_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt); 0102 % 0103 % See also GUROBI. 0104 0105 % MATPOWER 0106 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0107 % by Ray Zimmerman, PSERC Cornell 0108 % 0109 % $Id: miqps_gurobi.m 2661 2015-03-20 17:02:46Z ray $ 0110 % 0111 % This file is part of MATPOWER. 0112 % Covered by the 3-clause BSD License (see LICENSE file for details). 0113 % See http://www.pserc.cornell.edu/matpower/ for more info. 0114 0115 %%----- input argument handling ----- 0116 %% gather inputs 0117 if nargin == 1 && isstruct(H) %% problem struct 0118 p = H; 0119 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0120 if isfield(p, 'vtype'), vtype = p.vtype;else, vtype = []; end 0121 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0122 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0123 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0124 if isfield(p, 'u'), u = p.u; else, u = []; end 0125 if isfield(p, 'l'), l = p.l; else, l = []; end 0126 if isfield(p, 'A'), A = p.A; else, A = []; end 0127 if isfield(p, 'c'), c = p.c; else, c = []; end 0128 if isfield(p, 'H'), H = p.H; else, H = []; end 0129 else %% individual args 0130 if nargin < 10 0131 opt = []; 0132 if nargin < 9 0133 vtype = []; 0134 if nargin < 8 0135 x0 = []; 0136 if nargin < 7 0137 xmax = []; 0138 if nargin < 6 0139 xmin = []; 0140 end 0141 end 0142 end 0143 end 0144 end 0145 end 0146 0147 %% define nx, set default values for missing optional inputs 0148 if isempty(H) || ~any(any(H)) 0149 if isempty(A) && isempty(xmin) && isempty(xmax) 0150 error('miqps_gurobi: LP problem must include constraints or variable bounds'); 0151 else 0152 if ~isempty(A) 0153 nx = size(A, 2); 0154 elseif ~isempty(xmin) 0155 nx = length(xmin); 0156 else % if ~isempty(xmax) 0157 nx = length(xmax); 0158 end 0159 end 0160 else 0161 nx = size(H, 1); 0162 end 0163 if isempty(c) 0164 c = zeros(nx, 1); 0165 end 0166 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0167 (isempty(u) || all(u == Inf))) 0168 A = sparse(0,nx); %% no limits => no linear constraints 0169 end 0170 nA = size(A, 1); %% number of original linear constraints 0171 if isempty(u) %% By default, linear inequalities are ... 0172 u = Inf(nA, 1); %% ... unbounded above and ... 0173 end 0174 if isempty(l) 0175 l = -Inf(nA, 1); %% ... unbounded below. 0176 end 0177 if isempty(xmin) %% By default, optimization variables are ... 0178 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0179 end 0180 if isempty(xmax) 0181 xmax = Inf(nx, 1); %% ... unbounded above. 0182 end 0183 if isempty(x0) 0184 x0 = zeros(nx, 1); 0185 end 0186 0187 %% default options 0188 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0189 verbose = opt.verbose; 0190 else 0191 verbose = 0; 0192 end 0193 0194 %% set up options struct for Gurobi 0195 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0196 g_opt = gurobi_options(opt.grb_opt); 0197 else 0198 g_opt = gurobi_options; 0199 end 0200 if verbose > 1 0201 g_opt.LogToConsole = 1; 0202 g_opt.OutputFlag = 1; 0203 if verbose > 2 0204 g_opt.DisplayInterval = 1; 0205 else 0206 g_opt.DisplayInterval = 100; 0207 end 0208 else 0209 g_opt.LogToConsole = 0; 0210 g_opt.OutputFlag = 0; 0211 end 0212 0213 if ~issparse(A) 0214 A = sparse(A); 0215 end 0216 if issparse(c); 0217 c = full(c); 0218 end 0219 0220 %% split up linear constraints 0221 ieq = find( abs(u-l) <= eps ); %% equality 0222 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0223 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0224 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0225 0226 %% grab some dimensions 0227 nlt = length(ilt); %% number of upper bounded linear inequalities 0228 ngt = length(igt); %% number of lower bounded linear inequalities 0229 nbx = length(ibx); %% number of doubly bounded linear inequalities 0230 neq = length(ieq); %% number of equalities 0231 niq = nlt+ngt+2*nbx; %% number of inequalities 0232 0233 %% set up model 0234 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0235 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0236 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0237 m.lb = xmin; 0238 m.ub = xmax; 0239 m.obj = c'; 0240 if ~isempty(vtype) 0241 m.vtype = vtype; 0242 end 0243 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ... 0244 vtype == 'S' | vtype == 'N')) 0245 mi = 0; 0246 else 0247 mi = 1; 0248 end 0249 0250 %% call the solver 0251 if isempty(H) || ~any(any(H)) 0252 lpqp = 'LP'; 0253 else 0254 lpqp = 'QP'; 0255 if ~issparse(H) 0256 H = sparse(H); 0257 end 0258 m.Q = 0.5 * H; 0259 end 0260 if mi 0261 lpqp = ['MI' lpqp]; 0262 end 0263 if verbose 0264 methods = { 0265 'automatic', 0266 'primal simplex', 0267 'dual simplex', 0268 'interior point', 0269 'concurrent', 0270 'deterministic concurrent' 0271 }; 0272 vn = gurobiver; 0273 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0274 vn, methods{g_opt.Method+2}, lpqp); 0275 end 0276 results = gurobi(m, g_opt); 0277 switch results.status 0278 case 'LOADED', %% 1 0279 eflag = -1; 0280 case 'OPTIMAL', %% 2, optimal solution found 0281 eflag = 1; 0282 case 'INFEASIBLE', %% 3 0283 eflag = -3; 0284 case 'INF_OR_UNBD', %% 4 0285 eflag = -4; 0286 case 'UNBOUNDED', %% 5 0287 eflag = -5; 0288 case 'CUTOFF', %% 6 0289 eflag = -6; 0290 case 'ITERATION_LIMIT', %% 7 0291 eflag = -7; 0292 case 'NODE_LIMIT', %% 8 0293 eflag = -8; 0294 case 'TIME_LIMIT', %% 9 0295 eflag = -9; 0296 case 'SOLUTION_LIMIT', %% 10 0297 eflag = -10; 0298 case 'INTERRUPTED', %% 11 0299 eflag = -11; 0300 case 'NUMERIC', %% 12 0301 eflag = -12; 0302 case 'SUBOPTIMAL', %% 13 0303 eflag = -13; 0304 case 'INPROGRESS', %% 14 0305 eflag = -14; 0306 otherwise, 0307 eflag = 0; 0308 end 0309 output = results; 0310 0311 %% check for empty results (in case optimization failed) 0312 if ~isfield(results, 'x') || isempty(results.x) 0313 x = NaN(nx, 1); 0314 lam.lower = NaN(nx, 1); 0315 lam.upper = NaN(nx, 1); 0316 else 0317 x = results.x; 0318 lam.lower = zeros(nx, 1); 0319 lam.upper = zeros(nx, 1); 0320 end 0321 if ~isfield(results, 'objval') || isempty(results.objval) 0322 f = NaN; 0323 else 0324 f = results.objval; 0325 end 0326 if ~isfield(results, 'pi') || isempty(results.pi) 0327 pi = NaN(length(m.rhs), 1); 0328 else 0329 pi = results.pi; 0330 end 0331 if ~isfield(results, 'rc') || isempty(results.rc) 0332 rc = NaN(nx, 1); 0333 else 0334 rc = results.rc; 0335 end 0336 0337 kl = find(rc > 0); %% lower bound binding 0338 ku = find(rc < 0); %% upper bound binding 0339 lam.lower(kl) = rc(kl); 0340 lam.upper(ku) = -rc(ku); 0341 lam.eqlin = pi(1:neq); 0342 lam.ineqlin = pi(neq+(1:niq)); 0343 mu_l = zeros(nA, 1); 0344 mu_u = zeros(nA, 1); 0345 0346 %% repackage lambdas 0347 kl = find(lam.eqlin > 0); %% lower bound binding 0348 ku = find(lam.eqlin < 0); %% upper bound binding 0349 0350 mu_l(ieq(kl)) = lam.eqlin(kl); 0351 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0352 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0353 0354 mu_u(ieq(ku)) = -lam.eqlin(ku); 0355 mu_u(ilt) = -lam.ineqlin(1:nlt); 0356 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0357 0358 lambda = struct( ... 0359 'mu_l', mu_l, ... 0360 'mu_u', mu_u, ... 0361 'lower', lam.lower, ... 0362 'upper', lam.upper ... 0363 ); 0364 0365 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices) 0366 if verbose 0367 fprintf('--- Integer stage complete, starting price computation stage ---\n'); 0368 end 0369 tol = 1e-7; 0370 if length(vtype) == 1 0371 if vtype == 'I' || vtype == 'B' || vtype == 'N' 0372 k = (1:nx); 0373 elseif vtype == 'S' 0374 k = find(x == 0); 0375 end 0376 else %% length(vtype) == nx 0377 k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ... 0378 (vtype == 'S' & x' == 0)); 0379 end 0380 0381 x(k) = round(x(k)); 0382 xmin(k) = x(k); 0383 xmax(k) = x(k); 0384 x0 = x; 0385 % opt.grb_opt.Method = 0; %% primal simplex 0386 0387 [x_, f_, eflag_, output_, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt); 0388 if eflag ~= eflag_ 0389 error('miqps_gurobi: EXITFLAG from price computation stage = %d', eflag_); 0390 end 0391 if abs(f - f_)/max(abs(f), 1) > tol 0392 warning('miqps_gurobi: 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_gurobi: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k)); 0399 end 0400 output.price_stage = output_; 0401 end