Home > matpower7.1 > mp-opt-model > lib > miqps_gurobi.m

miqps_gurobi

PURPOSE ^

MIQPS_GUROBI Mixed Integer Quadratic Program Solver based on GUROBI.

SYNOPSIS ^

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

DESCRIPTION ^

MIQPS_GUROBI  Mixed Integer Quadratic Program Solver based on GUROBI.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIQPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_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
       VTYPE : character string of length NX (number of elements in X),
               or 1 (value applies to all variables in x),
               allowed values are 'C' (continuous), 'B' (binary),
               'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
       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
           skip_prices (0) - flag that specifies whether or not to
               skip the price computation stage, in which the problem
               is re-solved for only the continuous variables, with all
               others being constrained to their solved values
           price_stage_warn_tol (1e-7) - tolerance on the objective fcn
               value and primal variable relative match required to avoid
               mis-match warning message
           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, vtype, 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] = ...
           miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)

       x = miqps_gurobi(H, c, A, l, u)
       x = miqps_gurobi(H, c, A, l, u, xmin, xmax)
       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0)
       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype)
       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
       x = miqps_gurobi(problem), where problem is a struct with fields:
                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
                       all fields except 'c', 'A' and 'l' or 'u' are optional
       x = miqps_gurobi(...)
       [x, f] = miqps_gurobi(...)
       [x, f, exitflag] = miqps_gurobi(...)
       [x, f, exitflag, output] = miqps_gurobi(...)
       [x, f, exitflag, output, lambda] = miqps_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] = miqps_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt);

   See also MIQPS_MASTER, GUROBI.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0002 %MIQPS_GUROBI  Mixed Integer Quadratic Program Solver based on GUROBI.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIQPS_GUROBI(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_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 %       VTYPE : character string of length NX (number of elements in X),
0028 %               or 1 (value applies to all variables in x),
0029 %               allowed values are 'C' (continuous), 'B' (binary),
0030 %               'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
0031 %       OPT : optional options structure with the following fields,
0032 %           all of which are also optional (default values shown in
0033 %           parentheses)
0034 %           verbose (0) - controls level of progress output displayed
0035 %               0 = no progress output
0036 %               1 = some progress output
0037 %               2 = verbose progress output
0038 %               3 = even more verbose progress output
0039 %           skip_prices (0) - flag that specifies whether or not to
0040 %               skip the price computation stage, in which the problem
0041 %               is re-solved for only the continuous variables, with all
0042 %               others being constrained to their solved values
0043 %           price_stage_warn_tol (1e-7) - tolerance on the objective fcn
0044 %               value and primal variable relative match required to avoid
0045 %               mis-match warning message
0046 %           grb_opt - options struct for GUROBI, value in verbose
0047 %                   overrides these options
0048 %       PROBLEM : The inputs can alternatively be supplied in a single
0049 %           PROBLEM struct with fields corresponding to the input arguments
0050 %           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
0051 %
0052 %   Outputs:
0053 %       X : solution vector
0054 %       F : final objective function value
0055 %       EXITFLAG : GUROBI exit flag
0056 %           1 = converged
0057 %           0 or negative values = negative of GUROBI exit flag
0058 %           (see GUROBI documentation for details)
0059 %       OUTPUT : GUROBI output struct
0060 %           (see GUROBI documentation for details)
0061 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0062 %           multipliers on the constraints, with fields:
0063 %           mu_l - lower (left-hand) limit on linear constraints
0064 %           mu_u - upper (right-hand) limit on linear constraints
0065 %           lower - lower bound on optimization variables
0066 %           upper - upper bound on optimization variables
0067 %
0068 %   Note the calling syntax is almost identical to that of QUADPROG
0069 %   from MathWorks' Optimization Toolbox. The main difference is that
0070 %   the linear constraints are specified with A, L, U instead of
0071 %   A, B, Aeq, Beq.
0072 %
0073 %   Calling syntax options:
0074 %       [x, f, exitflag, output, lambda] = ...
0075 %           miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0076 %
0077 %       x = miqps_gurobi(H, c, A, l, u)
0078 %       x = miqps_gurobi(H, c, A, l, u, xmin, xmax)
0079 %       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0)
0080 %       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype)
0081 %       x = miqps_gurobi(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0082 %       x = miqps_gurobi(problem), where problem is a struct with fields:
0083 %                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
0084 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0085 %       x = miqps_gurobi(...)
0086 %       [x, f] = miqps_gurobi(...)
0087 %       [x, f, exitflag] = miqps_gurobi(...)
0088 %       [x, f, exitflag, output] = miqps_gurobi(...)
0089 %       [x, f, exitflag, output, lambda] = miqps_gurobi(...)
0090 %
0091 %
0092 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0093 %       H = [   1003.1  4.3     6.3     5.9;
0094 %               4.3     2.2     2.1     3.9;
0095 %               6.3     2.1     3.5     4.8;
0096 %               5.9     3.9     4.8     10  ];
0097 %       c = zeros(4,1);
0098 %       A = [   1       1       1       1;
0099 %               0.17    0.11    0.10    0.18    ];
0100 %       l = [1; 0.10];
0101 %       u = [1; Inf];
0102 %       xmin = zeros(4,1);
0103 %       x0 = [1; 0; 0; 1];
0104 %       opt = struct('verbose', 2);
0105 %       [x, f, s, out, lambda] = miqps_gurobi(H, c, A, l, u, xmin, [], x0, vtype, opt);
0106 %
0107 %   See also MIQPS_MASTER, GUROBI.
0108 
0109 %   MP-Opt-Model
0110 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0111 %   by Ray Zimmerman, PSERC Cornell
0112 %
0113 %   This file is part of MP-Opt-Model.
0114 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0115 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0116 
0117 %%----- input argument handling  -----
0118 %% gather inputs
0119 if nargin == 1 && isstruct(H)       %% problem struct
0120     p = H;
0121     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0122     if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
0123     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0124     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0125     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0126     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0127     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0128     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0129     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0130     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0131 else                                %% individual args
0132     if nargin < 10
0133         opt = [];
0134         if nargin < 9
0135             vtype = [];
0136             if nargin < 8
0137                 x0 = [];
0138                 if nargin < 7
0139                     xmax = [];
0140                     if nargin < 6
0141                         xmin = [];
0142                     end
0143                 end
0144             end
0145         end
0146     end
0147 end
0148 
0149 %% define nx, set default values for missing optional inputs
0150 if isempty(H) || ~any(any(H))
0151     if isempty(A) && isempty(xmin) && isempty(xmax)
0152         error('miqps_gurobi: LP problem must include constraints or variable bounds');
0153     else
0154         if ~isempty(A)
0155             nx = size(A, 2);
0156         elseif ~isempty(xmin)
0157             nx = length(xmin);
0158         else    % if ~isempty(xmax)
0159             nx = length(xmax);
0160         end
0161     end
0162 else
0163     nx = size(H, 1);
0164 end
0165 if isempty(c)
0166     c = zeros(nx, 1);
0167 end
0168 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0169                                  (isempty(u) || all(u == Inf)))
0170     A = sparse(0,nx);           %% no limits => no linear constraints
0171 end
0172 nA = size(A, 1);                %% number of original linear constraints
0173 if isempty(u)                   %% By default, linear inequalities are ...
0174     u = Inf(nA, 1);             %% ... unbounded above and ...
0175 end
0176 if isempty(l)
0177     l = -Inf(nA, 1);            %% ... unbounded below.
0178 end
0179 if isempty(xmin)                %% By default, optimization variables are ...
0180     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0181 end
0182 if isempty(xmax)
0183     xmax = Inf(nx, 1);          %% ... unbounded above.
0184 end
0185 if isempty(x0)
0186     x0 = zeros(nx, 1);
0187 end
0188 
0189 %% default options
0190 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0191     verbose = opt.verbose;
0192 else
0193     verbose = 0;
0194 end
0195 
0196 %% set up options struct for Gurobi
0197 if ~isempty(opt) && isfield(opt, 'grb_opt') && ~isempty(opt.grb_opt)
0198     g_opt = gurobi_options(opt.grb_opt);
0199 else
0200     g_opt = gurobi_options;
0201 end
0202 if verbose > 1
0203     g_opt.LogToConsole = 1;
0204     g_opt.OutputFlag = 1;
0205     if verbose > 2
0206         g_opt.DisplayInterval = 1;
0207     else
0208         g_opt.DisplayInterval = 100;
0209     end
0210 else
0211     g_opt.LogToConsole = 0;
0212     g_opt.OutputFlag = 0;
0213 end
0214 
0215 if ~issparse(A)
0216     A = sparse(A);
0217 end
0218 if issparse(c);
0219     c = full(c);
0220 end
0221 
0222 %% split up linear constraints
0223 ieq = find( abs(u-l) <= eps );          %% equality
0224 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0225 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0226 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0227 
0228 %% grab some dimensions
0229 nlt = length(ilt);      %% number of upper bounded linear inequalities
0230 ngt = length(igt);      %% number of lower bounded linear inequalities
0231 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0232 neq = length(ieq);      %% number of equalities
0233 niq = nlt+ngt+2*nbx;    %% number of inequalities
0234 
0235 %% set up model
0236 m.A     = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0237 m.rhs   = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)    ];
0238 m.sense = char([ double('=')*ones(1,neq) double('<')*ones(1,niq) ]);
0239 m.lb = xmin;
0240 m.ub = xmax;
0241 m.obj = c';
0242 if ~isempty(vtype)
0243     m.vtype = vtype;
0244 end
0245 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ...
0246         vtype == 'S' | vtype == 'N'))
0247     mi = 0;
0248 else
0249     mi = 1;
0250 end
0251 
0252 %% call the solver
0253 if isempty(H) || ~any(any(H))
0254     lpqp = 'LP';
0255 else
0256     lpqp = 'QP';
0257     if ~issparse(H)
0258         H = sparse(H);
0259     end
0260     m.Q = 0.5 * H;
0261 end
0262 if mi
0263     lpqp = ['MI' lpqp];
0264 end
0265 if verbose
0266     alg_names = {
0267         'automatic',
0268         'primal simplex',
0269         'dual simplex',
0270         'interior point',
0271         'concurrent',
0272         'deterministic concurrent'
0273     };
0274     vn = gurobiver;
0275     fprintf('Gurobi Version %s -- %s %s solver\n', ...
0276         vn, alg_names{g_opt.Method+2}, lpqp);
0277 end
0278 results = gurobi(m, g_opt);
0279 switch results.status
0280     case 'LOADED',          %% 1
0281         eflag = -1;
0282     case 'OPTIMAL',         %% 2, optimal solution found
0283         eflag = 1;
0284     case 'INFEASIBLE',      %% 3
0285         eflag = -3;
0286     case 'INF_OR_UNBD',     %% 4
0287         eflag = -4;
0288     case 'UNBOUNDED',       %% 5
0289         eflag = -5;
0290     case 'CUTOFF',          %% 6
0291         eflag = -6;
0292     case 'ITERATION_LIMIT', %% 7
0293         eflag = -7;
0294     case 'NODE_LIMIT',      %% 8
0295         eflag = -8;
0296     case 'TIME_LIMIT',      %% 9
0297         eflag = -9;
0298     case 'SOLUTION_LIMIT',  %% 10
0299         eflag = -10;
0300     case 'INTERRUPTED',     %% 11
0301         eflag = -11;
0302     case 'NUMERIC',         %% 12
0303         eflag = -12;
0304     case 'SUBOPTIMAL',      %% 13
0305         eflag = -13;
0306     case 'INPROGRESS',      %% 14
0307         eflag = -14;
0308     otherwise,
0309         eflag = 0;
0310 end
0311 output = results;
0312 
0313 %% check for empty results (in case optimization failed)
0314 if ~isfield(results, 'x') || isempty(results.x)
0315     x = NaN(nx, 1);
0316     lam.lower   = NaN(nx, 1);
0317     lam.upper   = NaN(nx, 1);
0318 else
0319     x = results.x;
0320     lam.lower   = zeros(nx, 1);
0321     lam.upper   = zeros(nx, 1);
0322 end
0323 if ~isfield(results, 'objval') || isempty(results.objval)
0324     f = NaN;
0325 else
0326     f = results.objval;
0327 end
0328 if ~isfield(results, 'pi') || isempty(results.pi)
0329     pi  = NaN(length(m.rhs), 1);
0330 else
0331     pi  = results.pi;
0332 end
0333 if ~isfield(results, 'rc') || isempty(results.rc)
0334     rc  = NaN(nx, 1);
0335 else
0336     rc  = results.rc;
0337 end
0338 
0339 kl = find(rc > 0);   %% lower bound binding
0340 ku = find(rc < 0);   %% upper bound binding
0341 lam.lower(kl)   =  rc(kl);
0342 lam.upper(ku)   = -rc(ku);
0343 lam.eqlin   = pi(1:neq);
0344 lam.ineqlin = pi(neq+(1:niq));
0345 mu_l        = zeros(nA, 1);
0346 mu_u        = zeros(nA, 1);
0347 
0348 %% repackage lambdas
0349 kl = find(lam.eqlin > 0);   %% lower bound binding
0350 ku = find(lam.eqlin < 0);   %% upper bound binding
0351 
0352 mu_l(ieq(kl)) = lam.eqlin(kl);
0353 mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
0354 mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0355 
0356 mu_u(ieq(ku)) = -lam.eqlin(ku);
0357 mu_u(ilt) = -lam.ineqlin(1:nlt);
0358 mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));
0359 
0360 lambda = struct( ...
0361     'mu_l', mu_l, ...
0362     'mu_u', mu_u, ...
0363     'lower', lam.lower, ...
0364     'upper', lam.upper ...
0365 );
0366 
0367 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
0368     if verbose
0369         fprintf('--- Integer stage complete, starting price computation stage ---\n');
0370     end
0371     if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
0372         tol = opt.price_stage_warn_tol;
0373     else
0374         tol = 1e-7;
0375     end
0376     if length(vtype) == 1
0377         if vtype == 'I' || vtype == 'B' || vtype == 'N'
0378             k = (1:nx);
0379         elseif vtype == 'S'
0380             k = find(x == 0);
0381         end
0382     else    %% length(vtype) == nx
0383         k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ...
0384                 (vtype == 'S' & x' == 0));
0385     end
0386     
0387     x(k) = round(x(k));
0388     xmin(k) = x(k);
0389     xmax(k) = x(k);
0390     x0 = x;
0391 %     opt.grb_opt.Method = 0;     %% primal simplex
0392     
0393     [x_, f_, eflag_, output_, lambda] = qps_gurobi(H, c, A, l, u, xmin, xmax, x0, opt);
0394     if eflag ~= eflag_
0395         error('miqps_gurobi: EXITFLAG from price computation stage = %d', eflag_);
0396     end
0397     if abs(f - f_)/max(abs(f), 1) > tol
0398         warning('miqps_gurobi: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
0399     end
0400     xn = x;
0401     xn(abs(xn)<1) = 1;
0402     [mx, k] = max(abs(x - x_) ./ xn);
0403     if mx > tol
0404         warning('miqps_gurobi: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
0405     end
0406     output.price_stage = output_;
0407 end

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