Home > matpower5.1 > 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 %   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 );

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005