Implementing an macroeconomic model in Python's GEKKO
This question is focussed somewhat on economic optimisation, and somewhat on python implementation, but maybe some in the community are able to help. I'm trying to implement a standard continuous-time macroeconomic savings model in Python's GEKKO platform, but haven't been able to get it to solve. I've taken the economic example provided in GEKKO's documentation, and adapted to the basic savings decision model, but things are not quite working out. The model maximises the sum of utility from consumption, where consumption + investment = output. E.g. max integral(U(y-i)). Output y = k^ALPHA. investment = dk/dt+delta*k.
Can anyone tell why my code can't be solved? Is the platform even capable of solving such a model? I haven't seen many examples of economists using this platform to solve models, but not sure if this is because the platform is not suited or otherwise. It's a great platform and really keen to make it work if possible. Thank you in advance.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
n=501
m.time = np.linspace(0,10,n)
ALPHA,DELTA = 0.333,0.99
i = m.MV(value=0)
i.STATUS = 1
i.DCOST = 0
x = m.Var(value=20,lb=0) # fish population
m.Equation(x.dt() == i-DELTA*x)
J = m.Var(value=0) # objective (profit)
Jf = m.FV() # final objective
Jf.STATUS = 1
m.Connection(Jf,J,pos2='end')
m.Equation(J.dt() == m.log(x**ALPHA-i))
m.Obj(-Jf) # maximize profit
m.options.IMODE = 6 # optimal control
m.options.NODES = 3 # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.solve(disp=True) # Solve
You are getting NaN in the equation dJ/dt = ln(x**ALPHA-i). When you include bounds i>0 and i<1, the solver finds a solution.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO()
n=501
m.time = np.linspace(0,10,n)
ALPHA,DELTA = 0.333,0.99
i = m.MV(value=0,lb=0,ub=1)
i.STATUS = 1
i.DCOST = 0
x = m.Var(value=20,lb=0) # fish population
m.Equation(x.dt() == i-DELTA*x)
J = m.Var(value=0) # objective (profit)
Jf = m.FV() # final objective
Jf.STATUS = 1
m.Connection(Jf,J,pos2='end')
m.Equation(J.dt() == m.log(x**ALPHA-i))
m.Obj(-Jf) # maximize profit
m.options.IMODE = 6 # optimal control
m.options.NODES = 3 # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.solve(disp=True) # Solve
plt.subplot(2,1,1)
plt.plot(m.time,x.value)
plt.ylabel('x')
plt.subplot(2,1,2)
plt.plot(m.time,i.value)
plt.ylabel('i')
plt.show()
Instead of m.Obj() (minimize) you can use the newer functions m.Minimize() or m.Maximize() to clarify the objective function intent. For example, you could switch to m.Maximize(Jf) to make it more readable.
There are also a couple other examples that may help you with integral objectives (see solution 2) and economic dynamic optimization.