3

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
Econner93
  • 33
  • 3

1 Answers1

1

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.

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.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • John - thanks for your time in answering my question. Just wondered if you know the reason why I need to bound the i variable? Best, Alex – Econner93 Mar 25 '20 at 20:50
  • 1
    The bounds help the solver find a solution. It helps the solver avoid guessing values such as `i=100` and then getting a `NaN` evaluation (and final objective) from the equation `J.dt() == m.log(x**ALPHA-i)`. You can adjust the bounds on `i` if you want it to be able to go negative. If you set the upper limit too high, it may get into `NaN` evaluations and the solver could get lost. – John Hedengren Mar 25 '20 at 20:56