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

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