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

nleqs_core

PURPOSE ^

NLEQS_CORE Core Nonlinear Equation Solver with arbitrary update function.

SYNOPSIS ^

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

DESCRIPTION ^

NLEQS_CORE  Core Nonlinear Equation Solver with arbitrary update function.
   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_CORE(SP, FCN, X0, OPT)
   A function providing the core code for a standardized interface for
   various methods of solving the nonlinear equation f(x) = 0, beginning
   from a starting point x0. Allows for an arbitrary update function.

   Inputs:
       SP : solver parameters, struct with the following fields (all required)
           alg : algorithm code, e.g. 'NEWTON', 'GS', etc.
           name : user-visible name to be followed by 'method' in output
               'Newton''s', 'Gauss-Seidel', etc.
           default_max_it : default value for max_it option, if not provided
           need_jac : 1 - compute Jacobian, 0 - do not compute Jacobian
           update_fcn : handle to function used to update x, with syntax:
               x = update_fcn(x, f)
               x = update_fcn(x, f, J)
       FCN : handle to function that evaluates the function f(x) to
           be solved and (optionally) its Jacobian, J(x). Calling syntax
           for this function is:
               f = FCN(x)
               [f, J] = FCN(x)
           If f and x are n x 1, then J is the n x n matrix of partial
           derivatives of f (rows) w.r.t. x (cols).
       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 (SP.default_max_it) - maximum number of iterations
           tol (1e-8) - tolerance on Inf-norm of f(x)

   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 (SP.alg)
           iterations - number of iterations performed
           hist - struct array with trajectories of the following:
                   normf
                   normdx
           message - exit message
       JAC : final Jacobian matrix, J

   Calling syntax options:
       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0);
       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);
       x = nleqs_core(...);
       [x, f] = nleqs_core(...);
       [x, f, exitflag] = nleqs_core(...);
       [x, f, exitflag, output] = nleqs_core(...);
       [x, f, exitflag, output, jac] = nleqs_core(...);

   Example: (problem from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/)
       function [f, J] = f1(x)
       f = [  x(1)   + x(2) - 1;
             -x(1)^2 + x(2) + 5    ];
       if nargout > 1
           J = [1 1; -2*x(1) 1];
       end

       fcn = @(x)f1(x);
       x0  = [0; 0];
       opt = struct('verbose', 2);
       sp = struct( ...
           'alg',              'NEWTON', ...
           'name',             'Newton''s', ...
           'default_max_it',   10, ...
           'need_jac',         1, ...
           'update_fcn',       @newton_update_fcn, ...
           )
       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);

   See also NLEQS_MASTER, NLEQS_NEWTON, NLEQS_GAUSS_SEIDEL.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, J] = nleqs_core(sp, fcn, x0, opt)
