Home > matpower7.1 > mp-opt-model > lib > 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)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_GUROBI(PROBLEM)
   A wrapper function providing a 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 https://v8doc.sas.com/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 QPS_MASTER, 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 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_GUROBI(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   GUROBI to solve the following QP (quadratic programming)
0008 %   problem:
0009 %
0010 %       min 1/2 X'*H*X + C'*X
0011 %        X
0012 %
0013 %   subject to
0014 %
0015 %       L <= A*X <= U       (linear constraints)
0016 %       XMIN <= X <= XMAX   (variable bounds)
0017 %
0018 %   Inputs (all optional except H, C, A and L):
0019 %       H : matrix (possibly sparse) of quadratic cost coefficients
0020 %       C : vector of linear cost coefficients
0021 %       A, L, U : define the optional linear constraints. Default
0022 %           values for the elements of L and U are -Inf and Inf,
0023 %           respectively.
0024 %       XMIN, XMAX : optional lower and upper bounds on the
0025 %           X variables, defaults are -Inf and Inf, respectively.
0026 %       X0 : optional starting value of optimization vector X
0027 %       OPT : optional options structure with the following fields,
0028 %           all of which are also optional (default values shown in
0029 %           parentheses)
0030 %           verbose (0) - controls level of progress output displayed
0031 %               0 = no progress output
0032 %               1 = some progress output
0033 %               2 = verbose progress output
0034 %               3 = even more verbose progress output
0035 %           grb_opt - options struct for GUROBI, value in verbose
0036 %                   overrides these options
0037 %       PROBLEM : The inputs can alternatively be supplied in a single
0038 %           PROBLEM struct with fields corresponding to the input arguments
0039 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0040 %
0041 %   Outputs:
0042 %       X : solution vector
0043 %       F : final objective function value
0044 %       EXITFLAG : GUROBI exit flag
0045 %           1 = converged
0046 %           0 or negative values = negative of GUROBI exit flag
0047 %           (see GUROBI documentation for details)
0048 %       OUTPUT : GUROBI output struct
0049 %           (see GUROBI documentation for details)
0050 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0051 %           multipliers on the constraints, with fields:
0052 %           mu_l - lower (left-hand) limit on linear constraints
0053 %           mu_u - upper (right-hand) limit on linear constraints
0054 %           lower - lower bound on optimization variables
0055 %           upper - upper bound on optimization variables
0056 %
0057 %   Note the calling syntax is almost identical to that of QUADPROG
0058 %   from MathWorks' Optimization Toolbox. The main difference is that
0059 %   the linear constraints are specified with A, L, U instead of
0060 %   A, B, Aeq, Beq.
0061 %
0062 %   Calling syntax options:
0063 %       [x, f, exitflag, output, lambda] = ...
0064 %           qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
0065 %
0066 %       x = qps_gurobi(H, c, A, l, u)
0067 %       x = qps_gurobi(H, c, A, l, u, xmin, xmax)
0068 %       x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0)
0069 %       x = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt)
0070 %       x = qps_gurobi(problem), where problem is a struct with fields:
0071 %                       H, c, A, l, u, xmin, xmax, x0, opt
0072 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0073 %       x = qps_gurobi(...)
0074 %       [x, f] = qps_gurobi(...)
0075 %       [x, f, exitflag] = qps_gurobi(...)
0076 %       [x, f, exitflag, output] = qps_gurobi(...)
0077 %       [x, f, exitflag, output, lambda] = qps_gurobi(...)
0078 %
0079 %
0080 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0081 %       H = [   1003.1  4.3     6.3     5.9;
0082 %               4.3     2.2     2.1     3.9;
0083 %               6.3     2.1     3.5     4.8;
0084 %               5.9     3.9     4.8     10  ];
0085 %       c = zeros(4,1);
0086 %       A = [   1       1       1       1;
0087 %               0.17    0.11    0.10    0.18    ];
0088 %       l = [1; 0.10];
0089 %       u = [1; Inf];
0090 %       xmin = zeros(4,1);
0091 %       x0 = [1; 0; 0; 1];
0092 %       opt = struct('verbose', 2);
0093 %       [x, f, s, out, lambda] = qps_gurobi(H, c, A, l, u, xmin, [], x0, opt);
0094 %
0095 %   See also QPS_MASTER, GUROBI.
0096 
0097 %   MP-Opt-Model
0098 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0099 %   by Ray Zimmerman, PSERC Cornell
0100 %
0101 %   This file is part of MP-Opt-Model.
0102 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0103 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0104 
0105 %%----- input argument handling  -----
0106 %% gather inputs
0107 if nargin == 1 && isstruct(H)       %% problem struct
0108     p = H;
0109     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0110     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0111     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0112     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0113     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0114     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0115     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0116     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0117     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0118 else                                %% individual args
0119     if nargin < 9
0120         opt = [];
0121         if nargin < 8
0122             x0 = [];
0123             if nargin < 7
0124                 xmax = [];
0125                 if nargin < 6
0126                     xmin = [];
0127                 end
0128             end
0129         end
0130     end
0131 end
0132 
0133 %% define nx, set default values for missing optional inputs
0134 if isempty(H) || ~any(any(H))
0135     if isempty(A) && isempty(xmin) && isempty(xmax)
0136         error('qps_gurobi: LP problem must include constraints or variable bounds');
0137     else
0138         if ~isempty(A)
0139             nx = size(A, 2);
0140         elseif ~isempty(xmin)
0141             nx = length(xmin);
0142         else    % if ~isempty(xmax)
0143             nx = length(xmax);
0144         end
0145     end
0146 else
0147     nx = size(H, 1);
0148 end
0149 if isempty(c)
0150     c = zeros(nx, 1);
0151 end
0152 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0153                                  (isempty(u) || all(u == Inf)))
0154     A = sparse(0,nx);           %% no limits => no linear constraints
0155 end
0156 nA = size(A, 1);                %% number of original linear constraints
0157 if isempty(u)                   %% By default, linear inequalities are ...
0158     u = Inf(nA, 1);             %% ... unbounded above and ...
0159 end
0160 if isempty(l)
0161     l = -Inf(nA, 1);            %% ... unbounded below.
0162 end
0163 if isempty(xmin)                %% By default, optimization variables are ...
0164     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0165 end
0166 if isempty(xmax)
0167     xmax = Inf(nx, 1);          %% ... unbounded above.
0168 end
0169 if isempty(x0)
0170     x0 = zeros(nx, 1);
0171 end
0172 
0173 %% default options
0174 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0175     verbose = opt.verbose;
0176 else
0177     verbose = 0;
0178 end
0179 
0180 %% set up options struct for Gurobi
0181 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt)
0182     g_opt = gurobi_options(opt.grb_opt);
0183 else
0184     g_opt = gurobi_options;
0185 end
0186 if verbose > 1
0187     g_opt.LogToConsole = 1;
0188     g_opt.OutputFlag = 1;
0189     if verbose > 2
0190         g_opt.DisplayInterval = 1;
0191     else
0192         g_opt.DisplayInterval = 100;
0193     end
0194 else
0195     g_opt.LogToConsole = 0;
0196     g_opt.OutputFlag = 0;
0197 end
0198 
0199 if ~issparse(A)
0200     A = sparse(A);
0201 end
0202 if issparse(c);
0203     c = full(c);
0204 end
0205 
0206 %% split up linear constraints
0207 ieq = find( abs(u-l) <= eps );          %% equality
0208 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0209 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0210 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0211 
0212 %% grab some dimensions
0213 nlt = length(ilt);      %% number of upper bounded linear inequalities
0214 ngt = length(igt);      %% number of lower bounded linear inequalities
0215 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0216 neq = length(ieq);      %% number of equalities
0217 niq = nlt+ngt+2*nbx;    %% number of inequalities
0218 
0219 %% set up model
0220 m.A     = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0221 m.rhs   = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)    ];
0222 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]);
0223 m.lb = xmin;
0224 m.ub = xmax;
0225 m.obj = c';
0226 
0227 %% call the solver
0228 if isempty(H) || ~any(any(H))
0229     lpqp = 'LP';
0230 else
0231     lpqp = 'QP';
0232     if ~issparse(H)
0233         H = sparse(H);
0234     end
0235     m.Q = 0.5 * H;
0236 end
0237 if verbose
0238     alg_names = {
0239         'automatic',
0240         'primal simplex',
0241         'dual simplex',
0242         'interior point',
0243         'concurrent',
0244         'deterministic concurrent'
0245     };
0246     vn = gurobiver;
0247     fprintf('Gurobi Version %s -- %s %s solver\n', ...
0248         vn, alg_names{g_opt.Method+2}, lpqp);
0249 end
0250 results = gurobi(m, g_opt);
0251 switch results.status
0252     case 'LOADED',          %% 1
0253         eflag = -1;
0254     case 'OPTIMAL',         %% 2, optimal solution found
0255         eflag = 1;
0256     case 'INFEASIBLE',      %% 3
0257         eflag = -3;
0258     case 'INF_OR_UNBD',     %% 4
0259         eflag = -4;
0260     case 'UNBOUNDED',       %% 5
0261         eflag = -5;
0262     case 'CUTOFF',          %% 6
0263         eflag = -6;
0264     case 'ITERATION_LIMIT', %% 7
0265         eflag = -7;
0266     case 'NODE_LIMIT',      %% 8
0267         eflag = -8;
0268     case 'TIME_LIMIT',      %% 9
0269         eflag = -9;
0270     case 'SOLUTION_LIMIT',  %% 10
0271         eflag = -10;
0272     case 'INTERRUPTED',     %% 11
0273         eflag = -11;
0274     case 'NUMERIC',         %% 12
0275         eflag = -12;
0276     case 'SUBOPTIMAL',      %% 13
0277         eflag = -13;
0278     case 'INPROGRESS',      %% 14
0279         eflag = -14;
0280     otherwise,
0281         eflag = 0;
0282 end
0283 output = results;
0284 
0285 %% check for empty results (in case optimization failed)
0286 if ~isfield(results, 'x') || isempty(results.x)
0287     x = NaN(nx, 1);
0288     lam.lower   = NaN(nx, 1);
0289     lam.upper   = NaN(nx, 1);
0290 else
0291     x = results.x;
0292     lam.lower   = zeros(nx, 1);
0293     lam.upper   = zeros(nx, 1);
0294 end
0295 if ~isfield(results, 'objval') || isempty(results.objval)
0296     f = NaN;
0297 else
0298     f = results.objval;
0299 end
0300 if ~isfield(results, 'pi') || isempty(results.pi)
0301     pi  = NaN(length(m.rhs), 1);
0302 else
0303     pi  = results.pi;
0304 end
0305 if ~isfield(results, 'rc') || isempty(results.rc)
0306     rc  = NaN(nx, 1);
0307 else
0308     rc  = results.rc;
0309 end
0310 
0311 kl = find(rc > 0);   %% lower bound binding
0312 ku = find(rc < 0);   %% upper bound binding
0313 lam.lower(kl)   =  rc(kl);
0314 lam.upper(ku)   = -rc(ku);
0315 lam.eqlin   = pi(1:neq);
0316 lam.ineqlin = pi(neq+(1:niq));
0317 mu_l        = zeros(nA, 1);
0318 mu_u        = zeros(nA, 1);
0319 
0320 %% repackage lambdas
0321 kl = find(lam.eqlin > 0);   %% lower bound binding
0322 ku = find(lam.eqlin < 0);   %% upper bound binding
0323 
0324 mu_l(ieq(kl)) = lam.eqlin(kl);
0325 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
0326 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0327 
0328 mu_u(ieq(ku)) = -lam.eqlin(ku);
0329 mu_u(ilt) = -lam.ineqlin(1:nlt);
0330 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));
0331 
0332 lambda = struct( ...
0333     'mu_l', mu_l, ...
0334     'mu_u', mu_u, ...
0335     'lower', lam.lower, ...
0336     'upper', lam.upper ...
0337 );

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005