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_MEX 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
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_MEX exit flag
1 = converged
0 or negative values = negative of GUROBI_MEX exit flag
(see GUROBI_MEX documentation for details)
OUTPUT : GUROBI_MEX output struct
(see GUROBI_MEX 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_MEX.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_MEX 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 % grb_opt - options struct for GUROBI, value in 0034 % verbose overrides these options 0035 % PROBLEM : The inputs can alternatively be supplied in a single 0036 % PROBLEM struct with fields corresponding to the input arguments 0037 % described above: H, c, A, l, u, xmin, xmax, x0, opt 0038 % 0039 % Outputs: 0040 % X : solution vector 0041 % F : final objective function value 0042 % EXITFLAG : GUROBI_MEX exit flag 0043 % 1 = converged 0044 % 0 or negative values = negative of GUROBI_MEX exit flag 0045 % (see GUROBI_MEX documentation for details) 0046 % OUTPUT : GUROBI_MEX output struct 0047 % (see GUROBI_MEX documentation for details) 0048 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0049 % multipliers on the constraints, with fields: 0050 % mu_l - lower (left-hand) limit on linear constraints 0051 % mu_u - upper (right-hand) limit on linear constraints 0052 % lower - lower bound on optimization variables 0053 % upper - upper bound on optimization variables 0054 % 0055 % Note the calling syntax is almost identical to that of QUADPROG 0056 % from MathWorks' Optimization Toolbox. The main difference is that 0057 % the linear constraints are specified with A, L, U instead of 0058 % A, B, Aeq, Beq. 0059 % 0060 % Calling syntax options: 0061 % [x, f, exitflag, output, lambda] = ... 0062 % qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt) 0063 % 0064 % x = qps_gurobi(H, c, A, l, u) 0065 % x = qps_gurobi(H, c, A, l, u, xmin, xmax) 0066 % x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0) 0067 % x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt) 0068 % x = qps_gurobi(problem), where problem is a struct with fields: 0069 % H, c, A, l, u, xmin, xmax, x0, opt 0070 % all fields except 'c', 'A' and 'l' or 'u' are optional 0071 % x = qps_gurobi(...) 0072 % [x, f] = qps_gurobi(...) 0073 % [x, f, exitflag] = qps_gurobi(...) 0074 % [x, f, exitflag, output] = qps_gurobi(...) 0075 % [x, f, exitflag, output, lambda] = qps_gurobi(...) 0076 % 0077 % 0078 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0079 % H = [ 1003.1 4.3 6.3 5.9; 0080 % 4.3 2.2 2.1 3.9; 0081 % 6.3 2.1 3.5 4.8; 0082 % 5.9 3.9 4.8 10 ]; 0083 % c = zeros(4,1); 0084 % A = [ 1 1 1 1; 0085 % 0.17 0.11 0.10 0.18 ]; 0086 % l = [1; 0.10]; 0087 % u = [1; Inf]; 0088 % xmin = zeros(4,1); 0089 % x0 = [1; 0; 0; 1]; 0090 % opt = struct('verbose', 2); 0091 % [x, f, s, out, lambda] = qps_gurobi(H, c, A, l, u, xmin, [], x0, opt); 0092 % 0093 % See also GUROBI_MEX. 0094 0095 % MATPOWER 0096 % $Id: qps_gurobi.m,v 1.3 2011/09/09 15:27:52 cvs Exp $ 0097 % by Ray Zimmerman, PSERC Cornell 0098 % Copyright (c) 2010-2011 by Power System Engineering Research Center (PSERC) 0099 % 0100 % This file is part of MATPOWER. 0101 % See http://www.pserc.cornell.edu/matpower/ for more info. 0102 % 0103 % MATPOWER is free software: you can redistribute it and/or modify 0104 % it under the terms of the GNU General Public License as published 0105 % by the Free Software Foundation, either version 3 of the License, 0106 % or (at your option) any later version. 0107 % 0108 % MATPOWER is distributed in the hope that it will be useful, 0109 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0110 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0111 % GNU General Public License for more details. 0112 % 0113 % You should have received a copy of the GNU General Public License 0114 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0115 % 0116 % Additional permission under GNU GPL version 3 section 7 0117 % 0118 % If you modify MATPOWER, or any covered work, to interface with 0119 % other modules (such as MATLAB code and MEX-files) available in a 0120 % MATLAB(R) or comparable environment containing parts covered 0121 % under other licensing terms, the licensors of MATPOWER grant 0122 % you additional permission to convey the resulting work. 0123 0124 %%----- input argument handling ----- 0125 %% gather inputs 0126 if nargin == 1 && isstruct(H) %% problem struct 0127 p = H; 0128 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0129 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0130 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0131 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0132 if isfield(p, 'u'), u = p.u; else, u = []; end 0133 if isfield(p, 'l'), l = p.l; else, l = []; end 0134 if isfield(p, 'A'), A = p.A; else, A = []; end 0135 if isfield(p, 'c'), c = p.c; else, c = []; end 0136 if isfield(p, 'H'), H = p.H; else, H = []; end 0137 else %% individual args 0138 if nargin < 9 0139 opt = []; 0140 if nargin < 8 0141 x0 = []; 0142 if nargin < 7 0143 xmax = []; 0144 if nargin < 6 0145 xmin = []; 0146 end 0147 end 0148 end 0149 end 0150 end 0151 0152 %% define nx, set default values for missing optional inputs 0153 if isempty(H) || ~any(any(H)) 0154 if isempty(A) && isempty(xmin) && isempty(xmax) 0155 error('qps_gurobi: LP problem must include constraints or variable bounds'); 0156 else 0157 if ~isempty(A) 0158 nx = size(A, 2); 0159 elseif ~isempty(xmin) 0160 nx = length(xmin); 0161 else % if ~isempty(xmax) 0162 nx = length(xmax); 0163 end 0164 end 0165 else 0166 nx = size(H, 1); 0167 end 0168 if isempty(c) 0169 c = zeros(nx, 1); 0170 end 0171 if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0172 (isempty(u) || all(u == Inf)) 0173 A = sparse(0,nx); %% no limits => no linear constraints 0174 end 0175 nA = size(A, 1); %% number of original linear constraints 0176 if isempty(u) %% By default, linear inequalities are ... 0177 u = Inf * ones(nA, 1); %% ... unbounded above and ... 0178 end 0179 if isempty(l) 0180 l = -Inf * ones(nA, 1); %% ... unbounded below. 0181 end 0182 if isempty(xmin) %% By default, optimization variables are ... 0183 xmin = -Inf * ones(nx, 1); %% ... unbounded below and ... 0184 end 0185 if isempty(xmax) 0186 xmax = Inf * ones(nx, 1); %% ... unbounded above. 0187 end 0188 if isempty(x0) 0189 x0 = zeros(nx, 1); 0190 end 0191 0192 %% default options 0193 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0194 verbose = opt.verbose; 0195 else 0196 verbose = 0; 0197 end 0198 % if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it) 0199 % max_it = opt.max_it; 0200 % else 0201 % max_it = 0; 0202 % end 0203 0204 %% set up options struct for Gurobi 0205 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt) 0206 g_opt = gurobi_options(opt.grb_opt); 0207 else 0208 g_opt = gurobi_options; 0209 end 0210 g_opt.Display = min(verbose, 3); 0211 if verbose 0212 g_opt.DisplayInterval = 1; 0213 else 0214 g_opt.DisplayInterval = Inf; 0215 end 0216 0217 if ~issparse(A) 0218 A = sparse(A); 0219 end 0220 0221 %% split up linear constraints 0222 ieq = find( abs(u-l) <= eps ); %% equality 0223 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0224 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0225 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0226 0227 %% grab some dimensions 0228 nlt = length(ilt); %% number of upper bounded linear inequalities 0229 ngt = length(igt); %% number of lower bounded linear inequalities 0230 nbx = length(ibx); %% number of doubly bounded linear inequalities 0231 neq = length(ieq); %% number of equalities 0232 niq = nlt+ngt+2*nbx; %% number of inequalities 0233 0234 AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0235 bb = [ u(ieq); u(ilt); -l(igt); u(ibx); -l(ibx) ]; 0236 contypes = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]); 0237 0238 %% call the solver 0239 if isempty(H) || ~any(any(H)) 0240 lpqp = 'LP'; 0241 else 0242 lpqp = 'QP'; 0243 [rr, cc, vv] = find(H); 0244 g_opt.QP.qrow = int32(rr' - 1); 0245 g_opt.QP.qcol = int32(cc' - 1); 0246 g_opt.QP.qval = 0.5 * vv'; 0247 end 0248 if verbose 0249 methods = { 0250 'primal simplex', 0251 'dual simplex', 0252 'interior point', 0253 'concurrent', 0254 'deterministic concurrent' 0255 }; 0256 fprintf('Gurobi Version %s -- %s %s solver\n', ... 0257 '<unknown>', methods{g_opt.Method+1}, lpqp); 0258 end 0259 [x, f, eflag, output, lambda] = ... 0260 gurobi_mex(c', 1, AA, bb, contypes, xmin, xmax, 'C', g_opt); 0261 pi = lambda.Pi; 0262 rc = lambda.RC; 0263 output.flag = eflag; 0264 if eflag == 2 0265 eflag = 1; %% optimal solution found 0266 else 0267 eflag = -eflag; %% failed somehow 0268 end 0269 0270 %% check for empty results (in case optimization failed) 0271 if isempty(x) 0272 x = NaN(nx, 1); 0273 lam.lower = NaN(nx, 1); 0274 lam.upper = NaN(nx, 1); 0275 else 0276 lam.lower = zeros(nx, 1); 0277 lam.upper = zeros(nx, 1); 0278 end 0279 if isempty(f) 0280 f = NaN; 0281 end 0282 if isempty(pi) 0283 pi = NaN(length(bb), 1); 0284 end 0285 0286 kl = find(rc > 0); %% lower bound binding 0287 ku = find(rc < 0); %% upper bound binding 0288 lam.lower(kl) = rc(kl); 0289 lam.upper(ku) = -rc(ku); 0290 lam.eqlin = pi(1:neq); 0291 lam.ineqlin = pi(neq+(1:niq)); 0292 mu_l = zeros(nA, 1); 0293 mu_u = zeros(nA, 1); 0294 0295 %% repackage lambdas 0296 kl = find(lam.eqlin > 0); %% lower bound binding 0297 ku = find(lam.eqlin < 0); %% upper bound binding 0298 0299 mu_l(ieq(kl)) = lam.eqlin(kl); 0300 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt)); 0301 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0302 0303 mu_u(ieq(ku)) = -lam.eqlin(ku); 0304 mu_u(ilt) = -lam.ineqlin(1:nlt); 0305 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx)); 0306 0307 lambda = struct( ... 0308 'mu_l', mu_l, ... 0309 'mu_u', mu_u, ... 0310 'lower', lam.lower, ... 0311 'upper', lam.upper ... 0312 );