QPS_GUROBI Quadratic Program Solver based on GUROBI.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, OPT)
[X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_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
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
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, 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] = ...
qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_gurobi(H, c, A, l, u)
x = qps_gurobi(H, c, A, l, u, xmin, xmax)
x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0)
x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_gurobi(problem), where problem is a struct with fields:
H, c, A, l, u, xmin, xmax, x0, opt
all fields except 'c', 'A' and 'l' or 'u' are optional
x = qps_gurobi(...)
[x, f] = qps_gurobi(...)
[x, f, exitflag] = qps_gurobi(...)
[x, f, exitflag, output] = qps_gurobi(...)
[x, f, exitflag, output, lambda] = qps_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] = qps_gurobi(H, c, A, l, u, xmin, [], x0, opt);
See also QPS_MASTER, GUROBI.0001 function [x, f, eflag, output, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_GUROBI Quadratic Program Solver based on GUROBI. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_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 % OPT : optional options structure with the following fields, 0028 % all of which are also optional (default values shown in 0029 % parentheses) 0030 % verbose (0) - controls level of progress output displayed 0031 % 0 = no progress output 0032 % 1 = some progress output 0033 % 2 = verbose progress output 0034 % 3 = even more verbose progress output 0035 % grb_opt - options struct for GUROBI, value in verbose 0036 % overrides these options 0037 % PROBLEM : The inputs can alternatively be supplied in a single 0038 % PROBLEM struct with fields corresponding to the input arguments 0039 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0040 % 0041 % Outputs: 0042 % X : solution vector 0043 % F : final objective function value 0044 % EXITFLAG : GUROBI exit flag 0045 % 1 = converged 0046 % 0 or negative values = negative of GUROBI exit flag 0047 % (see GUROBI documentation for details) 0048 % OUTPUT : GUROBI output struct 0049 % (see GUROBI documentation for details) 0050 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0051 % multipliers on the constraints, with fields: 0052 % mu_l - lower (left-hand) limit on linear constraints 0053 % mu_u - upper (right-hand) limit on linear constraints 0054 % lower - lower bound on optimization variables 0055 % upper - upper bound on optimization variables 0056 % 0057 % Note the calling syntax is almost identical to that of QUADPROG 0058 % from MathWorks' Optimization Toolbox. The main difference is that 0059 % the linear constraints are specified with A, L, U instead of 0060 % A, B, Aeq, Beq. 0061 % 0062 % Calling syntax options: 0063 % [x, f, exitflag, output, lambda] = ... 0064 % qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt) 0065 % 0066 % x = qps_gurobi(H, c, A, l, u) 0067 % x = qps_gurobi(H, c, A, l, u, xmin, xmax) 0068 % x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0) 0069 % x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt) 0070 % x = qps_gurobi(problem), where problem is a struct with fields: 0071 % H, c, A, l, u, xmin, xmax, x0, opt 0072 % all fields except 'c', 'A' and 'l' or 'u' are optional 0073 % x = qps_gurobi(...) 0074 % [x, f] = qps_gurobi(...) 0075 % [x, f, exitflag] = qps_gurobi(...) 0076 % [x, f, exitflag, output] = qps_gurobi(...) 0077 % [x, f, exitflag, output, lambda] = qps_gurobi(...) 0078 % 0079 % 0080 % Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm) 0081 % H = [ 1003.1 4.3 6.3 5.9; 0082 % 4.3 2.2 2.1 3.9; 0083 % 6.3 2.1 3.5 4.8; 0084 % 5.9 3.9 4.8 10 ]; 0085 % c = zeros(4,1); 0086 % A = [ 1 1 1 1; 0087 % 0.17 0.11 0.10 0.18 ]; 0088 % l = [1; 0.10]; 0089 % u = [1; Inf]; 0090 % xmin = zeros(4,1); 0091 % x0 = [1; 0; 0; 1]; 0092 % opt = struct('verbose', 2); 0093 % [x, f, s, out, lambda] = qps_gurobi(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also QPS_MASTER, GUROBI. 0096 0097 % MP-Opt-Model 0098 % Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % This file is part of MP-Opt-Model. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See https://github.com/MATPOWER/mp-opt-model for more info. 0104 0105 %%----- input argument handling ----- 0106 %% gather inputs 0107 if nargin == 1 && isstruct(H) %% problem struct 0108 p = H; 0109 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0110 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0111 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0112 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0113 if isfield(p, 'u'), u = p.u; else, u = []; end 0114 if isfield(p, 'l'), l = p.l; else, l = []; end 0115 if isfield(p, 'A'), A = p.A; else, A = []; end 0116 if isfield(p, 'c'), c = p.c; else, c = []; end 0117 if isfield(p, 'H'), H = p.H; else, H = []; end 0118 else %% individual args 0119 if nargin < 9 0120 opt = []; 0121 if nargin < 8 0122 x0 = []; 0123 if nargin < 7 0124 xmax = []; 0125 if nargin < 6 0126 xmin = []; 0127 end 0128 end 0129 end 0130 end 0131 end 0132 0133 %% define nx, set default values for missing optional inputs 0134 if isempty(H) || ~any(any(H)) 0135 if isempty(A) && isempty(xmin) && isempty(xmax) 0136 error('qps_gurobi: LP problem must include constraints or variable bounds'); 0137 else 0138 if ~isempty(A) 0139 nx = size(A, 2); 0140 elseif ~isempty(xmin) 0141 nx = length(xmin); 0142 else % if ~isempty(xmax) 0143 nx = length(xmax); 0144 end 0145 end 0146 else 0147 nx = size(H, 1); 0148 end 0149 if isempty(c) 0150 c = zeros(nx, 1); 0151 end 0152 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0153 (isempty(u) || all(u == Inf))) 0154 A = sparse(0,nx); %% no limits => no linear constraints 0155 end 0156 nA = size(A, 1); %% number of original linear constraints 0157 if isempty(u) %% By default, linear inequalities are ... 0158 u = Inf(nA, 1); %% ... unbounded above and ... 0159 end 0160 if isempty(l) 0161 l = -Inf(nA, 1); %% ... unbounded below. 0162 end 0163 if isempty(xmin) %% By default, optimization variables are ... 0164 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0165 end 0166 if isempty(xmax) 0167 xmax = Inf(nx, 1); %% ... unbounded above. 0168 end 0169 if isempty(x0) 0170 x0 = zeros(nx, 1); 0171 end 0172 0173 %% default options 0174 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0175 verbose = opt.verbose; 0176 else 0177 verbose = 0; 0178 end 0179 0180 %% set up options struct for Gurobi 0181 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0182 g_opt = gurobi_options(opt.grb_opt); 0183 else 0184 g_opt = gurobi_options; 0185 end 0186 if verbose > 1 0187 g_opt.LogToConsole = 1; 0188 g_opt.OutputFlag = 1; 0189 if verbose > 2 0190 g_opt.DisplayInterval = 1; 0191 else 0192 g_opt.DisplayInterval = 100; 0193 end 0194 else 0195 g_opt.LogToConsole = 0; 0196 g_opt.OutputFlag = 0; 0197 end 0198 0199 if ~issparse(A) 0200 A = sparse(A); 0201 end 0202 if issparse(c); 0203 c = full(c); 0204 end 0205 0206 %% split up linear constraints 0207 ieq = find( abs(u-l) <= eps ); %% equality 0208 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0209 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0210 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0211 0212 %% grab some dimensions 0213 nlt = length(ilt); %% number of upper bounded linear inequalities 0214 ngt = length(igt); %% number of lower bounded linear inequalities 0215 nbx = length(ibx); %% number of doubly bounded linear inequalities 0216 neq = length(ieq); %% number of equalities 0217 niq = nlt+ngt+2*nbx; %% number of inequalities 0218 0219 %% set up model 0220 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0221 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0222 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0223 m.lb = xmin; 0224 m.ub = xmax; 0225 m.obj = c'; 0226 0227 %% call the solver 0228 if isempty(H) || ~any(any(H)) 0229 lpqp = 'LP'; 0230 else 0231 lpqp = 'QP'; 0232 if ~issparse(H) 0233 H = sparse(H); 0234 end 0235 m.Q = 0.5 * H; 0236 end 0237 if verbose 0238 alg_names = { 0239 'automatic', 0240 'primal simplex', 0241 'dual simplex', 0242 'interior point', 0243 'concurrent', 0244 'deterministic concurrent' 0245 }; 0246 vn = gurobiver; 0247 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0248 vn, alg_names{g_opt.Method+2}, lpqp); 0249 end 0250 results = gurobi(m, g_opt); 0251 switch results.status 0252 case 'LOADED', %% 1 0253 eflag = -1; 0254 case 'OPTIMAL', %% 2, optimal solution found 0255 eflag = 1; 0256 case 'INFEASIBLE', %% 3 0257 eflag = -3; 0258 case 'INF_OR_UNBD', %% 4 0259 eflag = -4; 0260 case 'UNBOUNDED', %% 5 0261 eflag = -5; 0262 case 'CUTOFF', %% 6 0263 eflag = -6; 0264 case 'ITERATION_LIMIT', %% 7 0265 eflag = -7; 0266 case 'NODE_LIMIT', %% 8 0267 eflag = -8; 0268 case 'TIME_LIMIT', %% 9 0269 eflag = -9; 0270 case 'SOLUTION_LIMIT', %% 10 0271 eflag = -10; 0272 case 'INTERRUPTED', %% 11 0273 eflag = -11; 0274 case 'NUMERIC', %% 12 0275 eflag = -12; 0276 case 'SUBOPTIMAL', %% 13 0277 eflag = -13; 0278 case 'INPROGRESS', %% 14 0279 eflag = -14; 0280 otherwise, 0281 eflag = 0; 0282 end 0283 output = results; 0284 0285 %% check for empty results (in case optimization failed) 0286 if ~isfield(results, 'x') || isempty(results.x) 0287 x = NaN(nx, 1); 0288 lam.lower = NaN(nx, 1); 0289 lam.upper = NaN(nx, 1); 0290 else 0291 x = results.x; 0292 lam.lower = zeros(nx, 1); 0293 lam.upper = zeros(nx, 1); 0294 end 0295 if ~isfield(results, 'objval') || isempty(results.objval) 0296 f = NaN; 0297 else 0298 f = results.objval; 0299 end 0300 if ~isfield(results, 'pi') || isempty(results.pi) 0301 pi = NaN(length(m.rhs), 1); 0302 else 0303 pi = results.pi; 0304 end 0305 if ~isfield(results, 'rc') || isempty(results.rc) 0306 rc = NaN(nx, 1); 0307 else 0308 rc = results.rc; 0309 end 0310 0311 kl = find(rc > 0); %% lower bound binding 0312 ku = find(rc < 0); %% upper bound binding 0313 lam.lower(kl) = rc(kl); 0314 lam.upper(ku) = -rc(ku); 0315 lam.eqlin = pi(1:neq); 0316 lam.ineqlin = pi(neq+(1:niq)); 0317 mu_l = zeros(nA, 1); 0318 mu_u = zeros(nA, 1); 0319 0320 %% repackage lambdas 0321 kl = find(lam.eqlin > 0); %% lower bound binding 0322 ku = find(lam.eqlin < 0); %% upper bound binding 0323 0324 mu_l(ieq(kl)) = lam.eqlin(kl); 0325 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0326 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0327 0328 mu_u(ieq(ku)) = -lam.eqlin(ku); 0329 mu_u(ilt) = -lam.ineqlin(1:nlt); 0330 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0331 0332 lambda = struct( ... 0333 'mu_l', mu_l, ... 0334 'mu_u', mu_u, ... 0335 'lower', lam.lower, ... 0336 'upper', lam.upper ... 0337 );