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