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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005