diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 55a52d85fab571ef4b3998feef08c7ad3d8bc0a7..bf3c641cfb96b0808118d3f54c25950a1b4594c2 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -3015,6 +3015,19 @@ Finding the steady state with Dynare nonlinear solver
efficient because only the relevant elements are recomputed at
every iteration.
+ ``15``
+
+ Splits the model into recursive blocks and solves each
+ block in turn using the legacy implementation of the trust-region
+ solver with autoscaling. May occasionally work where the new
+ implementation fails and vice versa.
+
+ ``16``
+
+ Legacy implementation of the trust-region algorithm applied to the entire
+ model (same as value ``15``, but applied to the entire model, without splitting).
+ May occasionally work where the new implementation fails and vice versa.
+
|br| Default value is ``4``.
.. option:: homotopy_mode = INTEGER
diff --git a/license.txt b/license.txt
index e57a5163760ad0087309765fb3bf45e911ec9d44..570c1f3cdc10a9e401c10a1fe1d3bf8f7c6839b4 100644
--- a/license.txt
+++ b/license.txt
@@ -100,6 +100,11 @@ Copyright: 2001-2012 Nikolaus Hansen
2012-2017 Dynare Team
License: GPL-3+
+Files: matlab/optimization/trust_region_legacy.m
+Copyright: 2008-2012 VZLU Prague, a.s.
+ 2014-2023 Dynare Team
+License: GPL-3+
+
Files: matlab/optimization/solvopt.m
Copyright: 1997-2008 Alexei Kuntsevich and Franz Kappel
2008-2015 Giovanni Lombardo
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 292ef918c5294f20a65d9ef49eb65750b12ab043..fb629f4013c8c779255171600232832a89f85d33 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -39,8 +39,8 @@ function [ys,params,info] = evaluate_steady_state(ys_init,exo_ss,M_,options_,ste
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-if options_.solve_algo < 0 || options_.solve_algo > 14
- error('STEADY: solve_algo must be between 0 and 14')
+if options_.solve_algo < 0 || options_.solve_algo > 16
+ error('STEADY: solve_algo must be between 0 and 16')
end
if ~options_.bytecode && ~options_.block && options_.solve_algo > 4 && ...
diff --git a/matlab/optimization/dynare_solve.m b/matlab/optimization/dynare_solve.m
index b41f82e73ed4820d8031ccd09ce34ee428410c4e..7b6eab1d552a073e3b3ca759e95fa362f5590728 100644
--- a/matlab/optimization/dynare_solve.m
+++ b/matlab/optimization/dynare_solve.m
@@ -223,11 +223,13 @@ elseif ismember(options.solve_algo, [1, 12])
elseif options.solve_algo==9
[x, errorflag, errorcode] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, varargin{:});
[fvec, fjac] = feval(f, x, varargin{:});
-elseif ismember(options.solve_algo, [2, 4])
+elseif ismember(options.solve_algo, [2, 4, 15])
if options.solve_algo == 2
solver = @solve1;
- else
+ elseif options.solve_algo == 4
solver = @trust_region;
+ else
+ solver = @trust_region_legacy;
end
if ~jacobian_flag
fjac = zeros(nn,nn) ;
@@ -336,6 +338,9 @@ elseif ismember(options.solve_algo, [13, 14])
options.solve_algo == 13, ... % Only block-decompose with Dulmage-Mendelsohn for 13, not for 14
options.debug, varargin{:});
[fvec, fjac] = feval(f, x, varargin{:});
+elseif options.solve_algo==16
+ [x, errorflag, errorcode] = trust_region_legacy(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, varargin{:});
+ [fvec, fjac] = feval(f, x, varargin{:});
else
- error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
+ error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14,15,16]')
end
diff --git a/matlab/optimization/trust_region_legacy.m b/matlab/optimization/trust_region_legacy.m
new file mode 100644
index 0000000000000000000000000000000000000000..7bc646eb67a05330e9d9a2e104fe575c481d437c
--- /dev/null
+++ b/matlab/optimization/trust_region_legacy.m
@@ -0,0 +1,226 @@
+function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,factor,debug,varargin)
+% Solves systems of non linear equations of several variables, using a
+% trust-region method.
+%
+% INPUTS
+% fcn: name of the function to be solved
+% x0: guess values
+% j1: equations index for which the model is solved
+% j2: unknown variables index
+% jacobian_flag=true: jacobian given by the 'func' function
+% jacobian_flag=false: jacobian obtained numerically
+% gstep increment multiplier in numercial derivative
+% computation
+% tolf tolerance for residuals
+% tolx tolerance for solution variation
+% maxiter maximum number of iterations
+% factor real scalar, determines the initial step bound
+% debug debug flag
+% varargin: list of arguments following bad_cond_flag
+%
+% OUTPUTS
+% x: results
+% check=1: the model can not be solved
+% info: detailed exitcode
+% SPECIAL REQUIREMENTS
+% none
+
+% Copyright (C) 2008-2012 VZLU Prague, a.s.
+% Copyright (C) 2014-2021 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+%
+% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave
+
+if (ischar (fcn))
+ fcn = str2func (fcn);
+end
+
+n = length(j1);
+
+% These defaults are rather stringent. I think that normally, user
+% prefers accuracy to performance.
+
+macheps = eps (class (x0));
+
+niter = 1;
+
+x = x0;
+info = 0;
+
+% Initial evaluation.
+% Handle arbitrary shapes of x and f and remember them.
+fvec = fcn (x, varargin{:});
+fvec = fvec(j1);
+fn = norm (fvec);
+recompute_jacobian = true;
+
+% Outer loop.
+while (niter < maxiter && ~info)
+
+ % Calculate Jacobian (possibly via FD).
+ if recompute_jacobian
+ if jacobian_flag
+ [~, fjac] = fcn (x, varargin{:});
+ fjac = fjac(j1,j2);
+ else
+ dh = max(abs(x(j2)),gstep(1)*ones(n,1))*eps^(1/3);
+
+ for j = 1:n
+ xdh = x ;
+ xdh(j2(j)) = xdh(j2(j))+dh(j) ;
+ t = fcn(xdh,varargin{:});
+ fjac(:,j) = (t(j1) - fvec)./dh(j) ;
+ end
+ end
+ recompute_jacobian = false;
+ end
+
+ % Get column norms, use them as scaling factors.
+ jcn = sqrt(sum(fjac.*fjac))';
+ if (niter == 1)
+ dg = jcn;
+ dg(dg == 0) = 1;
+ else
+ % Rescale adaptively.
+ dg = max (dg, jcn);
+ end
+
+ if (niter == 1)
+ xn = norm (dg .* x(j2));
+ % FIXME: something better?
+ delta = max (xn, 1)*factor;
+ end
+
+ % Get trust-region model (dogleg) minimizer.
+ s = - dogleg (fjac, fvec, dg, delta);
+ w = fvec + fjac * s;
+
+ sn = norm (dg .* s);
+ if (niter == 1)
+ delta = min (delta, sn);
+ end
+
+ x2 = x;
+ x2(j2) = x2(j2) + s;
+ fvec1 = fcn (x2, varargin{:});
+ fvec1 = fvec1(j1);
+ fn1 = norm (fvec1);
+
+ if (fn1 < fn)
+ % Scaled actual reduction.
+ actred = 1 - (fn1/fn)^2;
+ else
+ actred = -1;
+ end
+
+ % Scaled predicted reduction, and ratio.
+ t = norm (w);
+ if (t < fn)
+ prered = 1 - (t/fn)^2;
+ ratio = actred / prered;
+ else
+ prered = 0;
+ ratio = 0;
+ end
+
+ % Update delta.
+ if (ratio < 0.1)
+ delta = 0.5*delta;
+ if (delta <= 1e1*macheps*xn)
+ % Trust region became uselessly small.
+ if (fn1 <= tolf)
+ info = 1;
+ else
+ info = -3;
+ end
+ break
+ end
+ elseif (abs (1-ratio) <= 0.1)
+ delta = 1.4142*sn;
+ elseif (ratio >= 0.5)
+ delta = max (delta, 1.4142*sn);
+ end
+
+ if (ratio >= 1e-4)
+ % Successful iteration.
+ x(j2) = x(j2) + s;
+ xn = norm (dg .* x(j2));
+ fvec = fvec1;
+ fn = fn1;
+ recompute_jacobian = true;
+ end
+
+ niter = niter + 1;
+
+
+ % Tests for termination condition
+ if (fn <= tolf)
+ info = 1;
+ end
+end
+if info==1
+ check = 0;
+else
+ check = 1;
+end
+end
+
+
+% Solve the double dogleg trust-region least-squares problem:
+% Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
+% x being a convex combination of the gauss-newton and scaled gradient.
+
+% TODO: error checks
+% TODO: handle singularity, or leave it up to mldivide?
+
+function x = dogleg (r, b, d, delta)
+% Get Gauss-Newton direction.
+if isoctave || matlab_ver_less_than('9.3')
+ % The decomposition() function does not exist in Octave and MATLAB < R2017b
+ x = r \ b;
+else
+ x = decomposition(r, 'CheckCondition', false) \ b;
+end
+xn = norm (d .* x);
+if (xn > delta)
+ % GN is too big, get scaled gradient.
+ s = (r' * b) ./ d;
+ sn = norm (s);
+ if (sn > 0)
+ % Normalize and rescale.
+ s = (s / sn) ./ d;
+ % Get the line minimizer in s direction.
+ tn = norm (r*s);
+ snm = (sn / tn) / tn;
+ if (snm < delta)
+ % Get the dogleg path minimizer.
+ bn = norm (b);
+ dxn = delta/xn; snmd = snm/delta;
+ t = (bn/sn) * (bn/xn) * snmd;
+ t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
+ alpha = dxn*(1-snmd^2) / t;
+ else
+ alpha = 0;
+ end
+ else
+ alpha = delta / xn;
+ snm = 0;
+ end
+ % Form the appropriate convex combination.
+ x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
+end
+end
\ No newline at end of file
diff --git a/meson.build b/meson.build
index 5e6cce3ae8e42ff37a91da735207e505b5c5f601..d4d1731af007b3b9cb5fa48fc996e93074c6c52e 100644
--- a/meson.build
+++ b/meson.build
@@ -1059,6 +1059,8 @@ mod_and_m_tests = [
{ 'test' : [ 'steady_state/multi_leads.mod' ] },
{ 'test' : [ 'steady_state/example1_trust_region.mod' ] },
{ 'test' : [ 'steady_state/example1_block_trust_region.mod' ] },
+ { 'test' : [ 'steady_state/example1_trust_region_legacy.mod' ] },
+ { 'test' : [ 'steady_state/example1_block_trust_region_legacy.mod' ] },
{ 'test' : [ 'steady_state/Gali_2015_chapter_6_4.mod' ] },
{ 'test' : [ 'steady_state/ramst_solve_algo_0.mod' ] },
{ 'test' : [ 'steady_state/fs2000_ssfile.mod' ],
diff --git a/preprocessor b/preprocessor
index e32025a76fb156a9af8101f9f43faf7eaa2f72ef..5fa91a08f68d04a8baac4f32dab9db62b5a3546f 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit e32025a76fb156a9af8101f9f43faf7eaa2f72ef
+Subproject commit 5fa91a08f68d04a8baac4f32dab9db62b5a3546f
diff --git a/tests/steady_state/example1_block_trust_region_legacy.mod b/tests/steady_state/example1_block_trust_region_legacy.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d8b2d1a688751463f6271e18cb4a426958d4c849
--- /dev/null
+++ b/tests/steady_state/example1_block_trust_region_legacy.mod
@@ -0,0 +1,47 @@
+// Test block trust region nonlinear solver (solve_algo=13)
+var y, c, k, a, h, b;
+varexo e, u;
+
+parameters beta, rho, alpha, delta, theta, psi, tau;
+
+alpha = 0.36;
+rho = 0.95;
+tau = 0.025;
+beta = 0.99;
+delta = 0.025;
+psi = 0;
+theta = 2.95;
+
+phi = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+ *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1);
+a = rho*a(-1)+tau*b(-1) + e;
+b = tau*a(-1)+rho*b(-1) + u;
+end;
+
+initval;
+y = 1;
+c = 0.8;
+h = 0.3;
+k = 10;
+a = 0;
+b = 0;
+e = 0;
+u = 0;
+end;
+
+options_.debug = true;
+steady(solve_algo=15);
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+end;
+
+stoch_simul(order=1,nomoments,irf=0);
diff --git a/tests/steady_state/example1_trust_region_legacy.mod b/tests/steady_state/example1_trust_region_legacy.mod
new file mode 100644
index 0000000000000000000000000000000000000000..c4eda2ad2d4d45f2f810a80ccabd183c76a1520f
--- /dev/null
+++ b/tests/steady_state/example1_trust_region_legacy.mod
@@ -0,0 +1,46 @@
+// Test trust region nonlinear solver (solve_algo=4)
+var y, c, k, a, h, b;
+varexo e, u;
+
+parameters beta, rho, alpha, delta, theta, psi, tau;
+
+alpha = 0.36;
+rho = 0.95;
+tau = 0.025;
+beta = 0.99;
+delta = 0.025;
+psi = 0;
+theta = 2.95;
+
+phi = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+ *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1);
+a = rho*a(-1)+tau*b(-1) + e;
+b = tau*a(-1)+rho*b(-1) + u;
+end;
+
+initval;
+y = 1;
+c = 0.8;
+h = 0.3;
+k = 10;
+a = 0;
+b = 0;
+e = 0;
+u = 0;
+end;
+
+steady(solve_algo=16);
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+end;
+
+stoch_simul(order=1,nomoments,irf=0);