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

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