用Pandas Datafram绘制Gekko模型

Question

我正在尝试将理论模型与历史数据进行比较。我的计划是使用Gekko将真实世界的数据从几个pandas数据帧输入到方程中(跳到# equation #):

import pandas as pd
import pylab as plt
import numpy as np
from gekko import GEKKO
#DATA#
#oil supply
data = pd.read_html('https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=pet&s=mcrfpus2&f=m')
s=data[4]
sply = s.loc[s['Year'] >= 1984]
supplyvalues = sply.filter(items = ['Jan', 'Feb', 'Mar', 'Apr', 'May','Jun','Jul','Aug','Sep', 'Oct', 'Nov','Dec'])
val2 = supplyvalues.dropna()
protoval3 = supplyvalues.shift(1)
val3=protoval3.dropna()
#price of oil
p = pd.read_html('https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=PET&s=EMA_EPM0_PTG_NUS_DPG&f=M')
Op=p[4]
prc = Op.loc[Op['Year'] >= 1984]
pricevalues = prc.filter(items = ['Jan', 'Feb', 'Mar', 'Apr', 'May','Jun','Jul','Aug','Sep', 'Oct', 'Nov','Dec'])
val1 = pricevalues.dropna()
#EQUATION#
m=GEKKO()
tm = np.linspace(0,10,50)
t=m.Param(value = tm)
m.time = tm
m.options.IMODE=4
#forming equation
a,b,c = m.Array(m.Param,3)
a.value = val1.values
b.value = val2.values
c.value = val3.values
x = m.Var()
m.Equation(a == c-a*b)
m.solve(disp=False)

问题是，每次我尝试，我得到这个错误:

Exception: Data arrays must have the same length, and match time discretization in dynamic problems

我检查了数据帧的形状，它们是相等的，所以我不确定出了什么问题。有点烦，因为我觉得我真的很亲近。我将非常感激你能给予的任何帮助!

p。s:我还是一个编程初学者，所以如果你能描述一下你的新代码行是做什么的，我会很感激的(如果你有的话)。这是完全可选的!

Answer 1

其中val1、val2、val3为38x12矩阵。

Jan    Feb    Mar    Apr    May  ...    Aug    Sep    Oct    Nov    Dec
1   0.906  0.902  0.907  0.929  0.934  ...  0.892  0.897  0.905  0.899  0.880
3   0.846  0.836  0.871  0.924  0.944  ...  0.940  0.919  0.908  0.917  0.919
4   0.893  0.805  0.654  0.591  0.638  ...  0.555  0.562  0.532  0.532  0.542
5   0.597  0.621  0.627  0.649  0.663  ...  0.716  0.705  0.697  0.694  0.674
6   0.649  0.633  0.625  0.660  0.684  ...  0.718  0.700  0.680  0.676  0.661

一个成功的解决方案需要做一些修改。

1. 仅选择一个月或将年度数据与时间离散点对齐。

#forming equation
a,b,c = m.Array(m.Param,3)
a.value = val1['Jan'].values
b.value = val2['Jan'].values
c.value = val3['Jan'].values

如果需要多个月的模拟，则将每个参数放入值向量中。

2. 将时间点与val1,val2和val3的输入对齐。

tm = np.linspace(0,10,38)

3. 核对公式

每个月都有38时间点，所以时间离散化也需要这么多时间点。这个方程看起来并不完整——动态模拟通常有像x.dt()==-x + a - b + c这样的微分方程。

4. 导入需要lxml包。安装:

pip install lxml

这是一个完整的脚本，解决成功，但不是预期的模拟。也许这是开发模拟问题的一个很好的起点。

import pandas as pd
import pylab as plt
import numpy as np
from gekko import GEKKO
#DATA#
#oil supply
data = pd.read_html('https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=pet&s=mcrfpus2&f=m')
s=data[4]
sply = s.loc[s['Year'] >= 1984]
supplyvalues = sply.filter(items = ['Jan', 'Feb', 'Mar', 'Apr', 'May','Jun','Jul','Aug','Sep', 'Oct', 'Nov','Dec'])
val2 = supplyvalues.dropna()
protoval3 = supplyvalues.shift(1)
val3=protoval3.dropna()
#price of oil
p = pd.read_html('https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=PET&s=EMA_EPM0_PTG_NUS_DPG&f=M')
Op=p[4]
prc = Op.loc[Op['Year'] >= 1984]
pricevalues = prc.filter(items = ['Jan', 'Feb', 'Mar', 'Apr', 'May','Jun','Jul','Aug','Sep', 'Oct', 'Nov','Dec'])
val1 = pricevalues.dropna()
print(np.shape(val1.values))
#EQUATION#
m=GEKKO(remote=False)
tm = np.linspace(0,10,38)
print(len(tm))
t=m.Param(value = tm)
m.time = tm
m.options.IMODE=4
#forming equation
a,b,c = m.Array(m.Param,3)
a.value = val1['Jan'].values
b.value = val2['Jan'].values
c.value = val3['Jan'].values
x = m.Var()
m.Equation(x == c-a*b)
m.solve(disp=True)