gekko - 如何设置 GEKKO 以从多个独立的数据集进行参数估计？

Question

我正在学习如何使用 GEKKO 进行基于实验室批式反应器数据的动力学参数估计，该数据基本上由 A、C 和 P 三种物质的浓度分布组成。就我的问题而言，我使用的是我之前的模型在与单个数据集的参数估计相关的问题中出现。

我的最终目标是能够使用多个实验运行进行参数估计，利用可能在不同温度、物种浓度等条件下收集的数据。由于单个间歇反应器实验的独立性，每个数据集都包含在不同温度下收集的样本时间点。这些不同的时间点（例如将来，不同的温度）对我来说很难实现到 GEKKO 模型中，因为我以前使用实验数据收集时间点作为 GEKKO 模型的 m.time 参数。（代码见帖子末尾）我过去用 gPROMS 和 Athena Visual Studio 解决了这样的问题。

为了说明我的问题，我通过将噪声引入物种浓度曲线并稍微改变实验时间点，从我的原始模型生成了一个人工数据集“实验”数据。然后，我将同一实验物种的所有数据集组合成具有多列的新数组。我这里的想法是，GEKKO 会使用数组对应的每一列的实验数据来进行参数估计，所以这times_comb[:,0]与 相关，A_comb[:,0]而times_comb[:,1]与 相关A_comb[:,1]。

当我尝试运行 GEKKO 模型时，系统确实获得了参数估计的解决方案，但我不清楚问题解决方案是否合理，因为我注意到 GEKKO 变量 A、B、C 和 P 为 34元素向量，它是每个实验数据集中元素的两倍。我认为 GEKKO 在模型设置期间以某种方式结合了时间列和参数向量，从而导致了这 34 个元素变量？我还担心在每个输入参数的列组合过程中，某个时间点与收集的物种信息之间的关系会丢失。

考虑到每个数据集的时间点可能不同，如何改进 GEKKO 可以同时用于参数估计的多个数据集的使用？我查看了 GEKKO 文档示例以及 APMonitor 网站，但我找不到可用于指导的具有多个数据集的示例，因为我对 GEKKO 包相当陌生。

感谢您花时间阅读我的问题以及您可能有的任何帮助/想法。

下面的代码：

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

#Experimental data
times  = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
                      0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
                      0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
                      0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
                      0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
                      0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
                      0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
                      0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
                      0.303198, 0.277822, 0.284194, 0.301471])

#Generate second set of 'experimental data'
times_new = times + np.random.uniform(0.0,0.01)
P_obs_noisy = P_obs+np.random.normal(0,0.05,P_obs.shape)
A_obs_noisy = A_obs+np.random.normal(0,0.05,A_obs.shape)
C_obs_noisy = A_obs+np.random.normal(0,0.05,C_obs.shape)

#Combine two data sets into multi-column arrays
times_comb = np.array([times, times_new]).T
P_comb = np.array([P_obs, P_obs_noisy]).T
A_comb = np.array([A_obs, A_obs_noisy]).T
C_comb = np.array([C_obs, C_obs_noisy]).T

m = GEKKO(remote=False)

t = m.time = times_comb #using two column time array

Am = m.Param(value=A_comb) #Using the two column data as observed parameter
Cm = m.Param(value=C_comb)
Pm = m.Param(value=P_comb)

A = m.Var(1, lb = 0)
B = m.Var(0, lb = 0)
C = m.Var(0, lb = 0)
P = m.Var(0, lb = 0)

k = m.Array(m.FV,6,value=1,lb=0)  
for ki in k:
    ki.STATUS = 1
k1,k2,k3,k4,k5,k6 = k

r1 = m.Var(0, lb = 0)
r2 = m.Var(0, lb = 0)
r3 = m.Var(0, lb = 0)
r4 = m.Var(0, lb = 0)
r5 = m.Var(0, lb = 0)
r6 = m.Var(0, lb = 0)
   
m.Equation(r1 == k1 * A)
m.Equation(r2 == k2 * A * B)
m.Equation(r3 == k3 * C * B)
m.Equation(r4 == k4 * A)
m.Equation(r5 == k5 * A)
m.Equation(r6 == k6 * A * B)

#mass balance diff eqs, function calls rxn function 
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6)
m.Equation(B.dt() ==  r1 - r2 - r3 - r6)
m.Equation(C.dt() ==  r2 - r3 + r4)
m.Equation(P.dt() ==  r3 + r5 + r6)

