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-2016, Power Systems Engineering Research Center (PSERC) 0098 % by Ray Zimmerman, PSERC Cornell 0099 % 0100 % This file is part of MATPOWER. 0101 % Covered by the 3-clause BSD License (see LICENSE file for details). 0102 % See http://www.pserc.cornell.edu/matpower/ for more info. 0103 0104 %%----- input argument handling ----- 0105 %% gather inputs 0106 if nargin == 1 && isstruct(H) %% problem struct 0107 p = H; 0108 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0109 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0110 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0111 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0112 if isfield(p, 'u'), u = p.u; else, u = []; end 0113 if isfield(p, 'l'), l = p.l; else, l = []; end 0114 if isfield(p, 'A'), A = p.A; else, A = []; end 0115 if isfield(p, 'c'), c = p.c; else, c = []; end 0116 if isfield(p, 'H'), H = p.H; else, H = []; end 0117 else %% individual args 0118 if nargin < 9 0119 opt = []; 0120 if nargin < 8 0121 x0 = []; 0122 if nargin < 7 0123 xmax = []; 0124 if nargin < 6 0125 xmin = []; 0126 end 0127 end 0128 end 0129 end 0130 end 0131 0132 %% define nx, set default values for missing optional inputs 0133 if isempty(H) || ~any(any(H)) 0134 if isempty(A) && isempty(xmin) && isempty(xmax) 0135 error('qps_gurobi: LP problem must include constraints or variable bounds'); 0136 else 0137 if ~isempty(A) 0138 nx = size(A, 2); 0139 elseif ~isempty(xmin) 0140 nx = length(xmin); 0141 else % if ~isempty(xmax) 0142 nx = length(xmax); 0143 end 0144 end 0145 else 0146 nx = size(H, 1); 0147 end 0148 if isempty(c) 0149 c = zeros(nx, 1); 0150 end 0151 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0152 (isempty(u) || all(u == Inf))) 0153 A = sparse(0,nx); %% no limits => no linear constraints 0154 end 0155 nA = size(A, 1); %% number of original linear constraints 0156 if isempty(u) %% By default, linear inequalities are ... 0157 u = Inf(nA, 1); %% ... unbounded above and ... 0158 end 0159 if isempty(l) 0160 l = -Inf(nA, 1); %% ... unbounded below. 0161 end 0162 if isempty(xmin) %% By default, optimization variables are ... 0163 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0164 end 0165 if isempty(xmax) 0166 xmax = Inf(nx, 1); %% ... unbounded above. 0167 end 0168 if isempty(x0) 0169 x0 = zeros(nx, 1); 0170 end 0171 0172 %% default options 0173 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0174 verbose = opt.verbose; 0175 else 0176 verbose = 0; 0177 end 0178 0179 %% set up options struct for Gurobi 0180 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0181 g_opt = gurobi_options(opt.grb_opt); 0182 else 0183 g_opt = gurobi_options; 0184 end 0185 if verbose > 1 0186 g_opt.LogToConsole = 1; 0187 g_opt.OutputFlag = 1; 0188 if verbose > 2 0189 g_opt.DisplayInterval = 1; 0190 else 0191 g_opt.DisplayInterval = 100; 0192 end 0193 else 0194 g_opt.LogToConsole = 0; 0195 g_opt.OutputFlag = 0; 0196 end 0197 0198 if ~issparse(A) 0199 A = sparse(A); 0200 end 0201 if issparse(c); 0202 c = full(c); 0203 end 0204 0205 %% split up linear constraints 0206 ieq = find( abs(u-l) <= eps ); %% equality 0207 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0208 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0209 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0210 0211 %% grab some dimensions 0212 nlt = length(ilt); %% number of upper bounded linear inequalities 0213 ngt = length(igt); %% number of lower bounded linear inequalities 0214 nbx = length(ibx); %% number of doubly bounded linear inequalities 0215 neq = length(ieq); %% number of equalities 0216 niq = nlt+ngt+2*nbx; %% number of inequalities 0217 0218 %% set up model 0219 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0220 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0221 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0222 m.lb = xmin; 0223 m.ub = xmax; 0224 m.obj = c'; 0225 0226 %% call the solver 0227 if isempty(H) || ~any(any(H)) 0228 lpqp = 'LP'; 0229 else 0230 lpqp = 'QP'; 0231 if ~issparse(H) 0232 H = sparse(H); 0233 end 0234 m.Q = 0.5 * H; 0235 end 0236 if verbose 0237 alg_names = { 0238 'automatic', 0239 'primal simplex', 0240 'dual simplex', 0241 'interior point', 0242 'concurrent', 0243 'deterministic concurrent' 0244 }; 0245 vn = gurobiver; 0246 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0247 vn, alg_names{g_opt.Method+2}, lpqp); 0248 end 0249 results = gurobi(m, g_opt); 0250 switch results.status 0251 case 'LOADED', %% 1 0252 eflag = -1; 0253 case 'OPTIMAL', %% 2, optimal solution found 0254 eflag = 1; 0255 case 'INFEASIBLE', %% 3 0256 eflag = -3; 0257 case 'INF_OR_UNBD', %% 4 0258 eflag = -4; 0259 case 'UNBOUNDED', %% 5 0260 eflag = -5; 0261 case 'CUTOFF', %% 6 0262 eflag = -6; 0263 case 'ITERATION_LIMIT', %% 7 0264 eflag = -7; 0265 case 'NODE_LIMIT', %% 8 0266 eflag = -8; 0267 case 'TIME_LIMIT', %% 9 0268 eflag = -9; 0269 case 'SOLUTION_LIMIT', %% 10 0270 eflag = -10; 0271 case 'INTERRUPTED', %% 11 0272 eflag = -11; 0273 case 'NUMERIC', %% 12 0274 eflag = -12; 0275 case 'SUBOPTIMAL', %% 13 0276 eflag = -13; 0277 case 'INPROGRESS', %% 14 0278 eflag = -14; 0279 otherwise, 0280 eflag = 0; 0281 end 0282 output = results; 0283 0284 %% check for empty results (in case optimization failed) 0285 if ~isfield(results, 'x') || isempty(results.x) 0286 x = NaN(nx, 1); 0287 lam.lower = NaN(nx, 1); 0288 lam.upper = NaN(nx, 1); 0289 else 0290 x = results.x; 0291 lam.lower = zeros(nx, 1); 0292 lam.upper = zeros(nx, 1); 0293 end 0294 if ~isfield(results, 'objval') || isempty(results.objval) 0295 f = NaN; 0296 else 0297 f = results.objval; 0298 end 0299 if ~isfield(results, 'pi') || isempty(results.pi) 0300 pi = NaN(length(m.rhs), 1); 0301 else 0302 pi = results.pi; 0303 end 0304 if ~isfield(results, 'rc') || isempty(results.rc) 0305 rc = NaN(nx, 1); 0306 else 0307 rc = results.rc; 0308 end 0309 0310 kl = find(rc > 0); %% lower bound binding 0311 ku = find(rc < 0); %% upper bound binding 0312 lam.lower(kl) = rc(kl); 0313 lam.upper(ku) = -rc(ku); 0314 lam.eqlin = pi(1:neq); 0315 lam.ineqlin = pi(neq+(1:niq)); 0316 mu_l = zeros(nA, 1); 0317 mu_u = zeros(nA, 1); 0318 0319 %% repackage lambdas 0320 kl = find(lam.eqlin > 0); %% lower bound binding 0321 ku = find(lam.eqlin < 0); %% upper bound binding 0322 0323 mu_l(ieq(kl)) = lam.eqlin(kl); 0324 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0325 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0326 0327 mu_u(ieq(ku)) = -lam.eqlin(ku); 0328 mu_u(ilt) = -lam.ineqlin(1:nlt); 0329 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0330 0331 lambda = struct( ... 0332 'mu_l', mu_l, ... 0333 'mu_u', mu_u, ... 0334 'lower', lam.lower, ... 0335 'upper', lam.upper ... 0336 );