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

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005