function [tt,xx,uu,val] = basicgrowthmodel(iv,N,T,constr,mpciter,myprint) % This is the basic growth model from Section 4.1 in L. Gruene, W. Semmler % and M. Stieler: "Using Nonlinear Model Predictive Control for Dynamic % Decision Problems in Economics", 2015 % This file needs nmpc.m. Both files must be contained in the same folder. % (c) Gruene, Semmler, Stieler 2015 % Input: iv - initial value of the system % N - optimization horizon % T - sampling time (needed for continuous-time models) % constr - state constraints given as [lower upper] % mpciter - number of iterations to be performed by the % algorithm % myprint - 0: no plot, 1: plot closed loop, 2: plot open and % closed loop % Output: tt - array of times % xx - closed loop trajectory % uu - mpc feedback % val - accumulated cost along the mpc closed loop % trajectory % Example calls: % This call will generate Figure 4.1 left % [tt,xx,uu,val] = basicgrowthmodel(5,5,1,[0.1 10],30,2); % % A smaller horizon (N=2), plot only closed loop % [tt,xx,uu,val] = basicgrowthmodel(5,2,1,[0.1 10],30,1); % % Or maybe you would like to try another initial value (x0 = 0.5) and % observe the turnpike property % [tt,xx,uu,val] = basicgrowthmodel(0.5,5,1,[0.1 10],30,2); % % If you're juste interested in the accumulated cost along the trajectory % [tt,xx,uu,val] = basicgrowthmodel(5,5,1,[0.1 10],30,0); % -val % Parameters of the model global A alpha beta; A = 5.0; alpha = 0.34; beta = 0.95; % Do not change the following three lines mpciterations = mpciter; t0 = 0.0; x0 = iv; % Initial guess for fmincon, sometimes needs to be adjusted u0 = ones(1,N); % Options for fmincon tol_opt = 1e-12; opt_option = 0; iprint = 5; % Model in discrete or continuous time? type = 'difference equation'; % 'differential equation' % Options for the ode solver (if used) atol_ode_real = 1e-12; rtol_ode_real = 1e-12; atol_ode_sim = 1e-4; rtol_ode_sim = 1e-4; global gN; gN = N; global gconstr; gconstr = constr; global gprint; gprint = myprint; global guu; guu = []; global gtt; gtt = [ t0 ]; global gxx; gxx = [ x0 ]; global gval; % Here, the NMPC algorithm is executed nmpc(@runningcosts, @terminalcosts, @constraints, ... @terminalconstraints, @linearconstraints, @system, ... mpciterations, N, T, t0, x0, u0, ... tol_opt, opt_option, ... type, atol_ode_real, rtol_ode_real, atol_ode_sim, rtol_ode_sim, ... iprint, @printHeader, @printClosedloopData, @plotTrajectories); tt = gtt; xx = gxx; uu = guu; val = gval; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Definition of the NMPC functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Definition of the payoff function g(x,u) function cost = runningcosts(t, x, u) global A alpha beta; cost = -beta^t*log(A*x(1)^alpha-u(1)); % The sign of the payoff has to be reversed since the program solves a % minimization problem end % Definition of a terminal cost function (not used in our paper) function cost = terminalcosts(t, x) cost = 0.0; end % State constraints at time t function [c,ceq] = constraints(t, x, u) global gconstr; c(1) = x(1) - gconstr(2); c(2) = -x(1) + gconstr(1); ceq = []; end % State constraints at time N function [c,ceq] = terminalconstraints(t, x) global gconstr; c(1) = x(1) - gconstr(2); c(2) = -x(1) + gconstr(1); ceq = []; end % Control constraints at time t function [A, b, Aeq, beq, lb, ub] = linearconstraints(t, x, u) A = []; b = []; Aeq = []; beq = []; lb = []; ub = []; end % The system (either a differential or a difference equation) function y = system(t, x, u, T) y(1) = u(1); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Definition of output format %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function printHeader() fprintf(' k | u(k) x(k) Time\n'); fprintf('--------------------------------------\n'); end function printClosedloopData(mpciter, u, x, t_Elapsed) fprintf(' %3d | %+11.6f %+11.6f %+6.3f', mpciter, ... u(1,1), x(1), t_Elapsed); end function plotTrajectories(dynamic, system, T, t0, x0, u, ... atol_ode, rtol_ode, type) global gN; global gprint; global gtt; global gxx; global guu; global gval; gtt = [gtt; t0+T]; if (gprint>=1) figure(1); hold on; end; x = x0; gval = 0; for k=1:gN gval = gval + runningcosts(t0+k-1,x,u(:,k)); [x,t_intermediate, x_intermediate] = ... dynamic(system, T, t0+k-1, x, u(:,k), ... atol_ode, rtol_ode, type); if (k==1) if (gprint>=1) plot(t_intermediate,x_intermediate(:,1),'-ok','LineWidth',2); end; gxx = [gxx; x]; guu = [guu, u(1,1)]; else if (gprint>=2) plot(t_intermediate,x_intermediate(:,1),'--k','LineWidth',2); end; end end gval = gval/gN; if (gprint>=1) if (gprint>=2) title(['open loop trajectories (dashed) and closed loop trajectory (solid)']); else title(['closed loop trajectory']); end; xlabel('k'); ylabel('x(k)'); grid on; end; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%