Home > matpower6.0 > mips.m

mips

PURPOSE ^

MIPS MATPOWER Interior Point Solver.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt)

DESCRIPTION ^

MIPS  MATPOWER Interior Point Solver.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIPS(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
   Primal-dual interior point method for NLP (nonlinear programming).
   Minimize a function F(X) beginning from a starting point X0, subject to
   optional linear and nonlinear constraints and variable bounds.

       min F(X)
        X

   subject to

       G(X) = 0            (nonlinear equalities)
       H(X) <= 0           (nonlinear inequalities)
       L <= A*X <= U       (linear constraints)
       XMIN <= X <= XMAX   (variable bounds)

   Inputs (all optional except F_FCN and X0):
       F_FCN : handle to function that evaluates the objective function,
           its gradients and Hessian for a given value of X. If there
           are nonlinear constraints, the Hessian information is
           provided by the HESS_FCN function passed in the 9th argument
           and is not required here. Calling syntax for this function:
               [F, DF, D2F] = F_FCN(X)
       X0 : starting value of optimization vector X
       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.
       GH_FCN : handle to function that evaluates the optional
           nonlinear constraints and their gradients for a given
           value of X. Calling syntax for this function is:
               [H, G, DH, DG] = GH_FCN(X)
       HESS_FCN : handle to function that computes the Hessian of the
           Lagrangian for given values of X, lambda and mu, where
           lambda and mu are the multipliers on the equality and
           inequality constraints, g and h, respectively. The calling
           syntax for this function is:
               LXX = HESS_FCN(X, LAM)
           where lambda = LAM.eqnonlin and mu = LAM.ineqnonlin.
       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
           linsolver ('') - linear system solver for solving update steps
               ''          = default solver (currently same as '\')
               '\'         = built-in \ operator
               'PARDISO'   = PARDISO solver (if available)
           feastol (1e-6) - termination tolerance for feasibility
               condition
           gradtol (1e-6) - termination tolerance for gradient
               condition
           comptol (1e-6) - termination tolerance for complementarity
               condition
           costtol (1e-6) - termination tolerance for cost condition
           max_it (150) - maximum number of iterations
           step_control (0) - set to 1 to enable step-size control
           sc.red_it (20) - maximum number of step-size reductions if
               step-control is on
           cost_mult (1) - cost multiplier used to scale the objective
               function for improved conditioning. Note: This value is
               also passed as the 3rd argument to the Hessian evaluation
               function so that it can appropriately scale the
               objective function term in the Hessian of the
               Lagrangian.
           xi (0.99995) - constant used in alpha updates*
           sigma (0.1) - centering parameter*
           z0 (1) - used to initialize slack variables*
           alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha
               parameter becomes smaller than this value*
           rho_min (0.95) - lower bound on rho_t*
           rho_max (1.05) - upper bound on rho_t*
           mu_threshold (1e-5) - KT multipliers smaller than this value
               for non-binding constraints are forced to zero
           max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of
               the reduced Newton step exceeds this value*
               * see the corresponding Appendix in the manual for details
       PROBLEM : The inputs can alternatively be supplied in a single
           PROBLEM struct with fields corresponding to the input arguments
           described above: f_fcn, x0, A, l, u, xmin, xmax,
                            gh_fcn, hess_fcn, opt

   Outputs:
       X : solution vector
       F : final objective function value
       EXITFLAG : exit flag
           1 = first order optimality conditions satisfied
           0 = maximum number of iterations reached
           -1 = numerically failed
       OUTPUT : output struct with fields:
           iterations - number of iterations performed
           hist - struct array with trajectories of the following:
                   feascond, gradcond, compcond, costcond, gamma,
                   stepsize, obj, alphap, alphad
           message - exit message
       LAMBDA : struct containing the Langrange and Kuhn-Tucker
           multipliers on the constraints, with fields:
           eqnonlin - nonlinear equality constraints
           ineqnonlin - nonlinear inequality constraints
           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 FMINCON
   from MathWorks' Optimization Toolbox. The main difference is that
   the linear constraints are specified with A, L, U instead of
   A, B, Aeq, Beq. The functions for evaluating the objective
   function, constraints and Hessian are identical.

   Calling syntax options:
       [x, f, exitflag, output, lambda] = ...
           mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);

       x = mips(f_fcn, x0);
       x = mips(f_fcn, x0, A, l);
       x = mips(f_fcn, x0, A, l, u);
       x = mips(f_fcn, x0, A, l, u, xmin);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
       x = mips(problem);
               where problem is a struct with fields:
                   f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt
                   all fields except 'f_fcn' and 'x0' are optional
       x = mips(...);
       [x, f] = mips(...);
       [x, f, exitflag] = mips(...);
       [x, f, exitflag, output] = mips(...);
       [x, f, exitflag, output, lambda] = mips(...);

   Example: (problem from http://en.wikipedia.org/wiki/Nonlinear_programming)
       function [f, df, d2f] = f2(x)
       f = -x(1)*x(2) - x(2)*x(3);
       if nargout > 1           %% gradient is required
           df = -[x(2); x(1)+x(3); x(2)];
           if nargout > 2       %% Hessian is required
               d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
           end                                 %% 'hess_fcn' is provided
       end
       
       function [h, g, dh, dg] = gh2(x)
       h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
       dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
       g = []; dg = [];
       
       function Lxx = hess2(x, lam, cost_mult)
       if nargin < 3, cost_mult = 1; end
       mu = lam.ineqnonlin;
       Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
               [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
       
       problem = struct( ...
           'f_fcn',    @(x)f2(x), ...
           'gh_fcn',   @(x)gh2(x), ...
           'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
           'x0',       [1; 1; 0], ...
           'opt',      struct('verbose', 2) ...
       );
       [x, f, exitflag, output, lambda] = mips(problem);

   Ported by Ray Zimmerman from C code written by H. Wang for his
   PhD dissertation:
     "On the Computation and Application of Multi-period
     Security-Constrained Optimal Power Flow for Real-time
     Electricity Market Operations", Cornell University, May 2007.

   See also:
     H. Wang, C. E. Murillo-Sanchez, R. D. Zimmerman, R. J. Thomas,
     "On Computational Issues of Market-Based Optimal Power Flow",
     IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
     pp. 1185-1193.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt)
0002 %MIPS  MATPOWER Interior Point Solver.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIPS(F_FCN, X0, A, L, U, XMIN, XMAX, GH_FCN, HESS_FCN, OPT)
0005 %   Primal-dual interior point method for NLP (nonlinear programming).
0006 %   Minimize a function F(X) beginning from a starting point X0, subject to
0007 %   optional linear and nonlinear constraints and variable bounds.
0008 %
0009 %       min F(X)
0010 %        X
0011 %
0012 %   subject to
0013 %
0014 %       G(X) = 0            (nonlinear equalities)
0015 %       H(X) <= 0           (nonlinear inequalities)
0016 %       L <= A*X <= U       (linear constraints)
0017 %       XMIN <= X <= XMAX   (variable bounds)
0018 %
0019 %   Inputs (all optional except F_FCN and X0):
0020 %       F_FCN : handle to function that evaluates the objective function,
0021 %           its gradients and Hessian for a given value of X. If there
0022 %           are nonlinear constraints, the Hessian information is
0023 %           provided by the HESS_FCN function passed in the 9th argument
0024 %           and is not required here. Calling syntax for this function:
0025 %               [F, DF, D2F] = F_FCN(X)
0026 %       X0 : starting value of optimization vector X
0027 %       A, L, U : define the optional linear constraints. Default
0028 %           values for the elements of L and U are -Inf and Inf,
0029 %           respectively.
0030 %       XMIN, XMAX : optional lower and upper bounds on the
0031 %           X variables, defaults are -Inf and Inf, respectively.
0032 %       GH_FCN : handle to function that evaluates the optional
0033 %           nonlinear constraints and their gradients for a given
0034 %           value of X. Calling syntax for this function is:
0035 %               [H, G, DH, DG] = GH_FCN(X)
0036 %       HESS_FCN : handle to function that computes the Hessian of the
0037 %           Lagrangian for given values of X, lambda and mu, where
0038 %           lambda and mu are the multipliers on the equality and
0039 %           inequality constraints, g and h, respectively. The calling
0040 %           syntax for this function is:
0041 %               LXX = HESS_FCN(X, LAM)
0042 %           where lambda = LAM.eqnonlin and mu = LAM.ineqnonlin.
0043 %       OPT : optional options structure with the following fields,
0044 %           all of which are also optional (default values shown in
0045 %           parentheses)
0046 %           verbose (0) - controls level of progress output displayed
0047 %           linsolver ('') - linear system solver for solving update steps
0048 %               ''          = default solver (currently same as '\')
0049 %               '\'         = built-in \ operator
0050 %               'PARDISO'   = PARDISO solver (if available)
0051 %           feastol (1e-6) - termination tolerance for feasibility
0052 %               condition
0053 %           gradtol (1e-6) - termination tolerance for gradient
0054 %               condition
0055 %           comptol (1e-6) - termination tolerance for complementarity
0056 %               condition
0057 %           costtol (1e-6) - termination tolerance for cost condition
0058 %           max_it (150) - maximum number of iterations
0059 %           step_control (0) - set to 1 to enable step-size control
0060 %           sc.red_it (20) - maximum number of step-size reductions if
0061 %               step-control is on
0062 %           cost_mult (1) - cost multiplier used to scale the objective
0063 %               function for improved conditioning. Note: This value is
0064 %               also passed as the 3rd argument to the Hessian evaluation
0065 %               function so that it can appropriately scale the
0066 %               objective function term in the Hessian of the
0067 %               Lagrangian.
0068 %           xi (0.99995) - constant used in alpha updates*
0069 %           sigma (0.1) - centering parameter*
0070 %           z0 (1) - used to initialize slack variables*
0071 %           alpha_min (1e-8) - returns EXITFLAG = -1 if either alpha
0072 %               parameter becomes smaller than this value*
0073 %           rho_min (0.95) - lower bound on rho_t*
0074 %           rho_max (1.05) - upper bound on rho_t*
0075 %           mu_threshold (1e-5) - KT multipliers smaller than this value
0076 %               for non-binding constraints are forced to zero
0077 %           max_stepsize (1e10) - returns EXITFLAG = -1 if the 2-norm of
0078 %               the reduced Newton step exceeds this value*
0079 %               * see the corresponding Appendix in the manual for details
0080 %       PROBLEM : The inputs can alternatively be supplied in a single
0081 %           PROBLEM struct with fields corresponding to the input arguments
0082 %           described above: f_fcn, x0, A, l, u, xmin, xmax,
0083 %                            gh_fcn, hess_fcn, opt
0084 %
0085 %   Outputs:
0086 %       X : solution vector
0087 %       F : final objective function value
0088 %       EXITFLAG : exit flag
0089 %           1 = first order optimality conditions satisfied
0090 %           0 = maximum number of iterations reached
0091 %           -1 = numerically failed
0092 %       OUTPUT : output struct with fields:
0093 %           iterations - number of iterations performed
0094 %           hist - struct array with trajectories of the following:
0095 %                   feascond, gradcond, compcond, costcond, gamma,
0096 %                   stepsize, obj, alphap, alphad
0097 %           message - exit message
0098 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0099 %           multipliers on the constraints, with fields:
0100 %           eqnonlin - nonlinear equality constraints
0101 %           ineqnonlin - nonlinear inequality constraints
0102 %           mu_l - lower (left-hand) limit on linear constraints
0103 %           mu_u - upper (right-hand) limit on linear constraints
0104 %           lower - lower bound on optimization variables
0105 %           upper - upper bound on optimization variables
0106 %
0107 %   Note the calling syntax is almost identical to that of FMINCON
0108 %   from MathWorks' Optimization Toolbox. The main difference is that
0109 %   the linear constraints are specified with A, L, U instead of
0110 %   A, B, Aeq, Beq. The functions for evaluating the objective
0111 %   function, constraints and Hessian are identical.
0112 %
0113 %   Calling syntax options:
0114 %       [x, f, exitflag, output, lambda] = ...
0115 %           mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0116 %
0117 %       x = mips(f_fcn, x0);
0118 %       x = mips(f_fcn, x0, A, l);
0119 %       x = mips(f_fcn, x0, A, l, u);
0120 %       x = mips(f_fcn, x0, A, l, u, xmin);
0121 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax);
0122 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn);
0123 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn);
0124 %       x = mips(f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt);
0125 %       x = mips(problem);
0126 %               where problem is a struct with fields:
0127 %                   f_fcn, x0, A, l, u, xmin, xmax, gh_fcn, hess_fcn, opt
0128 %                   all fields except 'f_fcn' and 'x0' are optional
0129 %       x = mips(...);
0130 %       [x, f] = mips(...);
0131 %       [x, f, exitflag] = mips(...);
0132 %       [x, f, exitflag, output] = mips(...);
0133 %       [x, f, exitflag, output, lambda] = mips(...);
0134 %
0135 %   Example: (problem from http://en.wikipedia.org/wiki/Nonlinear_programming)
0136 %       function [f, df, d2f] = f2(x)
0137 %       f = -x(1)*x(2) - x(2)*x(3);
0138 %       if nargout > 1           %% gradient is required
0139 %           df = -[x(2); x(1)+x(3); x(2)];
0140 %           if nargout > 2       %% Hessian is required
0141 %               d2f = -[0 1 0; 1 0 1; 0 1 0];   %% actually not used since
0142 %           end                                 %% 'hess_fcn' is provided
0143 %       end
0144 %
0145 %       function [h, g, dh, dg] = gh2(x)
0146 %       h = [ 1 -1 1; 1 1 1] * x.^2 + [-2; -10];
0147 %       dh = 2 * [x(1) x(1); -x(2) x(2); x(3) x(3)];
0148 %       g = []; dg = [];
0149 %
0150 %       function Lxx = hess2(x, lam, cost_mult)
0151 %       if nargin < 3, cost_mult = 1; end
0152 %       mu = lam.ineqnonlin;
0153 %       Lxx = cost_mult * [0 -1 0; -1 0 -1; 0 -1 0] + ...
0154 %               [2*[1 1]*mu 0 0; 0 2*[-1 1]*mu 0; 0 0 2*[1 1]*mu];
0155 %
0156 %       problem = struct( ...
0157 %           'f_fcn',    @(x)f2(x), ...
0158 %           'gh_fcn',   @(x)gh2(x), ...
0159 %           'hess_fcn', @(x, lam, cost_mult)hess2(x, lam, cost_mult), ...
0160 %           'x0',       [1; 1; 0], ...
0161 %           'opt',      struct('verbose', 2) ...
0162 %       );
0163 %       [x, f, exitflag, output, lambda] = mips(problem);
0164 %
0165 %   Ported by Ray Zimmerman from C code written by H. Wang for his
0166 %   PhD dissertation:
0167 %     "On the Computation and Application of Multi-period
0168 %     Security-Constrained Optimal Power Flow for Real-time
0169 %     Electricity Market Operations", Cornell University, May 2007.
0170 %
0171 %   See also:
0172 %     H. Wang, C. E. Murillo-Sanchez, R. D. Zimmerman, R. J. Thomas,
0173 %     "On Computational Issues of Market-Based Optimal Power Flow",
0174 %     IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
0175 %     pp. 1185-1193.
0176 
0177 %   MIPS
0178 %   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
0179 %   by Ray Zimmerman, PSERC Cornell
0180 %
0181 %   This file is part of MIPS.
0182 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0183 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0184 
0185 %%----- input argument handling  -----
0186 %% gather inputs
0187 if nargin == 1 && isstruct(f_fcn)       %% problem struct
0188     p = f_fcn;
0189     f_fcn = p.f_fcn;
0190     x0 = p.x0;
0191     nx = size(x0, 1);       %% number of optimization variables
0192     if isfield(p, 'opt'),       opt = p.opt;            else,   opt = [];       end
0193     if isfield(p, 'hess_fcn'),  hess_fcn = p.hess_fcn;  else,   hess_fcn = '';  end
0194     if isfield(p, 'gh_fcn'),    gh_fcn = p.gh_fcn;      else,   gh_fcn = '';    end
0195     if isfield(p, 'xmax'),      xmax = p.xmax;          else,   xmax = [];      end
0196     if isfield(p, 'xmin'),      xmin = p.xmin;          else,   xmin = [];      end
0197     if isfield(p, 'u'),         u = p.u;                else,   u = [];         end
0198     if isfield(p, 'l'),         l = p.l;                else,   l = [];         end
0199     if isfield(p, 'A'),         A = p.A;                else,   A=sparse(0,nx); end
0200 else                                    %% individual args
0201     nx = size(x0, 1);       %% number of optimization variables
0202     if nargin < 10
0203         opt = [];
0204         if nargin < 9
0205             hess_fcn = '';
0206             if nargin < 8
0207                 gh_fcn = '';
0208                 if nargin < 7
0209                     xmax = [];
0210                     if nargin < 6
0211                         xmin = [];
0212                         if nargin < 5
0213                             u = [];
0214                             if nargin < 4
0215                                 l = [];
0216                                 A = sparse(0,nx);
0217                             end
0218                         end
0219                     end
0220                 end
0221             end
0222         end
0223     end
0224 end
0225 %% set default argument values if missing
0226 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0227                                  (isempty(u) || all(u == Inf)))
0228     A = sparse(0,nx);           %% no limits => no linear constraints
0229 end
0230 nA = size(A, 1);                %% number of original linear constraints
0231 if isempty(u)                   %% By default, linear inequalities are ...
0232     u = Inf(nA, 1);             %% ... unbounded above and ...
0233 end
0234 if isempty(l)
0235     l = -Inf(nA, 1);            %% ... unbounded below.
0236 end
0237 if isempty(xmin)                %% By default, optimization variables are ...
0238     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0239 end
0240 if isempty(xmax)
0241     xmax = Inf(nx, 1);          %% ... unbounded above.
0242 end
0243 if isempty(gh_fcn)
0244     nonlinear = false;          %% no nonlinear constraints present
0245     gn = []; hn = [];
0246 else
0247     nonlinear = true;           %% we have some nonlinear constraints
0248 end
0249 
0250 %% default options
0251 if isempty(opt)
0252     opt = struct('verbose', 0);
0253 end
0254 if ~isfield(opt, 'linsolver') || isempty(opt.linsolver)
0255     opt.linsolver = '';
0256 end
0257 if ~isfield(opt, 'feastol') || isempty(opt.feastol) || opt.feastol == 0
0258     opt.feastol = 1e-6;
0259 end
0260 if ~isfield(opt, 'gradtol') || isempty(opt.gradtol) || opt.gradtol == 0
0261     opt.gradtol = 1e-6;
0262 end
0263 if ~isfield(opt, 'comptol') || isempty(opt.comptol) || opt.comptol == 0
0264     opt.comptol = 1e-6;
0265 end
0266 if ~isfield(opt, 'costtol') || isempty(opt.costtol) || opt.costtol == 0
0267     opt.costtol = 1e-6;
0268 end
0269 if ~isfield(opt, 'max_it') || isempty(opt.max_it)
0270     opt.max_it = 150;
0271 end
0272 if ~isfield(opt, 'sc') || ~isfield(opt.sc, 'red_it') || isempty(opt.sc.red_it)
0273     opt.sc.red_it = 20;
0274 end
0275 if ~isfield(opt, 'step_control') || isempty(opt.step_control)
0276     opt.step_control = 0;
0277 end
0278 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0279     opt.cost_mult = 1;
0280 end
0281 if ~isfield(opt, 'verbose') || isempty(opt.verbose)
0282     opt.verbose = 0;
0283 end
0284 %% algorithm constants
0285 if ~isfield(opt, 'xi') || isempty(opt.xi)
0286     opt.xi = 0.99995;           %% OPT_IPM_PHI
0287 end
0288 if ~isfield(opt, 'sigma') || isempty(opt.sigma)
0289     opt.sigma = 0.1;            %% OPT_IPM_SIGMA, centering parameter
0290 end
0291 if ~isfield(opt, 'z0') || isempty(opt.z0)
0292     opt.z0 = 1;                 %% OPT_IPM_INIT_SLACK
0293 end
0294 if ~isfield(opt, 'alpha_min') || isempty(opt.alpha_min)
0295     opt.alpha_min = 1e-8;       %% OPT_AP_AD_MIN
0296 end
0297 if ~isfield(opt, 'rho_min') || isempty(opt.rho_min)
0298     opt.rho_min = 0.95;         %% OPT_IPM_QUAD_LOWTHRESH
0299 end
0300 if ~isfield(opt, 'rho_max') || isempty(opt.rho_max)
0301     opt.rho_max = 1.05;         %% OPT_IPM_QUAD_HIGHTHRESH
0302 end
0303 if ~isfield(opt, 'mu_threshold') || isempty(opt.mu_threshold)
0304     opt.mu_threshold = 1e-5;    %% SCOPF_MULTIPLIERS_FILTER_THRESH
0305 end
0306 if ~isfield(opt, 'max_stepsize') || isempty(opt.max_stepsize)
0307     opt.max_stepsize = 1e10;
0308 end
0309 
0310 %% initialize history
0311 hist(opt.max_it+1) = struct('feascond', 0, 'gradcond', 0, 'compcond', 0, ...
0312     'costcond', 0, 'gamma', 0, 'stepsize', 0, 'obj', 0, ...
0313     'alphap', 0, 'alphad', 0);
0314 
0315 %%-----  set up problem  -----
0316 %% constants
0317 [xi, sigma, z0, alpha_min, rho_min, rho_max, mu_threshold, max_stepsize] = ...
0318     deal(opt.xi, opt.sigma, opt.z0, opt.alpha_min, ...
0319         opt.rho_min, opt.rho_max, opt.mu_threshold, opt.max_stepsize);
0320 if xi >= 1 || xi < 0.5
0321     error('mips: opt.xi (%g) must be a number slightly less than 1', opt.xi);
0322 end
0323 if sigma > 1 || sigma <= 0
0324     error('mips: opt.sigma (%g) must be a number between 0 and 1', opt.sigma);
0325 end
0326 
0327 %% initialize
0328 i = 0;                      %% iteration counter
0329 converged = 0;              %% flag
0330 eflag = 0;                  %% exit flag
0331 
0332 %% add var limits to linear constraints
0333 AA = [speye(nx); A];
0334 ll = [xmin; l];
0335 uu = [xmax; u];
0336 
0337 %% split up linear constraints
0338 ieq = find( abs(uu-ll) <= eps );            %% equality
0339 igt = find( uu >=  1e10 & ll > -1e10 );     %% greater than, unbounded above
0340 ilt = find( ll <= -1e10 & uu <  1e10 );     %% less than, unbounded below
0341 ibx = find( (abs(uu-ll) > eps) & (uu < 1e10) & (ll > -1e10) );
0342 Ae = AA(ieq, :);
0343 be = uu(ieq, 1);
0344 Ai  = [ AA(ilt, :); -AA(igt, :); AA(ibx, :); -AA(ibx, :) ];
0345 bi  = [ uu(ilt, 1); -ll(igt, 1); uu(ibx, 1); -ll(ibx, 1) ];
0346 
0347 %% evaluate cost f(x0) and constraints g(x0), h(x0)
0348 x = x0;
0349 [f, df] = f_fcn(x);             %% cost
0350 f = f * opt.cost_mult;
0351 df = df * opt.cost_mult;
0352 if nonlinear
0353     [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0354     h = [hn; Ai * x - bi];          %% inequality constraints
0355     g = [gn; Ae * x - be];          %% equality constraints
0356     dh = [dhn Ai'];                 %% 1st derivative of inequalities
0357     dg = [dgn Ae'];                 %% 1st derivative of equalities
0358 else
0359     h = Ai * x - bi;                %% inequality constraints
0360     g = Ae * x - be;                %% equality constraints
0361     dh = Ai';                       %% 1st derivative of inequalities
0362     dg = Ae';                       %% 1st derivative of equalities
0363 end
0364 
0365 %% grab some dimensions
0366 neq = size(g, 1);           %% number of equality constraints
0367 niq = size(h, 1);           %% number of inequality constraints
0368 neqnln = size(gn, 1);       %% number of nonlinear equality constraints
0369 niqnln = size(hn, 1);       %% number of nonlinear inequality constraints
0370 nlt = length(ilt);          %% number of upper bounded linear inequalities
0371 ngt = length(igt);          %% number of lower bounded linear inequalities
0372 nbx = length(ibx);          %% number of doubly bounded linear inequalities
0373 
0374 %% initialize gamma, lam, mu, z, e
0375 gamma = 1;                  %% barrier coefficient, r in Harry's code
0376 lam = zeros(neq, 1);
0377 z   = z0 * ones(niq, 1);
0378 mu  = z;
0379 k = find(h < -z0);
0380 z(k) = -h(k);
0381 k = find(gamma / z > z0);   %% (seems k is always empty if gamma = z0 = 1)
0382 if ~isempty(k)
0383     mu(k) = gamma / z(k);
0384 end
0385 e = ones(niq, 1);
0386 
0387 %% check tolerance
0388 f0 = f;
0389 if opt.step_control
0390     L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0391 end
0392 Lx = df + dg * lam + dh * mu;
0393 if isempty(h)
0394     maxh = zeros(1,0);
0395 else
0396     maxh = max(h);
0397 end
0398 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0399 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0400 compcond = (z' * mu) / (1 + norm(x, Inf));
0401 costcond = abs(f - f0) / (1 + abs(f0));
0402 %% save history
0403 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0404     'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0405     'stepsize', 0, 'obj', f/opt.cost_mult, 'alphap', 0, 'alphad', 0);
0406 if strcmp(upper(opt.linsolver), 'PARDISO')
0407     ls = 'PARDISO';
0408     if ~have_fcn('pardiso')
0409         warning('mips: PARDISO linear solver not available, using default');
0410         opt.linsolver = '';
0411         ls = 'built-in';
0412     end
0413 else
0414     ls = 'built-in';
0415 end
0416 if opt.verbose
0417     if opt.step_control, s = '-sc'; else, s = ''; end
0418     v = mipsver('all');
0419     fprintf('MATPOWER Interior Point Solver -- MIPS%s, Version %s, %s\n (using %s linear solver)', ...
0420         s, v.Version, v.Date, ls);
0421     if opt.verbose > 1
0422         fprintf('\n it    objective   step size   feascond     gradcond     compcond     costcond  ');
0423         fprintf('\n----  ------------ --------- ------------ ------------ ------------ ------------');
0424         fprintf('\n%3d  %12.8g %10s %12g %12g %12g %12g', ...
0425             i, f/opt.cost_mult, '', feascond, gradcond, compcond, costcond);
0426     end
0427 end
0428 if feascond < opt.feastol && gradcond < opt.gradtol && ...
0429                 compcond < opt.comptol && costcond < opt.costtol
0430     converged = 1;
0431     if opt.verbose
0432         fprintf('\nConverged!\n');
0433     end
0434 end
0435 
0436 %%-----  do Newton iterations  -----
0437 while (~converged && i < opt.max_it)
0438     %% update iteration counter
0439     i = i + 1;
0440 
0441     %% compute update step
0442     lambda = struct('eqnonlin', lam(1:neqnln), 'ineqnonlin', mu(1:niqnln));
0443     if nonlinear
0444         if isempty(hess_fcn)
0445             fprintf('mips: Hessian evaluation via finite differences not yet implemented.\n       Please provide your own hessian evaluation function.');
0446         end
0447         Lxx = hess_fcn(x, lambda, opt.cost_mult);
0448     else
0449         [f_, df_, d2f] = f_fcn(x);      %% cost
0450         Lxx = d2f * opt.cost_mult;
0451     end
0452     zinvdiag = sparse(1:niq, 1:niq, 1 ./ z, niq, niq);
0453     mudiag = sparse(1:niq, 1:niq, mu, niq, niq);
0454     dh_zinv = dh * zinvdiag;
0455     M = Lxx + dh_zinv * mudiag * dh';
0456     N = Lx + dh_zinv * (mudiag * h + gamma * e);
0457     dxdlam = mplinsolve([M dg; dg' sparse(neq, neq)], [-N; -g], opt.linsolver, []);
0458 %     AAA = [
0459 %         M  dg;
0460 %         dg'  sparse(neq, neq)
0461 %     ];
0462 %     rc = 1/condest(AAA);
0463 %     if rc < 1e-22
0464 %         fprintf('my RCOND = %g\n', rc);
0465 %         n = size(AAA, 1);
0466 %         AAA = AAA + 1e-3 * speye(n,n);
0467 %     end
0468 %     bbb = [-N; -g];
0469 %     dxdlam = AAA \ bbb;
0470     if any(isnan(dxdlam)) || norm(dxdlam) > max_stepsize
0471         if opt.verbose
0472             fprintf('\nNumerically Failed\n');
0473         end
0474         eflag = -1;
0475         break;
0476     end
0477     dx = dxdlam(1:nx, 1);
0478     dlam = dxdlam(nx+(1:neq), 1);
0479     dz = -h - z - dh' * dx;
0480     dmu = -mu + zinvdiag *(gamma*e - mudiag * dz);
0481 
0482     %% optional step-size control
0483     sc = 0;
0484     if opt.step_control
0485         x1 = x + dx;
0486 
0487         %% evaluate cost, constraints, derivatives at x1
0488         [f1, df1] = f_fcn(x1);          %% cost
0489         f1 = f1 * opt.cost_mult;
0490         df1 = df1 * opt.cost_mult;
0491         if nonlinear
0492             [hn1, gn1, dhn1, dgn1] = gh_fcn(x1); %% nonlinear constraints
0493             h1 = [hn1; Ai * x1 - bi];       %% inequality constraints
0494             g1 = [gn1; Ae * x1 - be];       %% equality constraints
0495             dh1 = [dhn1 Ai'];               %% 1st derivative of inequalities
0496             dg1 = [dgn1 Ae'];               %% 1st derivative of equalities
0497         else
0498             h1 = Ai * x1 - bi;              %% inequality constraints
0499             g1 = Ae * x1 - be;              %% equality constraints
0500             dh1 = dh;                       %% 1st derivative of inequalities
0501             dg1 = dg;                       %% 1st derivative of equalities
0502         end
0503 
0504         %% check tolerance
0505         Lx1 = df1 + dg1 * lam + dh1 * mu;
0506         if isempty(h1)
0507             maxh1 = zeros(1,0);
0508         else
0509             maxh1 = max(h1);
0510         end
0511         feascond1 = max([norm(g1, Inf), maxh1]) / (1 + max([ norm(x1, Inf), norm(z, Inf) ]));
0512         gradcond1 = norm(Lx1, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0513 
0514         if feascond1 > feascond && gradcond1 > gradcond
0515             sc = 1;
0516         end
0517     end
0518     if sc
0519         alpha = 1;
0520         for j = 1:opt.sc.red_it
0521             dx1 = alpha * dx;
0522             x1 = x + dx1;
0523             f1 = f_fcn(x1);                 %% cost
0524             f1 = f1 * opt.cost_mult;
0525             if nonlinear
0526                 [hn1, gn1] = gh_fcn(x1);    %% nonlinear constraints
0527                 h1 = [hn1; Ai * x1 - bi];   %% inequality constraints
0528                 g1 = [gn1; Ae * x1 - be];   %% equality constraints
0529             else
0530                 h1 = Ai * x1 - bi;          %% inequality constraints
0531                 g1 = Ae * x1 - be;          %% equality constraints
0532             end
0533             L1 = f1 + lam' * g1 + mu' * (h1+z) - gamma * sum(log(z));
0534             if opt.verbose > 2
0535                 fprintf('\n   %3d            %10g', -j, norm(dx1));
0536             end
0537             rho = (L1 - L) / (Lx' * dx1 + 0.5 * dx1' * Lxx * dx1);
0538             if rho > rho_min && rho < rho_max
0539                 break;
0540             else
0541                 alpha = alpha / 2;
0542             end
0543         end
0544         dx = alpha * dx;
0545         dz = alpha * dz;
0546         dlam = alpha * dlam;
0547         dmu = alpha * dmu;
0548     end
0549 
0550     %% do the update
0551     k = find(dz < 0);
0552     if isempty(k)
0553         alphap = 1;
0554     else
0555         alphap = min( [xi * min(z(k) ./ -dz(k)) 1] );
0556     end
0557     k = find(dmu < 0);
0558     if isempty(k)
0559         alphad = 1;
0560     else
0561         alphad = min( [xi * min(mu(k) ./ -dmu(k)) 1] );
0562     end
0563     x = x + alphap * dx;
0564     z = z + alphap * dz;
0565     lam = lam + alphad * dlam;
0566     mu  = mu  + alphad * dmu;
0567     if niq > 0
0568         gamma = sigma * (z' * mu) / niq;
0569     end
0570 
0571     %% evaluate cost, constraints, derivatives
0572     [f, df] = f_fcn(x);                 %% cost
0573     f = f * opt.cost_mult;
0574     df = df * opt.cost_mult;
0575     if nonlinear
0576         [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0577         h = [hn; Ai * x - bi];          %% inequality constraints
0578         g = [gn; Ae * x - be];          %% equality constraints
0579         dh = [dhn Ai'];                 %% 1st derivative of inequalities
0580         dg = [dgn Ae'];                 %% 1st derivative of equalities
0581     else
0582         h = Ai * x - bi;                %% inequality constraints
0583         g = Ae * x - be;                %% equality constraints
0584         %% 1st derivatives are constant, still dh = Ai', dg = Ae'
0585     end
0586 
0587     %% check tolerance
0588     Lx = df + dg * lam + dh * mu;
0589     if isempty(h)
0590         maxh = zeros(1,0);
0591     else
0592         maxh = max(h);
0593     end
0594     feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0595     gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0596     compcond = (z' * mu) / (1 + norm(x, Inf));
0597     costcond = abs(f - f0) / (1 + abs(f0));
0598     %% save history
0599     hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0600         'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0601         'stepsize', norm(dx), 'obj', f/opt.cost_mult, ...
0602         'alphap', alphap, 'alphad', alphad);
0603 
0604     if opt.verbose > 1
0605         fprintf('\n%3d  %12.8g %10.5g %12g %12g %12g %12g', ...
0606             i, f/opt.cost_mult, norm(dx), feascond, gradcond, compcond, costcond);
0607     end
0608     if feascond < opt.feastol && gradcond < opt.gradtol && ...
0609                     compcond < opt.comptol && costcond < opt.costtol
0610         converged = 1;
0611         if opt.verbose
0612             fprintf('\nConverged!\n');
0613         end
0614     else
0615         if any(isnan(x)) || alphap < alpha_min || alphad < alpha_min || ...
0616                 gamma < eps || gamma > 1/eps
0617             if opt.verbose
0618                 fprintf('\nNumerically Failed\n');
0619             end
0620             eflag = -1;
0621             break;
0622         end
0623         f0 = f;
0624         if opt.step_control
0625             L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0626         end
0627     end
0628 end
0629 
0630 if opt.verbose
0631     if ~converged
0632         fprintf('\nDid not converge in %d iterations.\n', i);
0633     end
0634 end
0635 
0636 %%-----  package up results  -----
0637 hist = hist(1:i+1);
0638 if eflag ~= -1
0639     eflag = converged;
0640 end
0641 output = struct('iterations', i, 'hist', hist, 'message', '');
0642 if eflag == 0
0643     output.message = 'Did not converge';
0644 elseif eflag == 1
0645     output.message = 'Converged';
0646 elseif eflag == -1
0647     output.message = 'Numerically failed';
0648 else
0649     output.message = 'Please hang up and dial again';
0650 end
0651 
0652 %% zero out multipliers on non-binding constraints
0653 mu(h < -opt.feastol & mu < mu_threshold) = 0;
0654 
0655 %% un-scale cost and prices
0656 f   = f   / opt.cost_mult;
0657 lam = lam / opt.cost_mult;
0658 mu  = mu  / opt.cost_mult;
0659 
0660 %% re-package multipliers into struct
0661 lam_lin = lam((neqnln+1):neq);              %% lambda for linear constraints
0662 mu_lin  = mu((niqnln+1):niq);               %% mu for linear constraints
0663 kl = find(lam_lin < 0);                     %% lower bound binding
0664 ku = find(lam_lin > 0);                     %% upper bound binding
0665 
0666 mu_l = zeros(nx+nA, 1);
0667 mu_l(ieq(kl)) = -lam_lin(kl);
0668 mu_l(igt) = mu_lin(nlt+(1:ngt));
0669 mu_l(ibx) = mu_lin(nlt+ngt+nbx+(1:nbx));
0670 
0671 mu_u = zeros(nx+nA, 1);
0672 mu_u(ieq(ku)) = lam_lin(ku);
0673 mu_u(ilt) = mu_lin(1:nlt);
0674 mu_u(ibx) = mu_lin(nlt+ngt+(1:nbx));
0675 
0676 fields = { ...
0677     'mu_l', mu_l((nx+1):end), ...
0678     'mu_u', mu_u((nx+1):end), ...
0679     'lower', mu_l(1:nx), ...
0680     'upper', mu_u(1:nx) ...
0681 };
0682 
0683 if niqnln > 0
0684     fields = { ...
0685         'ineqnonlin', mu(1:niqnln), ...
0686         fields{:} ...
0687     };
0688 end
0689 if neqnln > 0
0690     fields = { ...
0691         'eqnonlin', lam(1:neqnln), ...
0692         fields{:} ...
0693     };
0694 end
0695 
0696 lambda = struct(fields{:});
0697 
0698 % lambda = struct( ...
0699 %     'eqnonlin', lam(1:neqnln), ...
0700 %     'ineqnonlin', mu(1:niqnln), ...
0701 %     'mu_l', mu_l((nx+1):end), ...
0702 %     'mu_u', mu_u((nx+1):end), ...
0703 %     'lower', mu_l(1:nx), ...
0704 %     'upper', mu_u(1:nx) );

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005