python - 在 Python 中使用 GEKKO 的 ODE 系统中的数组变量

标签 python gekko

我之前曾发布过有关在参数估计领域使用 GEKKO 的问题。问题here对于理解 GEKKO 的一些工作原理非常有用，但我没有利用的一个方面是数组变量的使用。在重新审视这个问题时，我意识到我的“测量”数据是错误的，因为我将初始(和未测量)条件作为测量变量的一部分作为数组的第一个元素。 (我很抱歉，因为我不确定如何使 python 语法突出显示发生)


        #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])

在重新审视这个方面时，我遇到了一个问题，我必须向 GEKKO 指定一个积分周期，其中必须包括时间 = 0。模拟与测量数据的时间跨度的这种变化导致了我的问题由于项的长度不同，目标函数不起作用:

    m.Minimize((A-Am)**2) #Array A has 17 elements, array Am has 16
    m.Minimize((P-Pm)**2)
    m.Minimize((C-Cm)**2)

为了尝试解决这个问题，我认为我需要使用 m.Array() 来实现我的变量，这样我就可以按元素对它们进行索引，以期设置以下形式的目标函数:

    m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))

完整的代码显示在帖子的末尾，但我的问题似乎是在设置我的系统的 m.Equations() 方面。目前，我不使用索引来引用变量的各个元素:


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.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6 )

AttributeError: 'numpy.ndarray' object has no attribute 'dt'

当我尝试将它们设置为数组时，Python 卡住了，因为我认为我一定生成了太多方程或其他一些内存问题。也许它试图单独区分索引中的每个方程，这也是不正确的:


m.Equation([r1[k] == k1 * A[k] for k in range(len(r1))])
m.Equation([r2[k] == k2 * A[k] * B[k] for k in range(len(r2))])
m.Equation([r3[k] == k3 * C[k] * B[k] for k in range(len(r3))])
m.Equation([r4[k] == k4 * A[k] for k in range(len(r4))])
m.Equation([r5[k] == k5 * A[k] for k in range(len(r5))])
m.Equation([r6[k] == k6 * A[k] * B[k] for k in range(len(r6))])

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

如果有人能解释如何在方程、特别是 ODE 中使用数组变量，我将不胜感激，因为我找不到任何数组变量是 ODE 系统一部分的例子(我一定在许多例子中错过了它)例子)。此外，如果有另一种也许更简单的方法，用不同长度的变量编写我的目标函数(并且不必像我最初那样将它们指定为数组变量)，知道这一点也很好。

完整脚本如下:

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

