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