Implementing an economic model in Python's GEKKO

user_19324103 0
Viewed 56

I have a model

model

where I should find u, that maximizes total consumption usefulness. u = ln(c(t)), where c(t) is consumption.

dk/dt shows the dynamics of investments (production - consumption) we own when the u is optimized

I have a problem with implementing this model in Python's GEKKO

Here is code I've tried to write, it doesn't work. I don`t know where is my problem

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# create GEKKO model
m = GEKKO()
# time points
n=501
m.time = np.linspace(0,90,n)

# constants
koef = 0.1 # коефієнти 

# керування
lb_cal = np.log(100)
ub_cal = np.log(k)
k = m.Var(value=1000) # інвестиції
u = m.MV(value=101,lb=lb_cal,ub=ub_cal)
u.STATUS = 1
u.DCOST = 0

# investments rate
m.Equation(k.dt() == 10*k**(2/3)-koef*k-u)

J = m.Var(value=6.8) # objective (profit)
Jf = m.FV() # final objective
Jf.STATUS = 1

m.Connection(Jf,J,pos2='end')
m.Equation(J.dt() == np.exp(-m.time)*u)

m.Maximize(Jf) # maximize profit

m.options.IMODE = 6  # optimal control
m.options.NODES = 3  # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.solve(disp=False) # Solve

print('Мах загальної корисності: ' + str(Jf.value[0]))
plt.figure(1) # plot results
plt.subplot(2,1,1)
plt.plot(m.time,J.value,'r--',label='general utility')
plt.legend()
plt.subplot(2,2,1)
plt.plot(m.time,x.value,'b-',label='investments')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,u.value,'k--',label='rate')
plt.xlabel('Time (yr)')
plt.legend()
plt.show()

Would be extra grateful if someone could explain where is my problem

economics gekko mathematical-optimization modeling python
1 Answers

Reformulating the problem gives a successful solution without the upper bound for c(t). Use the gekko function m.log() instead of np.log(). The m.integral() function simplifies the problem. Additional details are available in the gekko documentation.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

# create GEKKO model
m = GEKKO()
# time points
n=101
m.time = np.linspace(0,90,n)
delta = 1
t = m.Param(m.time)

# constants
koef = 0.1 # коефієнти 

# керування
k = m.Var(value=1000) # інвестиції
c = m.MV(lb=100)
c.STATUS = 1; c.DCOST = 0

# investments rate
m.Equation(k.dt() == k - 0.1*k - c)

# maximize last value at T=90
p = np.zeros(n); p[-1] = 1
final = m.Param(p)
J = m.Intermediate(m.integral(m.exp(-delta*t)*m.log(c)))
m.Maximize(final * J) # maximize final profit

m.options.IMODE = 6  # optimal control
m.options.NODES = 3  # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.options.MAX_ITER = 1000
m.solve(disp=True) # Solve

try:
    # add inequality constraint
    m.Equation(c <= 10*k**(2/3) - 0.1*k)

    # solve again with inequality constraint
    m.solve(disp=True) # Solve
except:
    print('no feasible solution with inequality constraint')

print('Мах загальної корисності: ' + str(J.value[-1]))
plt.figure(1) # plot results
plt.subplot(3,1,1)
plt.plot(m.time,J.value,'r--',label='general utility')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time,k.value,'b-',label='consumption')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time,c.value,'k--',label='rate')
plt.xlabel('Time (yr)')
plt.legend()
plt.show()

When adding the upper bound for c(t) as m.Equation(c <= 10*k**(2/3) - 0.1*k), there is no feasible solution. Try adding a less restrictive constraint and tighten the constraint until the problem becomes infeasible. This helps to identify problems with enforcing the constraint.

One other issue to check is that m.exp(-delta*t) quickly becomes zero so the only contributions to the integral occur in the first few years. The value of delta may need to be updated to a lower value (maybe 0.01) to reflect the true depreciation rate.

user_23669410
Related