Python GEKKO实现消费效用最大化经济模型代码报错排查

问题背景 #

搭建的经济增长最优控制模型如下：

模型核心目标为求解控制变量u，实现跨期总消费效用最大化，基础变量定义：

• u = ln(c(t))，其中c(t)代表t时刻的消费量
• dk/dt为最优控制路径下的资本动态方程，物理意义为总产出扣除当期消费后的净投资。

使用Python的GEKKO库编写实现代码时运行失败，暂未定位错误点，原始代码如下：

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()

代码核心错误点 #

原始代码共存在6处会导致运行失败/结果错误的问题：

1. 变量定义顺序错误：ub_cal = np.log(k)写在k变量定义之前，运行时直接触发NameError，该行执行时k尚未被赋值。
2. 模型方程逻辑错误：资本动态方程扣除项错误。u是对数消费ln(c)，方程中需要扣除的是实际消费量c = exp(u)，而非u本身，属于核心模型设定错误。
3. 函数调用错误：折现项直接调用numpy的np.exp()作用于GEKKO符号变量，会触发类型不匹配错误，需要改用GEKKO内置的m.exp()构造符号表达式。
4. 初始值与边界设置错误：
   • 控制变量u初始值设为101，但初始资本为1000时对数消费最大值仅为ln(1000)≈6.9，初始点严重偏离可行域，求解器无法收敛
   • u的上界不能在初始化时设为静态值，k是随时间变化的状态变量，需要单独添加动态约束
   • 累计效用J初始值设为6.8不符合逻辑，t=0时累计效用应为0
5. 绘图代码错误：子图布局参数冲突（混用2行1列、2行2列的子图索引），还引用了未定义的变量x（实际对应资本变量k）。
6. 调试设置问题：关闭了求解器输出，出现收敛错误时无法快速定位原因，调试阶段建议打开求解器日志。

修正后可运行代码 #

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

# 初始化模型
m = GEKKO(remote=False)
# 时间网格
n = 501
m.time = np.linspace(0, 90, n)
t = m.Param(value=m.time)

# 常数参数
delta = 0.1  # 资本折旧率
rho = 0.03   # 效用折现率，可根据实际模型参数调整
c_min = 100  # 最低消费约束

# 变量定义
k = m.Var(value=1000, lb=1)  # 资本存量，下界设为1避免幂运算数值错误
u = m.MV(value=np.log(500), lb=np.log(c_min))  # 控制变量：对数消费
u.STATUS = 1
u.DCOST = 0

# 动态约束
c = m.exp(u)  # 实际消费量
m.Equation(k.dt() == 10*k**(2/3) - delta*k - c)  # 资本积累方程
m.Equation(u <= m.log(k))  # 消费动态上界约束

# 目标函数：折现总效用最大化
J = m.Var(value=0)  # 累计效用初始值为0
Jf = m.FV()
Jf.STATUS = 1
m.Connection(Jf, J, pos2='end')
m.Equation(J.dt() == m.exp(-rho*t) * u)
m.Maximize(Jf)

# 求解器设置
m.options.IMODE = 6   # 最优控制模式
m.options.NODES = 3   # 配置点数量
m.options.SOLVER = 3  # IPOPT求解器
m.solve(disp=True)

# 结果输出
print('最大总效用值: ' + str(Jf.value[0]))

# 结果可视化
plt.figure(figsize=(10,8))
plt.subplot(3,1,1)
plt.plot(m.time, J.value, 'r-', linewidth=2, label='累计总效用')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time, k.value, 'b-', linewidth=2, label='资本存量')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time, np.exp(u.value), 'k-', linewidth=2, label='消费量')
plt.xlabel('时间（年）')
plt.legend()
plt.tight_layout()
plt.show()

注：原代码折现因子写为exp(-t)对应折现率1，属于极高折现率，修正代码折现率取宏观模型常用的0.03，可根据实际需求调整。

内容的提问来源于stack exchange，提问作者Ліда Хань