0002 %NLEQS_CORE  Core Nonlinear Equation Solver with arbitrary update function.
0003 %   [X, F, EXITFLAG, OUTPUT, JAC] = NLEQS_CORE(SP, FCN, X0, OPT)
0004 %   A function providing the core code for a standardized interface for
0005 %   various methods of solving the nonlinear equation f(x) = 0, beginning
0006 %   from a starting point x0. Allows for an arbitrary update function.
0007 %
0008 %   Inputs:
0009 %       SP : solver parameters, struct with the following fields (all required)
0010 %           alg : algorithm code, e.g. 'NEWTON', 'GS', etc.
0011 %           name : user-visible name to be followed by 'method' in output
0012 %               'Newton''s', 'Gauss-Seidel', etc.
0013 %           default_max_it : default value for max_it option, if not provided
0014 %           need_jac : 1 - compute Jacobian, 0 - do not compute Jacobian
0015 %           update_fcn : handle to function used to update x, with syntax:
0016 %               x = update_fcn(x, f)
0017 %               x = update_fcn(x, f, J)
0018 %       FCN : handle to function that evaluates the function f(x) to
0019 %           be solved and (optionally) its Jacobian, J(x). Calling syntax
0020 %           for this function is:
0021 %               f = FCN(x)
0022 %               [f, J] = FCN(x)
0023 %           If f and x are n x 1, then J is the n x n matrix of partial
0024 %           derivatives of f (rows) w.r.t. x (cols).
0025 %       X0 : starting value, x0, of vector x
0026 %       OPT : optional options structure with the following fields,
0027 %           all of which are also optional (default values shown in
0028 %           parentheses)
0029 %           verbose (0) - controls level of progress output displayed
0030 %               0 = no progress output
0031 %               1 = some progress output
0032 %               2 = verbose progress output
0033 %           max_it (SP.default_max_it) - maximum number of iterations
0034 %           tol (1e-8) - tolerance on Inf-norm of f(x)
0035 %
0036 %   Outputs (all optional, except X):
0037 %       X : solution vector x
0038 %       F : final function value, f(x)
0039 %       EXITFLAG : exit flag
0040 %           1 = converged
0041 %           0 or negative values = solver specific failure codes
0042 %       OUTPUT : output struct with the following fields:
0043 %           alg - algorithm code of solver used (SP.alg)
0044 %           iterations - number of iterations performed
0045 %           hist - struct array with trajectories of the following:
0046 %                   normf
0047 %                   normdx
0048 %           message - exit message
0049 %       JAC : final Jacobian matrix, J
0050 %
0051 %   Calling syntax options:
0052 %       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0);
0053 %       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);
0054 %       x = nleqs_core(...);
0055 %       [x, f] = nleqs_core(...);
0056 %       [x, f, exitflag] = nleqs_core(...);
0057 %       [x, f, exitflag, output] = nleqs_core(...);
0058 %       [x, f, exitflag, output, jac] = nleqs_core(...);
0059 %
0060 %   Example: (problem from https://www.chilimath.com/lessons/advanced-algebra/systems-non-linear-equations/)
0061 %       function [f, J] = f1(x)
0062 %       f = [  x(1)   + x(2) - 1;
0063 %             -x(1)^2 + x(2) + 5    ];
0064 %       if nargout > 1
0065 %           J = [1 1; -2*x(1) 1];
0066 %       end
0067 %
0068 %       fcn = @(x)f1(x);
0069 %       x0  = [0; 0];
0070 %       opt = struct('verbose', 2);
0071 %       sp = struct( ...
0072 %           'alg',              'NEWTON', ...
0073 %           'name',             'Newton''s', ...
0074 %           'default_max_it',   10, ...
0075 %           'need_jac',         1, ...
0076 %           'update_fcn',       @newton_update_fcn, ...
0077 %           )
0078 %       [x, f, exitflag, output, jac] = nleqs_core(sp, fcn, x0, opt);
0079 %
0080 %   See also NLEQS_MASTER, NLEQS_NEWTON, NLEQS_GAUSS_SEIDEL.
0081 
0082 %   MP-Opt-Model
0083 %   Copyright (c) 1996-2020, Power Systems Engineering Research Center (PSERC)
0084 %   by Ray Zimmerman, PSERC Cornell
0085 %
0086 %   This file is part of MP-Opt-Model.
0087 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0088 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0089 
0090 %% set default options
0091 opt0 = struct(  'verbose', 0, ...
0092                 'max_it', sp.default_max_it, ...
0093                 'tol', 1e-8 );
0094 if isempty(opt)
0095     opt = opt0;
0096 end
0097 if isfield(opt, 'verbose') && ~isempty(opt.verbose)
0098     verbose = opt.verbose;
0099 else
0100     verbose = opt0.verbose;
0101 end
0102 if isfield(opt, 'max_it') && opt.max_it     %% not empty or zero
0103     max_it = opt.max_it;
0104 else
0105     max_it = opt0.max_it;
0106 end
0107 if isfield(opt, 'tol') && opt.tol           %% not empty or zero
0108     tol = opt.tol;
0109 else
0110     tol = opt0.tol;
0111 end
0112 
0113 %% initialize
0114 eflag = 0;
0115 i = 0;
0116 x = x0;
0117 hist(max_it+1) = struct('normf', 0, 'normdx', 0);
0118 
0119 %% evaluate f(x0)
0120 if sp.need_jac
0121     [f, J] = fcn(x);
0122 else
0123     f = fcn(x);
0124 end
0125 
0126 %% check tolerance
0127 normf = norm(f, inf);
0128 if verbose > 1
0129     fprintf('\n it    max residual        max ∆x');
0130     fprintf('\n----  --------------  --------------');
0131     fprintf('\n%3d     %10.3e           -    ', i, normf);
0132 end
0133 if normf < tol
0134     eflag = 1;
0135     msg = sprintf('%s method converged in %d iterations.', sp.name, i);
0136     if verbose > 1
0137         fprintf('\nConverged!\n');
0138     end
0139 end
0140 
0141 %% save history
0142 hist(i+1).normf = normf;
0143 
0144 %% do solver iterations
0145 while (~eflag && i < max_it)
0146     %% update iteration counter
0147     i = i + 1;
0148 
0149     %% save previous x
0150     xp = x;
0151 
0152     if sp.need_jac
0153         %% update x
0154         x = sp.update_fcn(x, f, J);
0155 
0156         %% evalute f(x) and J(x)
0157         [f, J] = fcn(x);
0158     else
0159         %% update x
0160         x = sp.update_fcn(x, f);
0161 
0162         %% evaluate f(x)
0163         f = fcn(x);
0164     end
0165 
0166     %% check for convergence
0167     normf  = norm(f, inf);
0168     normdx = norm(x-xp, inf);
0169     if verbose > 1
0170         fprintf('\n%3d     %10.3e      %10.3e', i, normf, normdx);
0171     end
0172     if normf < tol
0173         eflag = 1;
0174         msg = sprintf('%s method converged in %d iterations.', sp.name, i);
0175     end
0176 
0177     %% save history
0178     hist(i+1).normf  = normf;
0179     hist(i+1).normfx = normdx;
0180 end
0181 if eflag ~= 1
0182     msg = sprintf('%s method did not converge in %d iterations.', sp.name, i);
0183 end
0184 if verbose
0185     fprintf('\n%s\n', msg);
0186 end
0187 if nargout > 3
0188     output = struct('alg', sp.alg, ...
0189                     'iterations', i, ...
0190                     'hist', hist(1:i+1), ...
0191                     'message', msg  );
0192     if nargout > 4 && ~sp.need_jac
0193         J = [];
0194     end
0195 end

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