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

nleqs_fd_newton

PURPOSE ^

NLEQS_FD_NEWTON Nonlinear Equation Solver based on Newton's method.

SYNOPSIS ^

function [x, f, eflag, output, J] = nleqs_fd_newton(fcn, x0, opt)

DESCRIPTION ^

NLEQS_FD_NEWTON  Nonlinear Equation Solver based on Newton's method.
   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_FD_NEWTON(FCN, X0, OPT)
   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_FD_NEWTON(PROBLEM)
   A function providing a standardized interface for using Newton's
   method to solve the nonlinear equation f(x) = 0, beginning from a
   starting point x0.

   Inputs:
       FCN : handle to function that evaluates the function f(x) to
           be solved. Calling syntax for this function is:
               f = FCN(x)
       X0 : starting value, x0, of vector x
       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
           max_it (30) - maximum number of iterations for fast-decoupled
               Newton's method
           tol (1e-8) - tolerance on Inf-norm of f(x)
           fd_opt - options struct for Newton's method, with field:
               jac_approx_fcn (required) - handle to function with the
                   following calling syntax:
                       JJ = jac_approx_fcn();
                   where JJ is a cell array of matrices whose number of
                   rows sum to the dimension of f and number of columns
                   sum to the dimension of x. JJ{k} is a matrix
                   approximating the Jacobian of block k of f and x.
               labels (optional) - cell array of same length as JJ returned
                   by jac_approx_fcn(), with short string labels for
                   the decoupled blocks
       PROBLEM : The inputs can alternatively be supplied in a single
           PROBLEM struct with fields corresponding to the input arguments
           described above: fcn, x0, opt

   Outputs (all optional, except X):
       X : solution vector x
       F : final function value, f(x)
       EXITFLAG : exit flag
           1 = converged
           0 or negative values = solver specific failure codes
       OUTPUT : output struct with the following fields:
           alg - algorithm code of solver used ('FD')
           iterations - number of iterations performed
           hist - struct array with trajectories of the following:
                   normf
           message - exit message
       JAC : approximate Jacobian matrix

   Note the calling syntax is almost identical to that of FSOLVE from
   MathWorks' Optimization Toolbox. The function for evaluating the
   nonlinear function is identical.

   Calling syntax options:
       [x, f, exitflag, output, jac] = nleqs_fd_newton(fcn, x0);
       [x, f, exitflag, output, jac] = nleqs_fd_newton(fcn, x0, opt);
       x = nleqs_fd_newton(problem);
               where problem is a struct with fields: fcn, x0, opt
               and all fields except 'fcn' and 'x0' are optional
       x = nleqs_fd_newton(...);
       [x, f] = nleqs_fd_newton(...);
       [x, f, exitflag] = nleqs_fd_newton(...);
       [x, f, exitflag, output] = nleqs_fd_newton(...);
       [x, f, exitflag, output, jac] = nleqs_fd_newton(...);

   Example: (problem from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg)
       function f = f2(x)
       f = [  x(1)^2 +   x(1)*x(2)   - 10;
              x(2)   + 3*x(1)*x(2)^2 - 57  ];

       function JJ = jac_approx_fcn2()
       J = [7 2; 27 37];
       JJ = {J(1,1), J(2,2)};

       problem = struct( ...
           'fcn', @(x)f2(x), ...
           'x0',  [0; 0], ...
           'opt', struct( ...
               'verbose', 2, ...
               'fd_opt', struct( ...
                   'jac_approx_fcn', @jac_approx_fcn2 )));
       [x, f, exitflag, output, jac] = nleqs_fd_newton(problem);

   See also NLEQS_MASTER.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, J] = nleqs_fd_newton(fcn, x0, opt)