m.Minimize((A-Am)**2)
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)

m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.NODES = 6
m.solve()

k_opt = []
for ki in k:
    k_opt.append(ki.value[0])
print(k_opt)

plt.plot(t,A)
plt.plot(t,C)
plt.plot(t,P)
plt.plot(t,B)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
plt.plot(times_new, A_obs_noisy,'b*')
plt.plot(times_new, C_obs_noisy,'g*')
plt.plot(times_new, P_obs_noisy,'r*')
plt.show()

标签： gekko

Answer 1

要拥有具有不同时间和数据点的多个数据集，您可以将数据集连接为数据pandas框。这是一个简单的例子：

# data set 1
t_data1 = [0.0,  0.1,  0.2, 0.4, 0.8, 1.00]
x_data1 = [2.0,  1.6,  1.2, 0.7, 0.3, 0.15]

# data set 2
t_data2 = [0.0,  0.15, 0.25, 0.45, 0.85, 0.95]
x_data2 = [3.6,  2.25, 1.75, 1.00, 0.35, 0.20]

合并NaN后的数据有缺失数据的地方：

       x1    x2
Time           
0.00  2.0  3.60
0.10  1.6   NaN
0.15  NaN  2.25
0.20  1.2   NaN
0.25  NaN  1.75

记下数据缺失的地方，其中 1 = 已测量，0 = 未测量。

# indicate which points are measured
z1 = (data['x1']==data['x1']).astype(int) # 0 if NaN
z2 = (data['x2']==data['x2']).astype(int) # 1 if number

最后一步是设置 Gekko 变量、方程和目标以适应数据集。

xm = m.Array(m.Param,2)
zm = m.Array(m.Param,2)
for i in range(2):
    m.Equation(x[i].dt()== -k * x[i])  # differential equations
    m.Minimize(zm[i]*(x[i]-xm[i])**2)  # objectives

您还可以使用 计算初始条件m.free_initial(x[i])。k这为2 个数据集上的一个参数值 ( ) 提供了最佳解决方案。这种方法可以扩展到多个变量或具有不同时间的多个数据集。

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

# data set 1
t_data1 = [0.0,  0.1,  0.2, 0.4, 0.8, 1.00]
x_data1 = [2.0,  1.6,  1.2, 0.7, 0.3, 0.15]

# data set 2
t_data2 = [0.0,  0.15, 0.25, 0.45, 0.85, 0.95]
x_data2 = [3.6,  2.25, 1.75, 1.00, 0.35, 0.20]

# combine with dataframe join
data1 = pd.DataFrame({'Time':t_data1,'x1':x_data1})
data2 = pd.DataFrame({'Time':t_data2,'x2':x_data2})
data1.set_index('Time', inplace=True)
data2.set_index('Time', inplace=True)
data = data1.join(data2,how='outer')
print(data.head())

# indicate which points are measured
z1 = (data['x1']==data['x1']).astype(int) # 0 if NaN
z2 = (data['x2']==data['x2']).astype(int) # 1 if number

# replace NaN with any number (0)
data.fillna(0,inplace=True)

m = GEKKO(remote=False)

# measurements
xm = m.Array(m.Param,2)
xm[0].value = data['x1'].values
xm[1].value = data['x2'].values

# index for objective (0=not measured, 1=measured)
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z2

m.time = data.index
x = m.Array(m.Var,2)                   # fit to measurement
x[0].value=x_data1[0]; x[1].value=x_data2[0]

k = m.FV(); k.STATUS = 1               # adjustable parameter
for i in range(2):
    m.free_initial(x[i])               # calculate initial condition
    m.Equation(x[i].dt()== -k * x[i])  # differential equations
    m.Minimize(zm[i]*(x[i]-xm[i])**2)  # objectives

m.options.IMODE = 5   # dynamic estimation
m.options.NODES = 2   # collocation nodes
m.solve(disp=True)    # solve
k = k.value[0]
print('k = '+str(k))

# plot solution
plt.plot(m.time,x[0].value,'b.--',label='Predicted 1')
plt.plot(m.time,x[1].value,'r.--',label='Predicted 2')
plt.plot(t_data1,x_data1,'bx',label='Measured 1')
plt.plot(t_data2,x_data2,'rx',label='Measured 2')
plt.legend(); plt.xlabel('Time'); plt.ylabel('Value')
plt.xlabel('Time'); 
plt.show()