微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

Python GEKKO:求解模型时参数值的变化

如何解决Python GEKKO:求解模型时参数值的变化

我在使用 GEKKO 时遇到以下问题:在求解模型时,某些参数 (.Param) 发生了变化(其他则没有),我无法确定原因。

背景:我目前正在尝试将代码从 EViews(参见 gennaro.zezza.it)翻译成 Python。我使用 GEKKO 来模拟一个由 11 个方程组成的系统(目前)。我确实想使用参数(而不是似乎工作得很好的常量),因为我需要(“外生”)随时间改变它们的值(因此需要一个数组)。

示例:在以下示例中,“经济体系”对新的政府支出做出反应。在这里,我特别面临“m.alpha1”和“m.alpha2”的问题——如果它们被引入为“.Param”,那么在求解模型时,它们的值将变为 1.0(而不是 0.6 和 0.4)。我怎样才能阻止 GEKKO 这样做? (同样,我希望能够在时间 x 之后将 alpha1 更改为 0.7。例如,下限和上限在这里无济于事。)

感谢您的帮助!!

代码

const int *Get_ar() const

输出 #1:将 m.alpha1 和 m.alpha2 作为 .const

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go

# Initialize model
m = GEKKO(remote=False)
tstart = 1945
tend = 2000
tdur = tend-tstart+1
m.time = np.linspace(0,tend-tstart,tdur)

# Model parameters
m.t = m.Param(value=m.time)

# Exogenous parameters
alpha1_ex = 0.6
alpha2_ex = 0.4
theta_ex = 0.2
w_ex = 1

# -as .Const
m.alpha1 = m.Const(value=alpha1_ex,name='Propensity to consume out of income')
m.alpha2 = m.Const(value=alpha2_ex,name='Propensity to consume out of wealth')
#m.theta = m.Const(value=theta_ex,name='Tax rate')
#m.w = m.Const(value=w_ex,name='Wage rate')

# -as .Param: issues with alpha1 & alpha2
#m.alpha1 = m.Param(value=np.full(tdur,alpha1_ex),name='Propensity to consume out of income')
#m.alpha2 = m.Param(value=np.full(tdur,alpha2_ex),name='Propensity to consume out of wealth')
m.theta = m.Param(value=np.full(tdur,theta_ex),name='Tax rate')
m.w = m.Param(value=np.ones(tdur),name='Wage rate')

# no issues with g_d
m.g_d = m.Param(value=np.zeros(tdur),name='Government goods,demand')
m.g_d[1:] = 20

# Endogenous variables
m.c_d = m.Var(value=0,name='Consumption goods demand by households')
m.c_s = m.Var(value=0,name='Consumption goods supply')
m.g_s = m.Var(value=0,supply')
m.h_h = m.Var(value=0,name='Cash money held by households')
m.h_s = m.Var(value=0,name='Cash money supplied by government')
m.n_d = m.Var(value=0,name='Demand for labor')
m.n_s = m.Var(value=0,name='Supply for labor')
m.t_d = m.Var(value=0,name='Taxes,"demand"')
m.t_s = m.Var(value=0,"supply"')
m.y = m.Var(value=0,name='Income (=GDP)')
m.yd = m.Var(value=0,name='disposable income of households')

# Lag variables
m.h_h_lag = m.Var(value=0,name='Cash money held by households (t-1)')
m.delay(m.h_h,m.h_h_lag,1) # m.h_h_lag = m.h_h(t-1)
m.h_s_lag = m.Var(value=0,name='Cash money supplied by government (t-1)')
m.delay(m.h_s,m.h_s_lag,1)

# Equations
m.Equation(m.c_s == m.c_d)
m.Equation(m.g_s == m.g_d)
m.Equation(m.t_s == m.t_d)
m.Equation(m.n_s == m.n_d)
m.Equation(m.yd == m.w*m.n_s - m.t_s)
m.Equation(m.t_d == m.theta*m.w*m.n_s)
m.Equation(m.c_d == m.alpha1*m.yd + m.alpha2*m.h_h_lag)
m.Equation(m.h_s == m.h_s_lag + m.g_d - m.t_d)
m.Equation(m.h_h == m.h_h_lag + m.yd - m.c_d)
m.Equation(m.y == m.c_s + m.g_s)
m.Equation(m.n_d == m.y/m.w)

# Solve
m.options.IMODE = 4
m.solve(disp=False)

print("Alpha1 = ",m.alpha1.value)
print("Alpha2 = ",m.alpha2.value)
print("Theta = ",m.theta.value)
print("w = ",m.w.value)

# Plot results
fig,axes = plt.subplots(2,2,sharex=True,figsize=(8,7))
fig.canvas.manager.set_window_title('figures Chapter 3')
fig.suptitle('SIM Model - basic')
x_major_ticks = np.arange(0,tdur,5)

axes[0,0].plot(m.time,m.g_d.value,'-',color='black',linewidth=1)
axes[0,0].legend([m.g_d.name],loc=4,fontsize=7)
axes[0,0].grid()
axes[0,0].set_xticks(x_major_ticks)

axes[1,m.y.value,color='red',linewidth=1)
axes[1,0].legend([m.y.name],fontsize=7)
axes[1,0].grid()
axes[1,0].set_xlabel('Time (years)')
axes[1,0].set_xticks(x_major_ticks)

axes[0,1].plot(m.time,m.c_d.value,color='blue',linewidth=0.75)
axes[0,m.yd.value,color='green',1].legend([m.c_d.name,m.yd.name],1].grid()
axes[0,1].set_xticks(x_major_ticks)

ln1 = axes[1,m.h_h.value,color='purple',linewidth=0.75)
axes[1,1].tick_params(axis='y',labelcolor='purple')
ax2 = axes[1,1].twinx()
ln2 = ax2.plot(m.time,[a_i - b_i for a_i,b_i in zip(m.h_h,m.h_h_lag)],color='orange',linewidth=0.75)
ax2.tick_params(axis='y',labelcolor='orange')
lns = ln1+ln2
axes[1,1].legend(lns,[m.h_h.name,'Household savings'],1].grid()
axes[1,1].set_xticks(x_major_ticks)
axes[1,1].set_xlabel('Time (years)')

plt.show()

输出 #2:将 m.alpha1 作为 .param

Alpha1 =  0.6
Alpha2 =  0.4
Theta =  [0.2,0.2,0.2]
w =  [1.0,1.0,1.0]

解决方法

问题在于变量 name='Propensity to consume out of income' 的名称长度超过 25 个字符。

m.alpha1 = m.Param(value=np.full(tdur,alpha1_ex),name='Propensity to consume out of income')
m.alpha2 = m.Param(value=np.full(tdur,alpha2_ex),name='Propensity to consume out of wealth')

模型文件生成正确 (gk_model0.apm),但数据文件 (gk_model0.csv) 标题被截断为 25 个字符。这些文件可以通过 m.open_folder() 访问。错误出现在 gk_write_files.py 的这一行中,其中数字作为长度为 25 的字符串输出。

np.savetxt(os.path.join(self._path,file_name),csv_data.T,delimiter=",",fmt='%1.25s')

我已将此添加为 bug report with tracking on GitHub。一种解决方法是使用较短的变量名称或不使用变量名称。

m.alpha1 = m.Param(value=np.full(tdur,alpha1_ex)) # Propensity to consume out of income

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。