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.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