Home > matpower7.0 > 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, OPT)
   A wrapper function providing a MATPOWER 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 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, OPT)
0005 %   A wrapper function providing a MATPOWER standardized interface for using
0006 %   QUADPROG or LINPROG from the Optimization Toolbox to solve the
0007 %   following QP (quadratic programming) 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) or
0029 %               'I' (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 %           skip_prices (0) - flag that specifies whether or not to
0038 %               skip the price computation stage, in which the problem
0039 %               is re-solved for only the continuous variables, with all
0040 %               others being constrained to their solved values
0041 %           price_stage_warn_tol (1e-7) - tolerance on the objective fcn
0042 %               value and primal variable relative match required to avoid
0043 %               mis-match warning message
0044 %           intlinprog_opt - options struct for INTLINPROG, value in
0045 %               verbose overrides these options
0046 %           linprog_opt - options struct for LINPROG, value in
0047 %               verbose 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 : QUADPROG/LINPROG exit flag
0056 %           (see QUADPROG and LINPROG documentation for details)
0057 %       OUTPUT : QUADPROG/LINPROG output struct
0058 %           (see QUADPROG and LINPROG documentation for details)
0059 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0060 %           multipliers on the constraints, with fields:
0061 %           mu_l - lower (left-hand) limit on linear constraints
0062 %           mu_u - upper (right-hand) limit on linear constraints
0063 %           lower - lower bound on optimization variables
0064 %           upper - upper bound on optimization variables
0065 %
0066 %   Note the calling syntax is almost identical to that of QUADPROG
0067 %   from MathWorks' Optimization Toolbox. The main difference is that
0068 %   the linear constraints are specified with A, L, U instead of
0069 %   A, B, Aeq, Beq.
0070 %
0071 %   Calling syntax options:
0072 %       [x, f, exitflag, output, lambda] = ...
0073 %           miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0074 %
0075 %       x = miqps_ot(H, c, A, l, u)
0076 %       x = miqps_ot(H, c, A, l, u, xmin, xmax)
0077 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0)
0078 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype)
0079 %       x = miqps_ot(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0080 %       x = miqps_ot(problem), where problem is a struct with fields:
0081 %                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
0082 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0083 %       x = miqps_ot(...)
0084 %       [x, f] = miqps_ot(...)
0085 %       [x, f, exitflag] = miqps_ot(...)
0086 %       [x, f, exitflag, output] = miqps_ot(...)
0087 %       [x, f, exitflag, output, lambda] = miqps_ot(...)
0088 %
0089 %
0090 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0091 %       H = [   1003.1  4.3     6.3     5.9;
0092 %               4.3     2.2     2.1     3.9;
0093 %               6.3     2.1     3.5     4.8;
0094 %               5.9     3.9     4.8     10  ];
0095 %       c = zeros(4,1);
0096 %       A = [   1       1       1       1;
0097 %               0.17    0.11    0.10    0.18    ];
0098 %       l = [1; 0.10];
0099 %       u = [1; Inf];
0100 %       xmin = zeros(4,1);
0101 %       x0 = [1; 0; 0; 1];
0102 %       opt = struct('verbose', 2);
0103 %       [x, f, s, out, lambda] = miqps_ot(H, c, A, l, u, xmin, [], x0, vtype, opt);
0104 %
0105 %   See also INTLINPROG, QUADPROG, LINPROG.
0106 
0107 %   MATPOWER
0108 %   Copyright (c) 2010-2016, Power Systems Engineering Research Center (PSERC)
0109 %   by Ray Zimmerman, PSERC Cornell
0110 %
0111 %   This file is part of MATPOWER.
0112 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0113 %   See https://matpower.org for more info.
0114 
0115 %% check for Optimization Toolbox
0116 % if ~have_fcn('quadprog')
0117 %     error('miqps_ot: requires the Optimization Toolbox');
0118 % end
0119 
0120 %%----- input argument handling  -----
0121 %% gather inputs
0122 if nargin == 1 && isstruct(H)       %% problem struct
0123     p = H;
0124     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0125     if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
0126     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0127     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0128     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0129     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0130     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0131     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0132     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0133     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0134 else                                %% individual args
0135     if nargin < 10
0136         opt = [];
0137         if nargin < 9
0138             vtype = [];
0139             if nargin < 8
0140                 x0 = [];
0141                 if nargin < 7
0142                     xmax = [];
0143                     if nargin < 6
0144                         xmin = [];
0145                     end
0146                 end
0147             end
0148         end
0149     end
0150 end
0151 
0152 %% define nx, set default values for missing optional inputs
0153 if isempty(H) || ~any(any(H))
0154     if isempty(A) && isempty(xmin) && isempty(xmax)
0155         error('miqps_ot: LP problem must include constraints or variable bounds');
0156     else
0157         if ~isempty(A)
0158             nx = size(A, 2);
0159         elseif ~isempty(xmin)
0160             nx = length(xmin);
0161         else    % if ~isempty(xmax)
0162             nx = length(xmax);
0163         end
0164     end
0165 else
0166     nx = size(H, 1);
0167 end
0168 if isempty(c)
0169     c = zeros(nx, 1);
0170 end
0171 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0172                                  (isempty(u) || all(u == Inf)))
0173     A = sparse(0,nx);           %% no limits => no linear constraints
0174 end
0175 nA = size(A, 1);                %% number of original linear constraints
0176 if isempty(u)                   %% By default, linear inequalities are ...
0177     u = Inf(nA, 1);             %% ... unbounded above and ...
0178 end
0179 if isempty(l)
0180     l = -Inf(nA, 1);            %% ... unbounded below.
0181 end
0182 if isempty(xmin)                %% By default, optimization variables are ...
0183     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0184 end
0185 if isempty(xmax)
0186     xmax = Inf(nx, 1);          %% ... unbounded above.
0187 end
0188 if isempty(x0)
0189     x0 = zeros(nx, 1);
0190 end
0191 if isempty(H) || ~any(any(H))
0192     isLP = 1;   %% it's an LP
0193 else
0194     isLP = 0;   %% nope, it's a QP
0195 end
0196 
0197 %% default options
0198 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0199     verbose = opt.verbose;
0200 else
0201     verbose = 0;
0202 end
0203 
0204 %% split up linear constraints
0205 ieq = find( abs(u-l) <= eps );          %% equality
0206 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0207 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0208 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0209 Ae = A(ieq, :);
0210 be = u(ieq);
0211 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0212 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0213 
0214 %% grab some dimensions
0215 nlt = length(ilt);      %% number of upper bounded linear inequalities
0216 ngt = length(igt);      %% number of lower bounded linear inequalities
0217 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0218 
0219 %% mixed integer?
0220 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I'))
0221     mi = 0;
0222     vtype = repmat('C', 1, nx);
0223 else
0224     mi = 1;
0225     %% expand vtype to nx elements if necessary
0226     if length(vtype) == 1 && nx > 1
0227         vtype = char(vtype * ones(1, nx));
0228     end
0229 end
0230 %% check bounds for 'B' variables to make sure they are between 0 and 1
0231 k = find(vtype == 'B');
0232 if ~isempty(k)
0233     kk = find(xmax(k) > 1);
0234     xmax(k(kk)) = 1;
0235     kk = find(xmin(k) < 0);
0236     xmin(k(kk)) = 0;
0237 end
0238 intcon = find(vtype == 'B' | vtype == 'I');
0239 
0240 %% set up options
0241 if verbose > 1
0242     vrb = 'iter';       %% seems to be same as 'final'
0243 elseif verbose == 1
0244     vrb = 'final';
0245 else
0246     vrb = 'off';
0247 end
0248 
0249 if isLP
0250     if mi
0251         ot_opt = optimoptions('intlinprog');
0252         if ~isempty(opt) && isfield(opt, 'intlinprog_opt') && ~isempty(opt.intlinprog_opt)
0253             ot_opt = nested_struct_copy(ot_opt, opt.intlinprog_opt);
0254         end
0255     else
0256         ot_opt = optimoptions('linprog');
0257         if ~isempty(opt) && isfield(opt, 'linprog_opt') && ~isempty(opt.linprog_opt)
0258             ot_opt = nested_struct_copy(ot_opt, opt.linprog_opt);
0259         end
0260     end
0261 else
0262     if mi
0263         error('miqps_ot: MIQP problems not supported.');
0264     end
0265     ot_opt = optimoptions('quadprog');
0266     if have_fcn('quadprog_ls')
0267         ot_opt.Algorithm = 'interior-point-convex';
0268     else
0269         ot_opt.LargeScale = 'off';
0270     end
0271     if ~isempty(opt) && isfield(opt, 'quadprog_opt') && ~isempty(opt.quadprog_opt)
0272         ot_opt = nested_struct_copy(ot_opt, opt.quadprog_opt);
0273     end
0274 end
0275 ot_opt = optimoptions(ot_opt, 'Display', vrb);
0276 
0277 %% call the solver
0278 if isLP
0279     if mi
0280         [x, f, eflag, output] = ...
0281             intlinprog(c, intcon, Ai, bi, Ae, be, xmin, xmax, ot_opt);
0282         lam = [];
0283     else
0284         [x, f, eflag, output, lam] = ...
0285             linprog(c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0286     end
0287 else
0288     [x, f, eflag, output, lam] = ...
0289         quadprog(H, c, Ai, bi, Ae, be, xmin, xmax, x0, ot_opt);
0290 end
0291 
0292 %% repackage lambdas
0293 if isempty(x)
0294     x = NaN(nx, 1);
0295 end
0296 if isempty(lam) || (isempty(lam.eqlin) && isempty(lam.ineqlin) && ...
0297                     isempty(lam.lower) && isempty(lam.upper))
0298     lambda = struct( ...
0299         'mu_l', NaN(nA, 1), ...
0300         'mu_u', NaN(nA, 1), ...
0301         'lower', NaN(nx, 1), ...
0302         'upper', NaN(nx, 1) ...
0303     );
0304 else
0305     kl = find(lam.eqlin < 0);   %% lower bound binding
0306     ku = find(lam.eqlin > 0);   %% upper bound binding
0307 
0308     mu_l = zeros(nA, 1);
0309     mu_l(ieq(kl)) = -lam.eqlin(kl);
0310     mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0311     mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0312 
0313     mu_u = zeros(nA, 1);
0314     mu_u(ieq(ku)) = lam.eqlin(ku);
0315     mu_u(ilt) = lam.ineqlin(1:nlt);
0316     mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0317 
0318     lambda = struct( ...
0319         'mu_l', mu_l, ...
0320         'mu_u', mu_u, ...
0321         'lower', lam.lower(1:nx), ...
0322         'upper', lam.upper(1:nx) ...
0323     );
0324 end
0325 
0326 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
0327     if verbose
0328         fprintf('--- Integer stage complete, starting price computation stage ---\n');
0329     end
0330     if isfield(opt, 'price_stage_warn_tol') && ~isempty(opt.price_stage_warn_tol)
0331         tol = opt.price_stage_warn_tol;
0332     else
0333         tol = 1e-7;
0334     end
0335     k = intcon;
0336     x(k) = round(x(k));
0337     xmin(k) = x(k);
0338     xmax(k) = x(k);
0339     x0 = x;
0340 %     opt.linprog_opt.Algorithm = 'dual-simplex';     %% dual-simplex
0341     
0342     [x_, f_, eflag_, output_, lambda] = qps_ot(H, c, A, l, u, xmin, xmax, x0, opt);
0343 %     output
0344 %     output.message
0345     if eflag ~= eflag_
0346         error('miqps_ot: EXITFLAG from price computation stage = %d', eflag_);
0347     end
0348     if abs(f - f_)/max(abs(f), 1) > tol
0349         warning('miqps_ot: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
0350     end
0351     xn = x;
0352     xn(abs(xn)<1) = 1;
0353     [mx, k] = max(abs(x - x_) ./ xn);
0354     if mx > tol
0355         warning('miqps_ot: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
0356     end
0357     output.price_stage = output_;
0358 %     output_
0359 %     output_.message
0360 end

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