Home > matpower5.1 > mips.m

mips

PURPOSE ^

MIPS MATLAB 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  MATLAB 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  MATLAB 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-2015 by Power System Engineering Research Center (PSERC)
0179 %   by Ray Zimmerman, PSERC Cornell
0180 %
0181 %   $Id: mips.m 2661 2015-03-20 17:02:46Z ray $
0182 %
0183 %   This file is part of MIPS.
0184 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0185 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0186 
0187 %%----- input argument handling  -----
0188 %% gather inputs
0189 if nargin == 1 && isstruct(f_fcn)       %% problem struct
0190     p = f_fcn;
0191     f_fcn = p.f_fcn;
0192     x0 = p.x0;
0193     nx = size(x0, 1);       %% number of optimization variables
0194     if isfield(p, 'opt'),       opt = p.opt;            else,   opt = [];       end
0195     if isfield(p, 'hess_fcn'),  hess_fcn = p.hess_fcn;  else,   hess_fcn = '';  end
0196     if isfield(p, 'gh_fcn'),    gh_fcn = p.gh_fcn;      else,   gh_fcn = '';    end
0197     if isfield(p, 'xmax'),      xmax = p.xmax;          else,   xmax = [];      end
0198     if isfield(p, 'xmin'),      xmin = p.xmin;          else,   xmin = [];      end
0199     if isfield(p, 'u'),         u = p.u;                else,   u = [];         end
0200     if isfield(p, 'l'),         l = p.l;                else,   l = [];         end
0201     if isfield(p, 'A'),         A = p.A;                else,   A=sparse(0,nx); end
0202 else                                    %% individual args
0203     nx = size(x0, 1);       %% number of optimization variables
0204     if nargin < 10
0205         opt = [];
0206         if nargin < 9
0207             hess_fcn = '';
0208             if nargin < 8
0209                 gh_fcn = '';
0210                 if nargin < 7
0211                     xmax = [];
0212                     if nargin < 6
0213                         xmin = [];
0214                         if nargin < 5
0215                             u = [];
0216                             if nargin < 4
0217                                 l = [];
0218                                 A = sparse(0,nx);
0219                             end
0220                         end
0221                     end
0222                 end
0223             end
0224         end
0225     end
0226 end
0227 %% set default argument values if missing
0228 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0229                                  (isempty(u) || all(u == Inf)))
0230     A = sparse(0,nx);           %% no limits => no linear constraints
0231 end
0232 nA = size(A, 1);                %% number of original linear constraints
0233 if isempty(u)                   %% By default, linear inequalities are ...
0234     u = Inf(nA, 1);             %% ... unbounded above and ...
0235 end
0236 if isempty(l)
0237     l = -Inf(nA, 1);            %% ... unbounded below.
0238 end
0239 if isempty(xmin)                %% By default, optimization variables are ...
0240     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0241 end
0242 if isempty(xmax)
0243     xmax = Inf(nx, 1);          %% ... unbounded above.
0244 end
0245 if isempty(gh_fcn)
0246     nonlinear = false;          %% no nonlinear constraints present
0247     gn = []; hn = [];
0248 else
0249     nonlinear = true;           %% we have some nonlinear constraints
0250 end
0251 
0252 %% default options
0253 if isempty(opt)
0254     opt = struct('verbose', 0);
0255 end
0256 if ~isfield(opt, 'linsolver') || isempty(opt.linsolver)
0257     opt.linsolver = '';
0258 end
0259 if ~isfield(opt, 'feastol') || isempty(opt.feastol)
0260     opt.feastol = 1e-6;
0261 end
0262 if ~isfield(opt, 'gradtol') || isempty(opt.gradtol)
0263     opt.gradtol = 1e-6;
0264 end
0265 if ~isfield(opt, 'comptol') || isempty(opt.comptol)
0266     opt.comptol = 1e-6;
0267 end
0268 if ~isfield(opt, 'costtol') || isempty(opt.costtol)
0269     opt.costtol = 1e-6;
0270 end
0271 if ~isfield(opt, 'max_it') || isempty(opt.max_it)
0272     opt.max_it = 150;
0273 end
0274 if ~isfield(opt, 'sc') || ~isfield(opt.sc, 'red_it') || isempty(opt.sc.red_it)
0275     opt.sc.red_it = 20;
0276 end
0277 if ~isfield(opt, 'step_control') || isempty(opt.step_control)
0278     opt.step_control = 0;
0279 end
0280 if ~isfield(opt, 'cost_mult') || isempty(opt.cost_mult)
0281     opt.cost_mult = 1;
0282 end
0283 if ~isfield(opt, 'verbose') || isempty(opt.verbose)
0284     opt.verbose = 0;
0285 end
0286 %% algorithm constants
0287 if ~isfield(opt, 'xi') || isempty(opt.xi)
0288     opt.xi = 0.99995;           %% OPT_IPM_PHI
0289 end
0290 if ~isfield(opt, 'sigma') || isempty(opt.sigma)
0291     opt.sigma = 0.1;            %% OPT_IPM_SIGMA, centering parameter
0292 end
0293 if ~isfield(opt, 'z0') || isempty(opt.z0)
0294     opt.z0 = 1;                 %% OPT_IPM_INIT_SLACK
0295 end
0296 if ~isfield(opt, 'alpha_min') || isempty(opt.alpha_min)
0297     opt.alpha_min = 1e-8;       %% OPT_AP_AD_MIN
0298 end
0299 if ~isfield(opt, 'rho_min') || isempty(opt.rho_min)
0300     opt.rho_min = 0.95;         %% OPT_IPM_QUAD_LOWTHRESH
0301 end
0302 if ~isfield(opt, 'rho_max') || isempty(opt.rho_max)
0303     opt.rho_max = 1.05;         %% OPT_IPM_QUAD_HIGHTHRESH
0304 end
0305 if ~isfield(opt, 'mu_threshold') || isempty(opt.mu_threshold)
0306     opt.mu_threshold = 1e-5;    %% SCOPF_MULTIPLIERS_FILTER_THRESH
0307 end
0308 if ~isfield(opt, 'max_stepsize') || isempty(opt.max_stepsize)
0309     opt.max_stepsize = 1e10;
0310 end
0311 
0312 %% initialize history
0313 hist(opt.max_it+1) = struct('feascond', 0, 'gradcond', 0, 'compcond', 0, ...
0314     'costcond', 0, 'gamma', 0, 'stepsize', 0, 'obj', 0, ...
0315     'alphap', 0, 'alphad', 0);
0316 
0317 %%-----  set up problem  -----
0318 %% constants
0319 [xi, sigma, z0, alpha_min, rho_min, rho_max, mu_threshold, max_stepsize] = ...
0320     deal(opt.xi, opt.sigma, opt.z0, opt.alpha_min, ...
0321         opt.rho_min, opt.rho_max, opt.mu_threshold, opt.max_stepsize);
0322 if xi >= 1 || xi < 0.5
0323     error('mips: opt.xi (%g) must be a number slightly less than 1', opt.xi);
0324 end
0325 if sigma > 1 || sigma <= 0
0326     error('mips: opt.sigma (%g) must be a number between 0 and 1', opt.sigma);
0327 end
0328 
0329 %% initialize
0330 i = 0;                      %% iteration counter
0331 converged = 0;              %% flag
0332 eflag = 0;                  %% exit flag
0333 
0334 %% add var limits to linear constraints
0335 AA = [speye(nx); A];
0336 ll = [xmin; l];
0337 uu = [xmax; u];
0338 
0339 %% split up linear constraints
0340 ieq = find( abs(uu-ll) <= eps );            %% equality
0341 igt = find( uu >=  1e10 & ll > -1e10 );     %% greater than, unbounded above
0342 ilt = find( ll <= -1e10 & uu <  1e10 );     %% less than, unbounded below
0343 ibx = find( (abs(uu-ll) > eps) & (uu < 1e10) & (ll > -1e10) );
0344 Ae = AA(ieq, :);
0345 be = uu(ieq, 1);
0346 Ai  = [ AA(ilt, :); -AA(igt, :); AA(ibx, :); -AA(ibx, :) ];
0347 bi  = [ uu(ilt, 1); -ll(igt, 1); uu(ibx, 1); -ll(ibx, 1) ];
0348 
0349 %% evaluate cost f(x0) and constraints g(x0), h(x0)
0350 x = x0;
0351 [f, df] = f_fcn(x);             %% cost
0352 f = f * opt.cost_mult;
0353 df = df * opt.cost_mult;
0354 if nonlinear
0355     [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0356     h = [hn; Ai * x - bi];          %% inequality constraints
0357     g = [gn; Ae * x - be];          %% equality constraints
0358     dh = [dhn Ai'];                 %% 1st derivative of inequalities
0359     dg = [dgn Ae'];                 %% 1st derivative of equalities
0360 else
0361     h = Ai * x - bi;                %% inequality constraints
0362     g = Ae * x - be;                %% equality constraints
0363     dh = Ai';                       %% 1st derivative of inequalities
0364     dg = Ae';                       %% 1st derivative of equalities
0365 end
0366 
0367 %% grab some dimensions
0368 neq = size(g, 1);           %% number of equality constraints
0369 niq = size(h, 1);           %% number of inequality constraints
0370 neqnln = size(gn, 1);       %% number of nonlinear equality constraints
0371 niqnln = size(hn, 1);       %% number of nonlinear inequality constraints
0372 nlt = length(ilt);          %% number of upper bounded linear inequalities
0373 ngt = length(igt);          %% number of lower bounded linear inequalities
0374 nbx = length(ibx);          %% number of doubly bounded linear inequalities
0375 
0376 %% initialize gamma, lam, mu, z, e
0377 gamma = 1;                  %% barrier coefficient, r in Harry's code
0378 lam = zeros(neq, 1);
0379 z   = z0 * ones(niq, 1);
0380 mu  = z;
0381 k = find(h < -z0);
0382 z(k) = -h(k);
0383 k = find(gamma / z > z0);   %% (seems k is always empty if gamma = z0 = 1)
0384 if ~isempty(k)
0385     mu(k) = gamma / z(k);
0386 end
0387 e = ones(niq, 1);
0388 
0389 %% check tolerance
0390 f0 = f;
0391 if opt.step_control
0392     L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0393 end
0394 Lx = df + dg * lam + dh * mu;
0395 if isempty(h)
0396     maxh = zeros(1,0);
0397 else
0398     maxh = max(h);
0399 end
0400 feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0401 gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0402 compcond = (z' * mu) / (1 + norm(x, Inf));
0403 costcond = abs(f - f0) / (1 + abs(f0));
0404 %% save history
0405 hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0406     'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0407     'stepsize', 0, 'obj', f/opt.cost_mult, 'alphap', 0, 'alphad', 0);
0408 if strcmp(upper(opt.linsolver), 'PARDISO')
0409     ls = 'PARDISO';
0410     if ~have_fcn('pardiso')
0411         warning('mips: PARDISO linear solver not available, using default');
0412         opt.linsolver = '';
0413         ls = 'built-in';
0414     end
0415 else
0416     ls = 'built-in';
0417 end
0418 if opt.verbose
0419     if opt.step_control, s = '-sc'; else, s = ''; end
0420     v = mipsver('all');
0421     fprintf('MATLAB Interior Point Solver -- MIPS%s, Version %s, %s\n (using %s linear solver)', ...
0422         s, v.Version, v.Date, ls);
0423     if opt.verbose > 1
0424         fprintf('\n it    objective   step size   feascond     gradcond     compcond     costcond  ');
0425         fprintf('\n----  ------------ --------- ------------ ------------ ------------ ------------');
0426         fprintf('\n%3d  %12.8g %10s %12g %12g %12g %12g', ...
0427             i, f/opt.cost_mult, '', feascond, gradcond, compcond, costcond);
0428     end
0429 end
0430 if feascond < opt.feastol && gradcond < opt.gradtol && ...
0431                 compcond < opt.comptol && costcond < opt.costtol
0432     converged = 1;
0433     if opt.verbose
0434         fprintf('\nConverged!\n');
0435     end
0436 end
0437 
0438 %%-----  do Newton iterations  -----
0439 while (~converged && i < opt.max_it)
0440     %% update iteration counter
0441     i = i + 1;
0442 
0443     %% compute update step
0444     lambda = struct('eqnonlin', lam(1:neqnln), 'ineqnonlin', mu(1:niqnln));
0445     if nonlinear
0446         if isempty(hess_fcn)
0447             fprintf('mips: Hessian evaluation via finite differences not yet implemented.\n       Please provide your own hessian evaluation function.');
0448         end
0449         Lxx = hess_fcn(x, lambda, opt.cost_mult);
0450     else
0451         [f_, df_, d2f] = f_fcn(x);      %% cost
0452         Lxx = d2f * opt.cost_mult;
0453     end
0454     zinvdiag = sparse(1:niq, 1:niq, 1 ./ z, niq, niq);
0455     mudiag = sparse(1:niq, 1:niq, mu, niq, niq);
0456     dh_zinv = dh * zinvdiag;
0457     M = Lxx + dh_zinv * mudiag * dh';
0458     N = Lx + dh_zinv * (mudiag * h + gamma * e);
0459     dxdlam = mplinsolve([M dg; dg' sparse(neq, neq)], [-N; -g], opt.linsolver, []);
0460 %     AAA = [
0461 %         M  dg;
0462 %         dg'  sparse(neq, neq)
0463 %     ];
0464 %     rc = 1/condest(AAA);
0465 %     if rc < 1e-22
0466 %         fprintf('my RCOND = %g\n', rc);
0467 %         n = size(AAA, 1);
0468 %         AAA = AAA + 1e-3 * speye(n,n);
0469 %     end
0470 %     bbb = [-N; -g];
0471 %     dxdlam = AAA \ bbb;
0472     if any(isnan(dxdlam)) || norm(dxdlam) > max_stepsize
0473         if opt.verbose
0474             fprintf('\nNumerically Failed\n');
0475         end
0476         eflag = -1;
0477         break;
0478     end
0479     dx = dxdlam(1:nx, 1);
0480     dlam = dxdlam(nx+(1:neq), 1);
0481     dz = -h - z - dh' * dx;
0482     dmu = -mu + zinvdiag *(gamma*e - mudiag * dz);
0483 
0484     %% optional step-size control
0485     sc = 0;
0486     if opt.step_control
0487         x1 = x + dx;
0488 
0489         %% evaluate cost, constraints, derivatives at x1
0490         [f1, df1] = f_fcn(x1);          %% cost
0491         f1 = f1 * opt.cost_mult;
0492         df1 = df1 * opt.cost_mult;
0493         if nonlinear
0494             [hn1, gn1, dhn1, dgn1] = gh_fcn(x1); %% nonlinear constraints
0495             h1 = [hn1; Ai * x1 - bi];       %% inequality constraints
0496             g1 = [gn1; Ae * x1 - be];       %% equality constraints
0497             dh1 = [dhn1 Ai'];               %% 1st derivative of inequalities
0498             dg1 = [dgn1 Ae'];               %% 1st derivative of equalities
0499         else
0500             h1 = Ai * x1 - bi;              %% inequality constraints
0501             g1 = Ae * x1 - be;              %% equality constraints
0502             dh1 = dh;                       %% 1st derivative of inequalities
0503             dg1 = dg;                       %% 1st derivative of equalities
0504         end
0505 
0506         %% check tolerance
0507         Lx1 = df1 + dg1 * lam + dh1 * mu;
0508         if isempty(h1)
0509             maxh1 = zeros(1,0);
0510         else
0511             maxh1 = max(h1);
0512         end
0513         feascond1 = max([norm(g1, Inf), maxh1]) / (1 + max([ norm(x1, Inf), norm(z, Inf) ]));
0514         gradcond1 = norm(Lx1, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0515 
0516         if feascond1 > feascond && gradcond1 > gradcond
0517             sc = 1;
0518         end
0519     end
0520     if sc
0521         alpha = 1;
0522         for j = 1:opt.sc.red_it
0523             dx1 = alpha * dx;
0524             x1 = x + dx1;
0525             f1 = f_fcn(x1);                 %% cost
0526             f1 = f1 * opt.cost_mult;
0527             if nonlinear
0528                 [hn1, gn1] = gh_fcn(x1);    %% nonlinear constraints
0529                 h1 = [hn1; Ai * x1 - bi];   %% inequality constraints
0530                 g1 = [gn1; Ae * x1 - be];   %% equality constraints
0531             else
0532                 h1 = Ai * x1 - bi;          %% inequality constraints
0533                 g1 = Ae * x1 - be;          %% equality constraints
0534             end
0535             L1 = f1 + lam' * g1 + mu' * (h1+z) - gamma * sum(log(z));
0536             if opt.verbose > 2
0537                 fprintf('\n   %3d            %10g', -j, norm(dx1));
0538             end
0539             rho = (L1 - L) / (Lx' * dx1 + 0.5 * dx1' * Lxx * dx1);
0540             if rho > rho_min && rho < rho_max
0541                 break;
0542             else
0543                 alpha = alpha / 2;
0544             end
0545         end
0546         dx = alpha * dx;
0547         dz = alpha * dz;
0548         dlam = alpha * dlam;
0549         dmu = alpha * dmu;
0550     end
0551 
0552     %% do the update
0553     k = find(dz < 0);
0554     if isempty(k)
0555         alphap = 1;
0556     else
0557         alphap = min( [xi * min(z(k) ./ -dz(k)) 1] );
0558     end
0559     k = find(dmu < 0);
0560     if isempty(k)
0561         alphad = 1;
0562     else
0563         alphad = min( [xi * min(mu(k) ./ -dmu(k)) 1] );
0564     end
0565     x = x + alphap * dx;
0566     z = z + alphap * dz;
0567     lam = lam + alphad * dlam;
0568     mu  = mu  + alphad * dmu;
0569     if niq > 0
0570         gamma = sigma * (z' * mu) / niq;
0571     end
0572 
0573     %% evaluate cost, constraints, derivatives
0574     [f, df] = f_fcn(x);                 %% cost
0575     f = f * opt.cost_mult;
0576     df = df * opt.cost_mult;
0577     if nonlinear
0578         [hn, gn, dhn, dgn] = gh_fcn(x); %% nonlinear constraints
0579         h = [hn; Ai * x - bi];          %% inequality constraints
0580         g = [gn; Ae * x - be];          %% equality constraints
0581         dh = [dhn Ai'];                 %% 1st derivative of inequalities
0582         dg = [dgn Ae'];                 %% 1st derivative of equalities
0583     else
0584         h = Ai * x - bi;                %% inequality constraints
0585         g = Ae * x - be;                %% equality constraints
0586         %% 1st derivatives are constant, still dh = Ai', dg = Ae'
0587     end
0588 
0589     %% check tolerance
0590     Lx = df + dg * lam + dh * mu;
0591     if isempty(h)
0592         maxh = zeros(1,0);
0593     else
0594         maxh = max(h);
0595     end
0596     feascond = max([norm(g, Inf), maxh]) / (1 + max([ norm(x, Inf), norm(z, Inf) ]));
0597     gradcond = norm(Lx, Inf) / (1 + max([ norm(lam, Inf), norm(mu, Inf) ]));
0598     compcond = (z' * mu) / (1 + norm(x, Inf));
0599     costcond = abs(f - f0) / (1 + abs(f0));
0600     %% save history
0601     hist(i+1) = struct('feascond', feascond, 'gradcond', gradcond, ...
0602         'compcond', compcond, 'costcond', costcond, 'gamma', gamma, ...
0603         'stepsize', norm(dx), 'obj', f/opt.cost_mult, ...
0604         'alphap', alphap, 'alphad', alphad);
0605 
0606     if opt.verbose > 1
0607         fprintf('\n%3d  %12.8g %10.5g %12g %12g %12g %12g', ...
0608             i, f/opt.cost_mult, norm(dx), feascond, gradcond, compcond, costcond);
0609     end
0610     if feascond < opt.feastol && gradcond < opt.gradtol && ...
0611                     compcond < opt.comptol && costcond < opt.costtol
0612         converged = 1;
0613         if opt.verbose
0614             fprintf('\nConverged!\n');
0615         end
0616     else
0617         if any(isnan(x)) || alphap < alpha_min || alphad < alpha_min || ...
0618                 gamma < eps || gamma > 1/eps
0619             if opt.verbose
0620                 fprintf('\nNumerically Failed\n');
0621             end
0622             eflag = -1;
0623             break;
0624         end
0625         f0 = f;
0626         if opt.step_control
0627             L = f + lam' * g + mu' * (h+z) - gamma * sum(log(z));
0628         end
0629     end
0630 end
0631 
0632 if opt.verbose
0633     if ~converged
0634         fprintf('\nDid not converge in %d iterations.\n', i);
0635     end
0636 end
0637 
0638 %%-----  package up results  -----
0639 hist = hist(1:i+1);
0640 if eflag ~= -1
0641     eflag = converged;
0642 end
0643 output = struct('iterations', i, 'hist', hist, 'message', '');
0644 if eflag == 0
0645     output.message = 'Did not converge';
0646 elseif eflag == 1
0647     output.message = 'Converged';
0648 elseif eflag == -1
0649     output.message = 'Numerically failed';
0650 else
0651     output.message = 'Please hang up and dial again';
0652 end
0653 
0654 %% zero out multipliers on non-binding constraints
0655 mu(h < -opt.feastol & mu < mu_threshold) = 0;
0656 
0657 %% un-scale cost and prices
0658 f   = f   / opt.cost_mult;
0659 lam = lam / opt.cost_mult;
0660 mu  = mu  / opt.cost_mult;
0661 
0662 %% re-package multipliers into struct
0663 lam_lin = lam((neqnln+1):neq);              %% lambda for linear constraints
0664 mu_lin  = mu((niqnln+1):niq);               %% mu for linear constraints
0665 kl = find(lam_lin < 0);                     %% lower bound binding
0666 ku = find(lam_lin > 0);                     %% upper bound binding
0667 
0668 mu_l = zeros(nx+nA, 1);
0669 mu_l(ieq(kl)) = -lam_lin(kl);
0670 mu_l(igt) = mu_lin(nlt+(1:ngt));
0671 mu_l(ibx) = mu_lin(nlt+ngt+nbx+(1:nbx));
0672 
0673 mu_u = zeros(nx+nA, 1);
0674 mu_u(ieq(ku)) = lam_lin(ku);
0675 mu_u(ilt) = mu_lin(1:nlt);
0676 mu_u(ibx) = mu_lin(nlt+ngt+(1:nbx));
0677 
0678 fields = { ...
0679     'mu_l', mu_l((nx+1):end), ...
0680     'mu_u', mu_u((nx+1):end), ...
0681     'lower', mu_l(1:nx), ...
0682     'upper', mu_u(1:nx) ...
0683 };
0684 
0685 if niqnln > 0
0686     fields = { ...
0687         'ineqnonlin', mu(1:niqnln), ...
0688         fields{:} ...
0689     };
0690 end
0691 if neqnln > 0
0692     fields = { ...
0693         'eqnonlin', lam(1:neqnln), ...
0694         fields{:} ...
0695     };
0696 end
0697 
0698 lambda = struct(fields{:});
0699 
0700 % lambda = struct( ...
0701 %     'eqnonlin', lam(1:neqnln), ...
0702 %     'ineqnonlin', mu(1:niqnln), ...
0703 %     'mu_l', mu_l((nx+1):end), ...
0704 %     'mu_u', mu_u((nx+1):end), ...
0705 %     'lower', mu_l(1:nx), ...
0706 %     'upper', mu_u(1:nx) );

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