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