Home > matpower5.0 > qps_gurobi.m

qps_gurobi

PURPOSE ^

QPS_GUROBI Quadratic Program Solver based on GUROBI.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 );

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005