0002 %NLEQS_FD_NEWTON  Nonlinear Equation Solver based on Newton's method.
0003 %   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_FD_NEWTON(FCN, X0, OPT)
0004 %   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_FD_NEWTON(PROBLEM)
0005 %   A function providing a standardized interface for using Newton's
0006 %   method to solve the nonlinear equation f(x) = 0, beginning from a
0007 %   starting point x0.
0008 %
0009 %   Inputs:
0010 %       FCN : handle to function that evaluates the function f(x) to
0011 %           be solved. Calling syntax for this function is:
0012 %               f = FCN(x)
0013 %       X0 : starting value, x0, of vector x
0014 %       OPT : optional options structure with the following fields,
0015 %           all of which are also optional (default values shown in
0016 %           parentheses)
0017 %           verbose (0) - controls level of progress output displayed
0018 %               0 = no progress output
0019 %               1 = some progress output
0020 %               2 = verbose progress output
0021 %           max_it (30) - maximum number of iterations for fast-decoupled
0022 %               Newton's method
0023 %           tol (1e-8) - tolerance on Inf-norm of f(x)
0024 %           fd_opt - options struct for Newton's method, with field:
0025 %               jac_approx_fcn (required) - handle to function with the
0026 %                   following calling syntax:
0027 %                       JJ = jac_approx_fcn();
0028 %                   where JJ is a cell array of matrices whose number of
0029 %                   rows sum to the dimension of f and number of columns
0030 %                   sum to the dimension of x. JJ{k} is a matrix
0031 %                   approximating the Jacobian of block k of f and x.
0032 %               labels (optional) - cell array of same length as JJ returned
0033 %                   by jac_approx_fcn(), with short string labels for
0034 %                   the decoupled blocks
0035 %       PROBLEM : The inputs can alternatively be supplied in a single
0036 %           PROBLEM struct with fields corresponding to the input arguments
0037 %           described above: fcn, x0, opt
0038 %
0039 %   Outputs (all optional, except X):
0040 %       X : solution vector x
0041 %       F : final function value, f(x)
0042 %       EXITFLAG : exit flag
0043 %           1 = converged
0044 %           0 or negative values = solver specific failure codes
0045 %       OUTPUT : output struct with the following fields:
0046 %           alg - algorithm code of solver used ('FD')
0047 %           iterations - number of iterations performed
0048 %           hist - struct array with trajectories of the following:
0049 %                   normf
0050 %           message - exit message
0051 %       JAC : approximate Jacobian matrix
0052 %
0053 %   Note the calling syntax is almost identical to that of FSOLVE from
0054 %   MathWorks' Optimization Toolbox. The function for evaluating the
0055 %   nonlinear function is identical.
0056 %
0057 %   Calling syntax options:
0058 %       [x, f, exitflag, output, jac] = nleqs_fd_newton(fcn, x0);
0059 %       [x, f, exitflag, output, jac] = nleqs_fd_newton(fcn, x0, opt);
0060 %       x = nleqs_fd_newton(problem);
0061 %               where problem is a struct with fields: fcn, x0, opt
0062 %               and all fields except 'fcn' and 'x0' are optional
0063 %       x = nleqs_fd_newton(...);
0064 %       [x, f] = nleqs_fd_newton(...);
0065 %       [x, f, exitflag] = nleqs_fd_newton(...);
0066 %       [x, f, exitflag, output] = nleqs_fd_newton(...);
0067 %       [x, f, exitflag, output, jac] = nleqs_fd_newton(...);
0068 %
0069 %   Example: (problem from Christi Patton Luks, https://www.youtube.com/watch?v=pJG4yhtgerg)
0070 %       function f = f2(x)
0071 %       f = [  x(1)^2 +   x(1)*x(2)   - 10;
0072 %              x(2)   + 3*x(1)*x(2)^2 - 57  ];
0073 %
0074 %       function JJ = jac_approx_fcn2()
0075 %       J = [7 2; 27 37];
0076 %       JJ = {J(1,1), J(2,2)};
0077 %
0078 %       problem = struct( ...
0079 %           'fcn', @(x)f2(x), ...
0080 %           'x0',  [0; 0], ...
0081 %           'opt', struct( ...
0082 %               'verbose', 2, ...
0083 %               'fd_opt', struct( ...
0084 %                   'jac_approx_fcn', @jac_approx_fcn2 )));
0085 %       [x, f, exitflag, output, jac] = nleqs_fd_newton(problem);
0086 %
0087 %   See also NLEQS_MASTER.
0088 
0089 %   MP-Opt-Model
0090 %   Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC)
0091 %   by Ray Zimmerman, PSERC Cornell
0092 %
0093 %   This file is part of MP-Opt-Model.
0094 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0095 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0096 
0097 %%----- input argument handling  -----
0098 %% gather inputs
0099 if nargin == 1 && isstruct(fcn) %% problem struct
0100     p = fcn;
0101     fcn = p.fcn;
0102     x0 = p.x0;
0103     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0104 else                            %% individual args
0105     if nargin < 3
0106         opt = [];
0107     end
0108 end
0109 nx = size(x0, 1);           %% number of variables
0110 
0111 %% set default options
0112 opt0 = struct(  'verbose', 0, ...
0113                 'max_it', 30, ...
0114                 'tol', 1e-8 );
0115 if isempty(opt)
0116     opt = opt0;
0117 end
0118 if isfield(opt, 'verbose') && ~isempty(opt.verbose)
0119     verbose = opt.verbose;
0120 else
0121     verbose = opt0.verbose;
0122 end
0123 if isfield(opt, 'max_it') && opt.max_it     %% not empty or zero
0124     max_it = opt.max_it;
0125 else
0126     max_it = opt0.max_it;
0127 end
0128 if isfield(opt, 'tol') && opt.tol           %% not empty or zero
0129     tol = opt.tol;
0130 else
0131     tol = opt0.tol;
0132 end
0133 if isfield(opt, 'fd_opt') && isfield(opt.fd_opt, 'jac_approx_fcn')
0134     jac_approx_fcn = opt.fd_opt.jac_approx_fcn;
0135     if isfield(opt.fd_opt, 'labels')
0136         labels = opt.fd_opt.labels;
0137     else
0138         labels = {};
0139     end
0140 else
0141     error('nleqs_fd_newton: required ''fd_opt.jac_approx_fcn'' option missing');
0142 end
0143 
0144 %% initialize
0145 lu_vec = have_feature('lu_vec');
0146 eflag = 0;
0147 i = 0;
0148 x = x0;
0149 hist(max_it+1) = struct('normf', 0);
0150 
0151 %% get Jacobian approximation matrices
0152 JJ = jac_approx_fcn();
0153 
0154 %% get indices for each block
0155 nj = length(JJ);
0156 normf = zeros(nj, 1);   %% initialize
0157 b = 0;
0158 i1 = zeros(1, nj); iN = i1;
0159 for j = 1:nj
0160     [m, n] = size(JJ{j});
0161     if m ~= n
0162         error('nleqs_fd_newton: Jacobian approximation for block %d is not square (%d x %d)', j, m, n);
0163     end
0164     i1(j) = b+1;
0165     iN(j) = b+m;
0166     b = iN(j);
0167 end
0168 
0169 %% block labels
0170 if isempty(labels)  %% create default labels 'A', 'B', ...
0171     for j = 1:nj
0172         labels{j} = sprintf('%s', char('A'+j-1));
0173     end
0174 elseif length(labels) ~= nj
0175     error('nleqs_fd_newton: size of labels must be consistent with jac_approx_fcn');
0176 end
0177 
0178 %% evaluate f(x0)
0179 f = fcn(x);
0180 
0181 %% check tolerance
0182 for j = 1:nj
0183     normf(j) = norm(f(i1(j):iN(j)), inf);
0184 end
0185 if verbose > 1
0186     fprintf('\n iteration ');
0187     for j = 1:nj, fprintf('   max residual ');          end
0188     fprintf('\nblock    # ');
0189     for j = 1:nj
0190         n1 = length(labels{j});
0191         n2 = fix((14 - 3 - n1) / 2)+2;
0192         n3 = 13-n1-n2;
0193         lb = sprintf('%sf[%s]%s', repmat(' ', 1, n2), labels{j}, repmat(' ', 1, n3));
0194         fprintf('%16s', lb);
0195     end
0196     fprintf('\n------ ----');
0197     for j = 1:nj, fprintf('  --------------');          end
0198     fprintf('\n  -    %3d', i);
0199     for j = 1:nj, fprintf('      %10.3e', normf(j));    end
0200 end
0201 if max(normf) < tol
0202     eflag = 1;
0203     msg = sprintf('Fast-decoupled Newton''s method converged in %d iterations.', i);
0204     if verbose > 1
0205         fprintf('\nConverged!\n');
0206     end
0207 end
0208 
0209 %% factor Jacobian approximation matrices
0210 if lu_vec
0211     for j = nj:-1:1     %% backwards to init cell arrays to full-size at start
0212         if size(JJ{j}, 1) > 1
0213             [L{j}, U{j}, p{j}, q] = lu(JJ{j}, 'vector');
0214             [junk, iq{j}] = sort(q);
0215         else
0216             [L{j}, U{j}, p{j}] = lu(JJ{j}, 'vector');
0217             iq{j} = 1;
0218         end
0219     end
0220 else
0221     for j = nj:-1:1     %% backwards to init cell arrays to full-size at start
0222         [L{j}, U{j}, P{j}] = lu(J{j});
0223     end
0224 end
0225 
0226 %% save history
0227 hist(i+1).normf = normf;
0228 
0229 %% do Newton iterations
0230 while (~eflag && i < max_it)
0231     %% update iteration counter
0232     i = i + 1;
0233 
0234     %% do decoupled iterations
0235     for j = 1:nj
0236         ff = f(i1(j):iN(j));
0237         if lu_vec
0238             dx = -( U{j} \  (L{j} \ ff(p{j})) );
0239             dx = dx(iq{j});
0240         else
0241             dx = -( U{j} \  (L{j} \ (P{j} * ff)) );
0242         end
0243 
0244         %% update x
0245         x(i1(j):iN(j)) = x(i1(j):iN(j)) + dx;
0246 
0247         %% evalute f(x) and J(x)
0248         f = fcn(x);
0249 
0250         %% check for convergence
0251         for jj = 1:nj
0252             normf(jj) = norm(f(i1(jj):iN(jj)), inf);
0253         end
0254         if verbose > 1
0255             n1 = length(labels{j});
0256             n2 = fix((6 - n1)/2);
0257             n3 = 6-n1-n2;
0258             lb = sprintf('%s%s%s', repmat(' ', 1, n2), labels{j}, repmat(' ', 1, n3));
0259             fprintf('\n%s %3d', lb, i);
0260             for jj = 1:nj
0261                 fprintf('      %10.3e', normf(jj));
0262             end
0263         end
0264         if max(normf) < tol
0265             eflag = 1;
0266             msg = sprintf('Fast-decoupled Newton''s method converged in ');
0267             for jj = 1:nj
0268                 msg = sprintf('%s%d %s-', msg, i - (j<jj), labels{jj});
0269                 if jj ~= nj
0270                     msg = sprintf('%s and ', msg);
0271                 end
0272             end
0273             msg = sprintf('%siterations.', msg);
0274             break;
0275         end
0276     end
0277 
0278     %% save history
0279     hist(i+1).normf = normf;
0280 end
0281 if eflag ~= 1
0282     msg = sprintf('Fast-decoupled Newton''s method did not converge in %d iterations.', i);
0283 end
0284 if verbose
0285     fprintf('\n%s\n', msg);
0286 end
0287 if nargout > 3
0288     output = struct('alg', 'FD', ...
0289                     'iterations', i, ...
0290                     'hist', hist(1:i+1), ...
0291                     'message', msg  );
0292     if nargout > 4
0293         nx = length(x);
0294         if issparse(JJ{1})
0295             J = sparse(nx, nx);
0296         else
0297             J = zeros(nx, nx);
0298         end
0299         for j = 1:nj
0300             J(i1(j):iN(j), i1(j):iN(j)) = JJ{j};
0301         end
0302     end
0303 end

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