data = {'times':[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':[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':[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':[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]}

df = pd.DataFrame(data)
df.set_index('times',inplace = True)



m = GEKKO(remote=False)

time_obs = df.index

t = m.time = np.append(0,time_obs) #adding a time = 0 start value to time discretization for integration

Am= m.Array(m.Param,(len(time_obs)))
Cm= m.Array(m.Param,(len(time_obs)))
Pm= m.Array(m.Param,(len(time_obs)))

for i in range(len(time_obs)):
    
    Am[i].value = df.iloc[i]['A_obs']
    Cm[i].value = df.iloc[i]['C_obs']
    Pm[i].value = df.iloc[i]['P_obs']

A = m.Array(m.Var,(len(t)),lb=0,value=0.0)
A[0].value=1.0
B = m.Array(m.Var,(len(t)),lb=0,value=0.0)
C = m.Array(m.Var,(len(t)),lb=0,value=0.0)
P = m.Array(m.Var,(len(t)),lb=0,value=0.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.Array(m.Var,(len(t)),lb=0,value=0.0)
r2 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r3 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r4 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r5 = m.Array(m.Var,(len(t)),lb=0,value=0.0)
r6 = m.Array(m.Var,(len(t)),lb=0,value=0.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.Equation([r1[k] == k1 * A[k] for k in range(len(r1))])
m.Equation([r2[k] == k2 * A[k] * B[k] for k in range(len(r2))])
m.Equation([r3[k] == k3 * C[k] * B[k] for k in range(len(r3))])
m.Equation([r4[k] == k4 * A[k] for k in range(len(r4))])
m.Equation([r5[k] == k5 * A[k] for k in range(len(r5))])
m.Equation([r6[k] == k6 * A[k] * B[k] for k in range(len(r6))])

#mass balance diff eqs, function calls rxn function 
m.Equation([A[j].dt() == - r1[j] - r2[j] - r4[j] - r5[j] - r6[j] for j in range(len(A))] )
m.Equation([B[j].dt() ==  r1[j] - r2[j] - r3[j] - r6[j] for j in range(len(B))] )
m.Equation([C[j].dt() ==  r2[j] - r3[j] + r4[j] for j in range(len(C))])
m.Equation([P[j].dt() ==  r3[j] + r5[j] + r6[j] for j in range(len(P)) ])
'''
m.Obj(np.sum([(A[i+1]-Am[i])**2 for i in range(len(Am))]))
m.Obj(np.sum([(C[i+1]-Cm[i])**2 for i in range(len(Cm))]))
m.Obj(np.sum([(P[i+1]-Pm[i])**2 for i in range(len(Pm))]))

m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-6
m.options.OTOL = 1E-6
m.options.NODES = 6
m.options.DIAGLEVEL = 10
#m.open_folder()
m.solve()

最佳答案

使用数组变量没有问题，但这需要 collocation equations 的自定义实现将导数与变量联系起来。为缺失的测量放入占位符并创建新参数 mobj 来忽略某些测量(例如初始测量)要容易得多。

# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)

数组对于压缩代码很有用，例如变量的定义。

A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)

解决方案如下:

Solution 

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

data = {'times':[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':[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':[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':[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]}

df = pd.DataFrame(data)
df.set_index('times',inplace = True)

m = GEKKO(remote=False)

time_obs = df.index

t = m.time = time_obs

# ignore first measurement
tobj = np.ones(len(t))
tobj[0] = 0
mobj = m.Param(tobj)

Am= m.Param(df['A_obs'].values)
Cm= m.Param(df['C_obs'].values)
Pm= m.Param(df['P_obs'].values)

A,B,C,P = m.Array(m.Var,4,lb=0,fixed_initial=False)

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

r = m.Array(m.Var,6,lb=0,value=0.0)
r1,r2,r3,r4,r5,r6 = r
   
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(mobj*(A-Am)**2)
m.Minimize(mobj*(C-Cm)**2)
m.Minimize(mobj*(P-Pm)**2)

m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-6
m.options.OTOL = 1E-6
m.options.NODES = 3
m.options.DIAGLEVEL = 0
#m.open_folder()
m.solve()

import matplotlib.pyplot as plt
plt.plot(m.time,Am.value,'ro')
plt.plot(m.time,A.value,'r-',label='A')
plt.plot(m.time,Cm.value,'bo')
plt.plot(m.time,C.value,'b-',label='C')
plt.plot(m.time,Pm.value,'go')
plt.plot(m.time,P.value,'g-',label='P')
plt.xlabel('Time')
plt.legend()
plt.show()

关于python - 在 Python 中使用 GEKKO 的 ODE 系统中的数组变量，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/68927734/

上一篇：flutter - 如何解决构建运行程序 : "This requires the ' non-nullable' language feature to be enabled. 的问题？

下一篇：python - 在pygame中旋转一个矩形

相关文章：

python - gekko中定义最大函数时出错

python - 在 Python 中使用 Gekko 求解函数

gekko - 如何将 GEKKO 返回的整数转换回 python 可识别的整数？

python-3.x - 评估 gekko 中的表达式

python - 使用 python 在有限定义区间上数值求解非线性方程组

python - 为什么要在 Python 中设置部分对象的 __doc__？

python - 为什么 ast.literal_eval ('5 * 7' ) 失败了？

python - 如何在Python中使用 `is None`编写一行三元表达式？

python 使用 bool 表达式通过 While 循环遍历列表的元素

python - Django:不推荐使用对 url() 的字符串 View 参数的支持，并将在 Django 1.10 中删除

©2024 IT工具网  联系我们