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

miqps_glpk

PURPOSE ^

MIQPS_GLPK Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit.

SYNOPSIS ^

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

DESCRIPTION ^

MIQPS_GLPK  Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIQPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GLPK(PROBLEM)
   A wrapper function providing a standardized interface for using
   GLKP to solve the following LP (linear programming) problem:

       min 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 : dummy matrix (possibly sparse) of quadratic cost coefficients
           for QP problems, which GLPK does not handle
       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 (NOT USED)
       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) or
               'I' (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
           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
           glpk_opt - options struct for GLPK, 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 : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other
       OUTPUT : struct with errnum and status fields from GLPK output args
       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 GLPK. 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_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)

       x = miqps_glpk(H, c, A, l, u)
       x = miqps_glpk(H, c, A, l, u, xmin, xmax)
       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0)
       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype)
       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
       x = miqps_glpk(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_glpk(...)
       [x, f] = miqps_glpk(...)
       [x, f, exitflag] = miqps_glpk(...)
       [x, f, exitflag, output] = miqps_glpk(...)
       [x, f, exitflag, output, lambda] = miqps_glpk(...)


   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_glpk(H, c, A, l, u, xmin, [], x0, vtype, opt);

   See also MIQPS_MASTER, GLPK.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0002 %MIQPS_GLPK  Mixed Integer Linear Program Solver based on GLPK - GNU Linear Programming Kit.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIQPS_GLPK(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_GLPK(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   GLKP to solve the following LP (linear programming) problem:
0008 %
0009 %       min 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 : dummy matrix (possibly sparse) of quadratic cost coefficients
0019 %           for QP problems, which GLPK does not handle
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 (NOT USED)
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) or
0030 %               'I' (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 %           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 %           glpk_opt - options struct for GLPK, value in verbose
0046 %                   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 : exit flag, 1 - optimal, <= 0 - infeasible, unbounded or other
0055 %       OUTPUT : struct with errnum and status fields from GLPK output args
0056 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0057 %           multipliers on the constraints, with fields:
0058 %           mu_l - lower (left-hand) limit on linear constraints
0059 %           mu_u - upper (right-hand) limit on linear constraints
0060 %           lower - lower bound on optimization variables
0061 %           upper - upper bound on optimization variables
0062 %
0063 %   Note the calling syntax is almost identical to that of GLPK. The main
0064 %   difference is that the linear constraints are specified with A, L, U
0065 %   instead of A, B, Aeq, Beq.
0066 %
0067 %   Calling syntax options:
0068 %       [x, f, exitflag, output, lambda] = ...
0069 %           miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0070 %
0071 %       x = miqps_glpk(H, c, A, l, u)
0072 %       x = miqps_glpk(H, c, A, l, u, xmin, xmax)
0073 %       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0)
0074 %       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype)
0075 %       x = miqps_glpk(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0076 %       x = miqps_glpk(problem), where problem is a struct with fields:
0077 %                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
0078 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0079 %       x = miqps_glpk(...)
0080 %       [x, f] = miqps_glpk(...)
0081 %       [x, f, exitflag] = miqps_glpk(...)
0082 %       [x, f, exitflag, output] = miqps_glpk(...)
0083 %       [x, f, exitflag, output, lambda] = miqps_glpk(...)
0084 %
0085 %
0086 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0087 %       H = [   1003.1  4.3     6.3     5.9;
0088 %               4.3     2.2     2.1     3.9;
0089 %               6.3     2.1     3.5     4.8;
0090 %               5.9     3.9     4.8     10  ];
0091 %       c = zeros(4,1);
0092 %       A = [   1       1       1       1;
0093 %               0.17    0.11    0.10    0.18    ];
0094 %       l = [1; 0.10];
0095 %       u = [1; Inf];
0096 %       xmin = zeros(4,1);
0097 %       x0 = [1; 0; 0; 1];
0098 %       opt = struct('verbose', 2);
0099 %       [x, f, s, out, lambda] = miqps_glpk(H, c, A, l, u, xmin, [], x0, vtype, opt);
0100 %
0101 %   See also MIQPS_MASTER, GLPK.
0102 
0103 %   MP-Opt-Model
0104 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0105 %   by Ray Zimmerman, PSERC Cornell
0106 %
0107 %   This file is part of MP-Opt-Model.
0108 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0109 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0110 
0111 %% check for Optimization Toolbox
0112 % if ~have_feature('quadprog')
0113 %     error('miqps_glpk: requires the MEX interface to GLPK');
0114 % end
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_glpk: 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     error('miqps_glpk: GLPK handles only LP problems, not QP problems');
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 %% split up linear constraints
0197 ieq = find( abs(u-l) <= eps );          %% equality
0198 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0199 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0200 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0201 AA = [ A(ieq, :); A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0202 bb = [ u(ieq);    u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0203 
0204 %% grab some dimensions
0205 nlt = length(ilt);      %% number of upper bounded linear inequalities
0206 ngt = length(igt);      %% number of lower bounded linear inequalities
0207 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0208 neq = length(ieq);      %% number of equalities
0209 nie = nlt+ngt+2*nbx;    %% number of inequalities
0210 
0211 ctype = [repmat('S', neq, 1); repmat('U', nlt+ngt+2*nbx, 1)];
0212 
0213 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I'))
0214     mi = 0;
0215     vtype = repmat('C', nx, 1);
0216 else
0217     mi = 1;
0218     %% expand vtype to nx elements if necessary
0219     if length(vtype) == 1 && nx > 1
0220         vtype = char(vtype * ones(nx, 1));
0221     elseif size(vtype, 2) > 1   %% make sure it's a col vector
0222         vtype = vtype';
0223     end
0224 end
0225 %% convert 'B' variables to 'I' and clip bounds to [0, 1]
0226 k = find(vtype == 'B');
0227 if ~isempty(k)
0228     kk = find(xmax(k) > 1);
0229     xmax(k(kk)) = 1;
0230     kk = find(xmin(k) < 0);
0231     xmin(k(kk)) = 0;
0232     vtype(k) = 'I';
0233 end
0234 
0235 %% set options struct for GLPK
0236 if ~isempty(opt) && isfield(opt, 'glpk_opt') && ~isempty(opt.glpk_opt)
0237     glpk_opt = glpk_options(opt.glpk_opt);
0238 else
0239     glpk_opt = glpk_options;
0240 end
0241 glpk_opt.msglev = verbose;
0242 
0243 %% call the solver
0244 [x, f, errnum, extra] = ...
0245     glpk(c, AA, bb, xmin, xmax, ctype, vtype, 1, glpk_opt);
0246 
0247 %% set exit flag
0248 if isfield(extra, 'status')             %% status found in extra.status
0249     output.errnum = errnum;
0250     output.status = extra.status;
0251     eflag = -errnum;
0252     if eflag == 0 && extra.status == 5
0253         eflag = 1;
0254     end
0255 else                                    %% status found in errnum
0256     output.errnum = [];
0257     output.status = errnum;
0258     if have_feature('octave')
0259         if errnum == 180 || errnum == 151 || errnum == 171
0260             eflag = 1;
0261         else
0262             eflag = 0;
0263         end
0264     else
0265         if errnum == 5
0266             eflag = 1;
0267         else
0268             eflag = 0;
0269         end
0270     end
0271 end
0272 
0273 %% repackage lambdas
0274 if isempty(extra) || ~isfield(extra, 'lambda') || isempty(extra.lambda)
0275     lambda = struct( ...
0276         'mu_l', zeros(nA, 1), ...
0277         'mu_u', zeros(nA, 1), ...
0278         'lower', zeros(nx, 1), ...
0279         'upper', zeros(nx, 1) ...
0280     );
0281 else
0282     lam.eqlin = extra.lambda(1:neq);
0283     lam.ineqlin = extra.lambda(neq+(1:nie));
0284     lam.lower = extra.redcosts;
0285     lam.upper = -extra.redcosts;
0286     lam.lower(lam.lower < 0) = 0;
0287     lam.upper(lam.upper < 0) = 0;
0288     kl = find(lam.eqlin > 0);   %% lower bound binding
0289     ku = find(lam.eqlin < 0);   %% upper bound binding
0290 
0291     mu_l = zeros(nA, 1);
0292     mu_l(ieq(kl)) = lam.eqlin(kl);
0293     mu_l(igt) = -lam.ineqlin(nlt+(1:ngt));
0294     mu_l(ibx) = -lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0295 
0296     mu_u = zeros(nA, 1);
0297     mu_u(ieq(ku)) = -lam.eqlin(ku);
0298     mu_u(ilt) = -lam.ineqlin(1:nlt);
0299     mu_u(ibx) = -lam.ineqlin(nlt+ngt+(1:nbx));
0300 
0301     lambda = struct( ...
0302         'mu_l', mu_l, ...
0303         'mu_u', mu_u, ...
0304         'lower', lam.lower(1:nx), ...
0305         'upper', lam.upper(1:nx) ...
0306     );
0307 end
0308 
0309 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
0310     if verbose
0311         fprintf('--- Integer stage complete, starting price computation stage ---\n');
0312     end
0313     if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
0314         tol = opt.price_stage_warn_tol;
0315     else
0316         tol = 1e-7;
0317     end
0318     k = find(vtype == 'I' | vtype == 'B');
0319     x(k) = round(x(k));
0320     xmin(k) = x(k);
0321     xmax(k) = x(k);
0322     x0 = x;
0323     opt.glpk_opt.lpsolver = 1;      %% simplex
0324     opt.glpk_opt.dual = 0;          %% primal simplex
0325     if have_feature('octave') && have_feature('octave', 'vnum') >= 3.007
0326         opt.glpk_opt.dual = 1;      %% primal simplex
0327     end
0328     
0329     [x_, f_, eflag_, output_, lambda] = qps_glpk(H, c, A, l, u, xmin, xmax, x0, opt);
0330     if eflag ~= eflag_
0331         error('miqps_glpk: EXITFLAG from price computation stage = %d', eflag_);
0332     end
0333     if abs(f - f_)/max(abs(f), 1) > tol
0334         warning('miqps_glpk: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
0335     end
0336     xn = x;
0337     xn(abs(xn)<1) = 1;
0338     [mx, k] = max(abs(x - x_) ./ xn);
0339     if mx > tol
0340         warning('miqps_glpk: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
0341     end
0342     output.price_stage = output_;
0343 end

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