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 % $Id: qps_gurobi.m 2294 2014-03-10 19:00:05Z ray $ 0098 % by Ray Zimmerman, PSERC Cornell 0099 % Copyright (c) 2010-2013 by Power System Engineering Research Center (PSERC) 0100 % 0101 % This file is part of MATPOWER. 0102 % See http://www.pserc.cornell.edu/matpower/ for more info. 0103 % 0104 % MATPOWER is free software: you can redistribute it and/or modify 0105 % it under the terms of the GNU General Public License as published 0106 % by the Free Software Foundation, either version 3 of the License, 0107 % or (at your option) any later version. 0108 % 0109 % MATPOWER is distributed in the hope that it will be useful, 0110 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0111 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0112 % GNU General Public License for more details. 0113 % 0114 % You should have received a copy of the GNU General Public License 0115 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0116 % 0117 % Additional permission under GNU GPL version 3 section 7 0118 % 0119 % If you modify MATPOWER, or any covered work, to interface with 0120 % other modules (such as MATLAB code and MEX-files) available in a 0121 % MATLAB(R) or comparable environment containing parts covered 0122 % under other licensing terms, the licensors of MATPOWER grant 0123 % you additional permission to convey the resulting work. 0124 0125 %%----- input argument handling ----- 0126 %% gather inputs 0127 if nargin == 1 && isstruct(H) %% problem struct 0128 p = H; 0129 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0130 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0131 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0132 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0133 if isfield(p, 'u'), u = p.u; else, u = []; end 0134 if isfield(p, 'l'), l = p.l; else, l = []; end 0135 if isfield(p, 'A'), A = p.A; else, A = []; end 0136 if isfield(p, 'c'), c = p.c; else, c = []; end 0137 if isfield(p, 'H'), H = p.H; else, H = []; end 0138 else %% individual args 0139 if nargin < 9 0140 opt = []; 0141 if nargin < 8 0142 x0 = []; 0143 if nargin < 7 0144 xmax = []; 0145 if nargin < 6 0146 xmin = []; 0147 end 0148 end 0149 end 0150 end 0151 end 0152 0153 %% define nx, set default values for missing optional inputs 0154 if isempty(H) || ~any(any(H)) 0155 if isempty(A) && isempty(xmin) && isempty(xmax) 0156 error('qps_gurobi: LP problem must include constraints or variable bounds'); 0157 else 0158 if ~isempty(A) 0159 nx = size(A, 2); 0160 elseif ~isempty(xmin) 0161 nx = length(xmin); 0162 else % if ~isempty(xmax) 0163 nx = length(xmax); 0164 end 0165 end 0166 else 0167 nx = size(H, 1); 0168 end 0169 if isempty(c) 0170 c = zeros(nx, 1); 0171 end 0172 if isempty(A) || ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0173 (isempty(u) || all(u == Inf)) 0174 A = sparse(0,nx); %% no limits => no linear constraints 0175 end 0176 nA = size(A, 1); %% number of original linear constraints 0177 if isempty(u) %% By default, linear inequalities are ... 0178 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0179 end 0180 if isempty(l) 0181 l = -Inf * ones(nA, 1); %% ... unbounded below. 0182 end 0183 if isempty(xmin) %% By default, optimization variables are ... 0184 xmin = -Inf * ones(nx, 1); %% ... unbounded below and ... 0185 end 0186 if isempty(xmax) 0187 xmax = Inf * ones(nx, 1); %% ... unbounded above. 0188 end 0189 if isempty(x0) 0190 x0 = zeros(nx, 1); 0191 end 0192 0193 %% default options 0194 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0195 verbose = opt.verbose; 0196 else 0197 verbose = 0; 0198 end 0199 0200 %% set up options struct for Gurobi 0201 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0202 g_opt = gurobi_options(opt.grb_opt); 0203 else 0204 g_opt = gurobi_options; 0205 end 0206 if verbose > 1 0207 g_opt.LogToConsole = 1; 0208 g_opt.OutputFlag = 1; 0209 if verbose > 2 0210 g_opt.DisplayInterval = 1; 0211 else 0212 g_opt.DisplayInterval = 100; 0213 end 0214 else 0215 g_opt.LogToConsole = 0; 0216 g_opt.OutputFlag = 0; 0217 end 0218 0219 if ~issparse(A) 0220 A = sparse(A); 0221 end 0222 if issparse(c); 0223 c = full(c); 0224 end 0225 0226 %% split up linear constraints 0227 ieq = find( abs(u-l) <= eps ); %% equality 0228 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0229 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0230 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0231 0232 %% grab some dimensions 0233 nlt = length(ilt); %% number of upper bounded linear inequalities 0234 ngt = length(igt); %% number of lower bounded linear inequalities 0235 nbx = length(ibx); %% number of doubly bounded linear inequalities 0236 neq = length(ieq); %% number of equalities 0237 niq = nlt+ngt+2*nbx; %% number of inequalities 0238 0239 %% set up model 0240 m.A = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0241 m.rhs = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0242 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0243 m.lb = xmin; 0244 m.ub = xmax; 0245 m.obj = c'; 0246 0247 %% call the solver 0248 if isempty(H) || ~any(any(H)) 0249 lpqp = 'LP'; 0250 else 0251 lpqp = 'QP'; 0252 if ~issparse(H) 0253 H = sparse(H); 0254 end 0255 m.Q = 0.5 * H; 0256 end 0257 if verbose 0258 methods = { 0259 'automatic', 0260 'primal simplex', 0261 'dual simplex', 0262 'interior point', 0263 'concurrent', 0264 'deterministic concurrent' 0265 }; 0266 vn = gurobiver; 0267 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0268 vn, methods{g_opt.Method+2}, lpqp); 0269 end 0270 results = gurobi(m, g_opt); 0271 switch results.status 0272 case 'LOADED', %% 1 0273 eflag = -1; 0274 case 'OPTIMAL', %% 2, optimal solution found 0275 eflag = 1; 0276 case 'INFEASIBLE', %% 3 0277 eflag = -3; 0278 case 'INF_OR_UNBD', %% 4 0279 eflag = -4; 0280 case 'UNBOUNDED', %% 5 0281 eflag = -5; 0282 case 'CUTOFF', %% 6 0283 eflag = -6; 0284 case 'ITERATION_LIMIT', %% 7 0285 eflag = -7; 0286 case 'NODE_LIMIT', %% 8 0287 eflag = -8; 0288 case 'TIME_LIMIT', %% 9 0289 eflag = -9; 0290 case 'SOLUTION_LIMIT', %% 10 0291 eflag = -10; 0292 case 'INTERRUPTED', %% 11 0293 eflag = -11; 0294 case 'NUMERIC', %% 12 0295 eflag = -12; 0296 case 'SUBOPTIMAL', %% 13 0297 eflag = -13; 0298 case 'INPROGRESS', %% 14 0299 eflag = -14; 0300 otherwise, 0301 eflag = 0; 0302 end 0303 output = results; 0304 0305 %% check for empty results (in case optimization failed) 0306 if ~isfield(results, 'x') || isempty(results.x) 0307 x = NaN(nx, 1); 0308 lam.lower = NaN(nx, 1); 0309 lam.upper = NaN(nx, 1); 0310 else 0311 x = results.x; 0312 lam.lower = zeros(nx, 1); 0313 lam.upper = zeros(nx, 1); 0314 end 0315 if ~isfield(results, 'objval') || isempty(results.objval) 0316 f = NaN; 0317 else 0318 f = results.objval; 0319 end 0320 if ~isfield(results, 'pi') || isempty(results.pi) 0321 pi = NaN(length(m.rhs), 1); 0322 else 0323 pi = results.pi; 0324 end 0325 if ~isfield(results, 'rc') || isempty(results.rc) 0326 rc = NaN(nx, 1); 0327 else 0328 rc = results.rc; 0329 end 0330 0331 kl = find(rc > 0); %% lower bound binding 0332 ku = find(rc < 0); %% upper bound binding 0333 lam.lower(kl) = rc(kl); 0334 lam.upper(ku) = -rc(ku); 0335 lam.eqlin = pi(1:neq); 0336 lam.ineqlin = pi(neq+(1:niq)); 0337 mu_l = zeros(nA, 1); 0338 mu_u = zeros(nA, 1); 0339 0340 %% repackage lambdas 0341 kl = find(lam.eqlin > 0); %% lower bound binding 0342 ku = find(lam.eqlin < 0); %% upper bound binding 0343 0344 mu_l(ieq(kl)) = lam.eqlin(kl); 0345 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0346 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0347 0348 mu_u(ieq(ku)) = -lam.eqlin(ku); 0349 mu_u(ilt) = -lam.ineqlin(1:nlt); 0350 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0351 0352 lambda = struct( ... 0353 'mu_l', mu_l, ... 0354 'mu_u', mu_u, ... 0355 'lower', lam.lower, ... 0356 'upper', lam.upper ... 0357 );