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

miqps_ot

PURPOSE ^

MIQPS_OT Mixed Integer Linear Program Solver based on INTLINPROG.

SYNOPSIS ^

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

DESCRIPTION ^

MIQPS_OT  Mixed Integer Linear Program Solver based on INTLINPROG.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIQPS_OT(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_OT(PROBLEM)
   A wrapper function providing a standardized interface for using
   QUADPROG or LINPROG from the Optimization Toolbox 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) 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
           intlinprog_opt - options struct for INTLINPROG, value in verbose
                   overrides these options
           linprog_opt - options struct for LINPROG, 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 : QUADPROG/LINPROG exit flag
           (see QUADPROG and LINPROG documentation for details)
       OUTPUT : QUADPROG/LINPROG output struct
           (see QUADPROG and LINPROG 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_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)

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


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

   See also MIQPS_MASTER, INTLINPROG, QUADPROG, LINPROG.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0002 %MIQPS_OT  Mixed Integer Linear Program Solver based on INTLINPROG.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIQPS_OT(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = MIQPS_OT(PROBLEM)
0006 %   A wrapper function providing a standardized interface for using
0007 %   QUADPROG or LINPROG from the Optimization Toolbox to solve the
0008 %   following QP (quadratic programming) 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) 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 %           intlinprog_opt - options struct for INTLINPROG, value in verbose
0046 %                   overrides these options
0047 %           linprog_opt - options struct for LINPROG, value in verbose
0048 %                   overrides these options
0049 %       PROBLEM : The inputs can alternatively be supplied in a single
0050 %           PROBLEM struct with fields corresponding to the input arguments
0051 %           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
0052 %
0053 %   Outputs:
0054 %       X : solution vector
0055 %       F : final objective function value
0056 %       EXITFLAG : QUADPROG/LINPROG exit flag
0057 %           (see QUADPROG and LINPROG documentation for details)
0058 %       OUTPUT : QUADPROG/LINPROG output struct
0059 %           (see QUADPROG and LINPROG 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_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0075 %
0076 %       x = miqps_ot(H, c, A, l, u)
0077 %       x = miqps_ot(H, c, A, l, u, xmin, xmax)
0078 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0)
0079 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype)
0080 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0081 %       x = miqps_ot(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_ot(...)
0085 %       [x, f] = miqps_ot(...)
0086 %       [x, f, exitflag] = miqps_ot(...)
0087 %       [x, f, exitflag, output] = miqps_ot(...)
0088 %       [x, f, exitflag, output, lambda] = miqps_ot(...)
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_ot(H, c, A, l, u, xmin, [], x0, vtype, opt);
0105 %
0106 %   See also MIQPS_MASTER, INTLINPROG, QUADPROG, LINPROG.
0107 
0108 %   MP-Opt-Model
0109 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0110 %   by Ray Zimmerman, PSERC Cornell
0111 %
0112 %   This file is part of MP-Opt-Model.
0113 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0114 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0115 
0116 %% check for Optimization Toolbox
0117 % if ~have_feature('quadprog')
0118 %     error('miqps_ot: requires the Optimization Toolbox');
0119 % end
0120 
0121 %%----- input argument handling  -----
0122 %% gather inputs
0123 if nargin == 1 && isstruct(H)       %% problem struct
0124     p = H;
0125     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0126     if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
0127     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0128     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0129     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0130     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0131     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0132     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0133     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0134     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0135 else                                %% individual args
0136     if nargin < 10
0137         opt = [];
0138         if nargin < 9
0139             vtype = [];
0140             if nargin < 8
0141                 x0 = [];
0142                 if nargin < 7
0143                     xmax = [];
0144                     if nargin < 6
0145                         xmin = [];
0146                     end
0147                 end
0148             end
0149         end
0150     end
0151 end
0152 
0153 %% define nx, set default values for missing optional inputs
0154 if isempty(H) || ~any(any(H))
0155     if isempty(A) && isempty(xmin) && isempty(xmax)
0156         error('miqps_ot: LP problem must include constraints or variable bounds');
0157     else
0158         if ~isempty(A)
0159             nx = size(A, 2);
0160         elseif ~isempty(xmin)
0161             nx = length(xmin);
0162         else    % if ~isempty(xmax)
0163             nx = length(xmax);
0164         end
0165     end
0166 else
0167     nx = size(H, 1);
0168 end
0169 if isempty(c)
0170     c = zeros(nx, 1);
0171 end
0172 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0173                                  (isempty(u) || all(u == Inf)))
0174     A = sparse(0,nx);           %% no limits => no linear constraints
0175 end
0176 nA = size(A, 1);                %% number of original linear constraints
0177 if isempty(u)                   %% By default, linear inequalities are ...
0178     u = Inf(nA, 1);             %% ... unbounded above and ...
0179 end
0180 if isempty(l)
0181     l = -Inf(nA, 1);            %% ... unbounded below.
0182 end
0183 if isempty(xmin)                %% By default, optimization variables are ...
0184     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0185 end
0186 if isempty(xmax)
0187     xmax = Inf(nx, 1);          %% ... unbounded above.
0188 end
0189 if isempty(x0)
0190     x0 = zeros(nx, 1);
0191 end
0192 if isempty(H) || ~any(any(H))
0193     isLP = 1;   %% it's an LP
0194 else
0195     isLP = 0;   %% nope, it's a QP
0196 end
0197 
0198 %% default options
0199 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0200     verbose = opt.verbose;
0201 else
0202     verbose = 0;
0203 end
0204 
0205 %% split up linear constraints
0206 ieq = find( abs(u-l) <= eps );          %% equality
0207 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0208 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0209 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0210 Ae = A(ieq, :);
0211 be = u(ieq);
0212 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0213 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0214 
0215 %% grab some dimensions
0216 nlt = length(ilt);      %% number of upper bounded linear inequalities
0217 ngt = length(igt);      %% number of lower bounded linear inequalities
0218 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0219 
0220 %% mixed integer?
0221 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I'))
0222     mi = 0;
0223     vtype = repmat('C', 1, nx);
0224 else
0225     mi = 1;
0226     %% expand vtype to nx elements if necessary
0227     if length(vtype) == 1 && nx > 1
0228         vtype = char(vtype * ones(1, nx));
0229     end
0230 end
0231 %% check bounds for 'B' variables to make sure they are between 0 and 1
0232 k = find(vtype == 'B');
0233 if ~isempty(k)
0234     kk = find(xmax(k) > 1);
0235     xmax(k(kk)) = 1;
0236     kk = find(xmin(k) < 0);
0237     xmin(k(kk)) = 0;
0238 end
0239 intcon = find(vtype == 'B' | vtype == 'I');
0240 
0241 %% set up options
0242 if verbose > 1
0243     vrb = 'iter';       %% seems to be same as 'final'
0244 elseif verbose == 1
0245     vrb = 'final';
0246 else
0247     vrb = 'off';
0248 end
0249 
0250 if isLP
0251     if mi
0252         ot_opt = optimoptions('intlinprog');
0253         if ~isempty(opt) && isfield(opt, 'intlinprog_opt') && ~isempty(opt.intlinprog_opt)
0254             ot_opt = nested_struct_copy(ot_opt, opt.intlinprog_opt);
0255         end
0256     else
0257         ot_opt = optimoptions('linprog');
0258         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0259             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0260         end
0261     end
0262 else
0263     if mi
0264         error('miqps_ot: MIQP problems not supported.');
0265     end
0266     ot_opt = optimoptions('quadprog');
0267     if have_feature('quadprog_ls')
0268         ot_opt.Algorithm = 'interior-point-convex';
0269     else
0270         ot_opt.LargeScale = 'off';
0271     end
0272     if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0273         ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0274     end
0275 end
0276 ot_opt = optimoptions(ot_opt, 'Display', vrb);
0277 
0278 %% call the solver
0279 if isLP
0280     if mi
0281         [x, f, eflag, output] = ...
0282             intlinprog(c, intcon, Ai, bi, Ae, be, xmin, xmax, ot_opt);
0283         lam = [];
0284     else
0285         [x, f, eflag, output, lam] = ...
0286             linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0287     end
0288 else
0289     [x, f, eflag, output, lam] = ...
0290         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0291 end
0292 
0293 %% repackage lambdas
0294 if isempty(x)
0295     x = NaN(nx, 1);
0296 end
0297 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ...
0298                     isempty(lam.lower) && isempty(lam.upper))
0299     lambda = struct( ...
0300         'mu_l', NaN(nA, 1), ...
0301         'mu_u', NaN(nA, 1), ...
0302         'lower', NaN(nx, 1), ...
0303         'upper', NaN(nx, 1) ...
0304     );
0305 else
0306     kl = find(lam.eqlin < 0);   %% lower bound binding
0307     ku = find(lam.eqlin > 0);   %% upper bound binding
0308 
0309     mu_l = zeros(nA, 1);
0310     mu_l(ieq(kl)) = -lam.eqlin(kl);
0311     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0312     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0313 
0314     mu_u = zeros(nA, 1);
0315     mu_u(ieq(ku)) = lam.eqlin(ku);
0316     mu_u(ilt) = lam.ineqlin(1:nlt);
0317     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0318 
0319     lambda = struct( ...
0320         'mu_l', mu_l, ...
0321         'mu_u', mu_u, ...
0322         'lower', lam.lower(1:nx), ...
0323         'upper', lam.upper(1:nx) ...
0324     );
0325 end
0326 
0327 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
0328     if verbose
0329         fprintf('--- Integer stage complete, starting price computation stage ---\n');
0330     end
0331     if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
0332         tol = opt.price_stage_warn_tol;
0333     else
0334         tol = 1e-7;
0335     end
0336     k = intcon;
0337     x(k) = round(x(k));
0338     xmin(k) = x(k);
0339     xmax(k) = x(k);
0340     x0 = x;
0341 %     opt.linprog_opt.Algorithm = 'dual-simplex';     %% dual-simplex
0342     
0343     [x_, f_, eflag_, output_, lambda] = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt);
0344 %     output
0345 %     output.message
0346     if eflag ~= eflag_
0347         error('miqps_ot: EXITFLAG from price computation stage = %d', eflag_);
0348     end
0349     if abs(f - f_)/max(abs(f), 1) > tol
0350         warning('miqps_ot: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
0351     end
0352     xn = x;
0353     xn(abs(xn)<1) = 1;
0354     [mx, k] = max(abs(x - x_) ./ xn);
0355     if mx > tol
0356         warning('miqps_ot: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
0357     end
0358     output.price_stage = output_;
0359 %     output_
0360 %     output_.message
0361 end

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