带有额外模拟点的 Gekko 中的不频繁测量

如何解决带有额外模拟点的 Gekko 中的不频繁测量

我正在解决一个更复杂的动态数据协调问题，该问题类似于 APMonitor dynamic optimization course 中发布的示例问题。与这个问题类似，我在不规则的时间间隔很少进行测量。我不想使用更准确的 NODES=5（为了提高准确性），而是希望在 NODES=2 处使用 [0,0.1,0.2,0.3,0.4,0.5,0.6.0.7,0.8,0.9,1.0] 进行更频繁的模拟点，其中还包括 {{} 处目标的不频繁测量1}}。

[0,1]

造成这种情况的原因之一是模拟模型需要比测量更频繁地求解，以保持 from gekko import GEKKO t_data = [0,1] x_data = [2.0,1.6,1.2,0.7,0.15] m = GEKKO(remote=False) m.time = t_data x = m.Var(value=x_data) # fit to measurement xd = m.Param(x_data) k = m.FV(); k.STATUS = 1 # adjustable parameter m.Equation(x.dt()== -k * x) # differential equation m.Minimize((x-xd)**2) m.options.IMODE = 5 # dynamic estimation m.options.NODES = 2 # collocation nodes m.solve(disp=False) # display solver output k = k.value[0] 的积分精度。如果只在测量点求解，精确解与数值解有偏差。

NODES=2

如何在 Gekko 中使用具有更频繁模拟点和不频繁测量的低阶搭配模式？

解决方法

您可以使用 join() 函数将 Pandas 中不频繁的测量与更频繁的模拟步骤结合起来：

from gekko import GEKKO
import numpy as np
import pandas as pd
# measurements
t_data = [0,0.1,0.2,0.4,0.8,1]
x_data = [2.0,1.6,1.2,0.7,0.3,0.15]
df1 = pd.DataFrame({'time':t_data,'x':x_data})
df1.set_index('time',inplace=True)
# simulation time points
df2 = pd.DataFrame({'time':np.linspace(0,1,51)})
df2.set_index('time',inplace=True)
# merge dataframes
df = df2.join(df1,how='outer')
# get True (1) or False (0) for measurement
df['meas'] = (df['x'].values==df['x'].values).astype(int)
# replace NaN with zeros
df0 = df.fillna(value=0)
print(df.head(15))

这给出了缺失数据点的 NaN 组合集。

        x  meas
time           
0.00  2.0     1
0.02  NaN     0
0.04  NaN     0
0.06  NaN     0
0.08  NaN     0
0.10  1.6     1
0.12  NaN     0
0.14  NaN     0
0.16  NaN     0
0.18  NaN     0
0.20  1.2     1
0.22  NaN     0
0.24  NaN     0
0.26  NaN     0
0.28  NaN     0

然后，您可以创建修改后的目标函数，该函数使用 meas 来控制何时仅在定义的点处最小化模拟和测量。

m = GEKKO(remote=False)
m.time = df.index.values
x = m.Var(value=x_data[0])     # fit to measurement
xd = m.Param(df0['x'].values)
meas = m.Param(df0['meas'].values)
k = m.FV(); k.STATUS = 1    # adjustable parameter
m.Equation(x.dt()== -k * x) # differential equation
m.Minimize(meas*(x-xd)**2)
m.options.SOLVER = 1  # APOPT solver
m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 2   # collocation nodes
m.solve(disp=True)   # display solver output
k = k.value[0]

绘图显示了结果。我包含了更多数据点，因为 NODES=2 不是很准确。切换到 NODES=3 比添加更多点更有助于提高准确性。

import numpy as np
import matplotlib.pyplot as plt  # plot solution
plt.plot(m.time,x.value,'bo',\
         label='Predicted (k='+str(np.round(k,2))+')')
plt.plot(m.time,df['x'].values,'rs',label='Measured')
# plot exact solution
t = np.linspace(0,1); xe = 2*np.exp(-k*t)
plt.plot(t,xe,'k:',label='Exact Solution')
plt.legend()
plt.xlabel('Time'),plt.ylabel('Value')
plt.ylim([-0.2,2.2])
plt.show()