1function [x,info,fvec,fjac] = dynare_solve(func,x,options,varargin) 2% function [x,info,fvec,fjac] = dynare_solve(func,x,options,varargin) 3% proposes different solvers 4% 5% INPUTS 6% func: name of the function to be solved 7% x: guess values 8% options: struct of Dynare options 9% varargin: list of arguments following jacobian_flag 10% 11% OUTPUTS 12% x: solution 13% info=1: the model can not be solved 14% fvec: Function value (used for debugging when check=1) 15% fjac: Jacobian (used for debugging when check=1) 16% 17% SPECIAL REQUIREMENTS 18% none 19 20% Copyright (C) 2001-2017 Dynare Team 21% 22% This file is part of Dynare. 23% 24% Dynare is free software: you can redistribute it and/or modify 25% it under the terms of the GNU General Public License as published by 26% the Free Software Foundation, either version 3 of the License, or 27% (at your option) any later version. 28% 29% Dynare is distributed in the hope that it will be useful, 30% but WITHOUT ANY WARRANTY; without even the implied warranty of 31% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 32% GNU General Public License for more details. 33% 34% You should have received a copy of the GNU General Public License 35% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 36 37% jacobian_flag=true: jacobian given by the 'func' function 38% jacobian_flag=false: jacobian obtained numerically 39jacobian_flag = options.jacobian_flag; 40 41% Set tolerance parameter depending the the caller function. 42stack = dbstack; 43if isoctave 44 [pathstr,name,ext]=fileparts(stack(2).file); 45 caller_file_name=[name,ext]; 46else 47 caller_file_name=stack(2).file; 48end 49if strcmp(caller_file_name,'solve_stacked_problem.m') 50 tolf = options.dynatol.f; 51 tolx = options.dynatol.x; 52else 53 tolf = options.solve_tolf; 54 tolx = options.solve_tolx; 55end 56 57if strcmp(caller_file_name,'dyn_ramsey_static.m') 58 maxit = options.ramsey.maxit; 59else 60 maxit = options.steady.maxit; 61end 62 63info = 0; 64nn = size(x,1); 65 66% Get status of the initial guess (default values?) 67if any(x) 68 % The current initial guess is not the default for all the variables. 69 idx = find(x); % Indices of the variables with default initial guess values. 70 in0 = length(idx); 71else 72 % The current initial guess is the default for all the variables. 73 idx = transpose(1:nn); 74 in0 = nn; 75end 76 77% checking initial values 78if jacobian_flag 79 [fvec, fjac] = feval(func, x, varargin{:}); 80 wrong_initial_guess_flag = false; 81 if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) ... 82 || any(~isreal(fvec)) || any(~isreal(fjac(:))) 83 if max(abs(fvec)) < tolf %return if initial value solves problem 84 info = 0; 85 return; 86 end 87 disp('Randomize initial guess...') 88 % Let's try random numbers for the variables initialized with the default value. 89 wrong_initial_guess_flag = true; 90 % First try with positive numbers. 91 tentative_number = 0; 92 while wrong_initial_guess_flag && tentative_number<=in0*10 93 tentative_number = tentative_number+1; 94 x(idx) = rand(in0, 1)*10; 95 [fvec, fjac] = feval(func, x, varargin{:}); 96 wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))); 97 end 98 % If all previous attempts failed, try with real numbers. 99 tentative_number = 0; 100 while wrong_initial_guess_flag && tentative_number<=in0*10 101 tentative_number = tentative_number+1; 102 x(idx) = randn(in0, 1)*10; 103 [fvec, fjac] = feval(func, x, varargin{:}); 104 wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))); 105 end 106 % Last tentative, ff all previous attempts failed, try with negative numbers. 107 tentative_number = 0; 108 while wrong_initial_guess_flag && tentative_number<=in0*10 109 tentative_number = tentative_number+1; 110 x(idx) = -rand(in0, 1)*10; 111 [fvec, fjac] = feval(func, x, varargin{:}); 112 wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))); 113 end 114 end 115else 116 fvec = feval(func,x,varargin{:}); 117 fjac = zeros(nn,nn); 118 wrong_initial_guess_flag = false; 119 if ~all(isfinite(fvec)) 120 % Let's try random numbers for the variables initialized with the default value. 121 wrong_initial_guess_flag = true; 122 % First try with positive numbers. 123 tentative_number = 0; 124 while wrong_initial_guess_flag && tentative_number<=in0*10 125 tentative_number = tentative_number+1; 126 x(idx) = rand(in0, 1)*10; 127 fvec = feval(func, x, varargin{:}); 128 wrong_initial_guess_flag = ~all(isfinite(fvec)); 129 end 130 % If all previous attempts failed, try with real numbers. 131 tentative_number = 0; 132 while wrong_initial_guess_flag && tentative_number<=in0*10 133 tentative_number = tentative_number+1; 134 x(idx) = randn(in0, 1)*10; 135 fvec = feval(func, x, varargin{:}); 136 wrong_initial_guess_flag = ~all(isfinite(fvec)); 137 end 138 % Last tentative, ff all previous attempts failed, try with negative numbers. 139 tentative_number = 0; 140 while wrong_initial_guess_flag && tentative_number<=in0*10 141 tentative_number = tentative_number+1; 142 x(idx) = -rand(in0, 1)*10; 143 fvec = feval(func, x, varargin{:}); 144 wrong_initial_guess_flag = ~all(isfinite(fvec)); 145 end 146 end 147end 148 149% Exit with error if no initial guess has been found. 150if wrong_initial_guess_flag 151 info=1; 152 x = NaN(size(fvec)); 153 return 154end 155 156% this test doesn't check complementarity conditions and is not used for 157% mixed complementarity problems 158if (~ismember(options.solve_algo,[10,11])) && (max(abs(fvec)) < tolf) 159 return ; 160end 161 162if options.solve_algo == 0 163 if ~isoctave 164 if ~user_has_matlab_license('optimization_toolbox') 165 error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox') 166 end 167 end 168 options4fsolve=optimset('fsolve'); 169 options4fsolve.MaxFunEvals = 50000; 170 options4fsolve.MaxIter = maxit; 171 options4fsolve.TolFun = tolf; 172 if options.debug==1 173 options4fsolve.Display = 'final'; 174 else 175 options4fsolve.Display = 'off'; 176 end 177 if jacobian_flag 178 options4fsolve.Jacobian = 'on'; 179 else 180 options4fsolve.Jacobian = 'off'; 181 end 182 if ~isoctave 183 [x,fval,exitval,output] = fsolve(func,x,options4fsolve,varargin{:}); 184 else 185 % Under Octave, use a wrapper, since fsolve() does not have a 4th arg 186 if ischar(func) 187 func2 = str2func(func); 188 else 189 func2 = func; 190 end 191 func = @(x) func2(x, varargin{:}); 192 % The Octave version of fsolve does not converge when it starts from the solution 193 fvec = feval(func,x); 194 if max(abs(fvec)) >= tolf 195 [x,fval,exitval,output] = fsolve(func,x,options4fsolve); 196 else 197 exitval = 3; 198 end 199 end 200 201 if exitval == 1 202 info = 0; 203 elseif exitval > 1 204 if ischar(func) 205 func2 = str2func(func); 206 else 207 func2 = func; 208 end 209 func = @(x) func2(x, varargin{:}); 210 fvec = feval(func,x); 211 if max(abs(fvec)) >= tolf 212 info = 1; 213 else 214 info = 0; 215 end 216 else 217 info = 1; 218 end 219elseif options.solve_algo == 1 220 [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options.gstep, ... 221 tolf,tolx, ... 222 maxit,options.debug,varargin{:}); 223elseif options.solve_algo == 9 224 [x,info]=trust_region(func,x,1:nn,1:nn,jacobian_flag,options.gstep, ... 225 tolf,tolx, ... 226 maxit,options.debug,varargin{:}); 227elseif options.solve_algo == 2 || options.solve_algo == 4 228 229 if options.solve_algo == 2 230 solver = @solve1; 231 else 232 solver = @trust_region; 233 end 234 235 if ~jacobian_flag 236 fjac = zeros(nn,nn) ; 237 dh = max(abs(x),options.gstep(1)*ones(nn,1))*eps^(1/3); 238 for j = 1:nn 239 xdh = x ; 240 xdh(j) = xdh(j)+dh(j) ; 241 fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ; 242 end 243 end 244 245 [j1,j2,r,s] = dmperm(fjac); 246 247 if options.debug 248 disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]); 249 end 250 251 for i=length(r)-1:-1:1 252 if options.debug 253 disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]); 254 end 255 [x,info]=solver(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ... 256 options.gstep, ... 257 tolf,tolx, ... 258 maxit,options.debug,varargin{:}); 259 if info 260 return 261 end 262 end 263 fvec = feval(func,x,varargin{:}); 264 if max(abs(fvec)) > tolf 265 [x,info]=solver(func,x,1:nn,1:nn,jacobian_flag, ... 266 options.gstep, tolf,tolx, ... 267 maxit,options.debug,varargin{:}); 268 end 269elseif options.solve_algo == 3 270 if jacobian_flag 271 [x,info] = csolve(func,x,func,1e-6,500,varargin{:}); 272 else 273 [x,info] = csolve(func,x,[],1e-6,500,varargin{:}); 274 end 275 [fvec, fjac] = feval(func, x, varargin{:}); 276elseif options.solve_algo == 10 277 % LMMCP 278 olmmcp = options.lmmcp; 279 280 [x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:}); 281 if exitflag == 1 282 info = 0; 283 else 284 info = 1; 285 end 286elseif options.solve_algo == 11 287 % PATH mixed complementary problem 288 % PATH linear mixed complementary problem 289 if ~exist('mcppath') 290 error(['PATH can''t be provided with Dynare. You need to install it ' ... 291 'yourself and add its location to Matlab/Octave path before ' ... 292 'running Dynare']) 293 end 294 omcppath = options.mcppath; 295 global mcp_data 296 mcp_data.func = func; 297 mcp_data.args = varargin; 298 info=0; 299 try 300 [x,fval,jac,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0); 301 catch 302 info = 1; 303 end 304else 305 error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11]') 306end 307