Python使用GEKKO求解ODE系统时数组变量与目标函数适配问题

问题核心原因与解决方案 #

你混淆了GEKKO动态变量的时间序列属性和数组变量的用途，实际上完全不需要定义数组变量就能解决长度不匹配问题：

• GEKKO中每个普通动态变量（不定义为数组）本身就对应m.time全长度的序列值，支持类似numpy的切片操作，写目标函数时取变量的[1:]切片去掉0点，刚好和16个测量值匹配即可。
• 你之前将A/B/C等定义为和m.time等长的变量数组，相当于创建了17个独立的无时间维度的静态变量，自然没有dt（时间导数）属性，后续试图给每个静态变量写导数方程、逐个写反应方程，会生成上百个冗余方程，自然会报错或者卡死。
• 数组变量的设计用途是定义多主体/多组分的平行变量（比如多个反应器的浓度、多种组分的参数），每个数组元素仍然是一个完整的动态变量，而非某个动态变量的单个时间点值。

修正后可直接运行的代码 #

import numpy as np
import pandas as pd
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)
# 时间轴保留0点，共17个点
m.time = np.append(0, df.index.values)

# 测量值参数，长度16，对应1~16的时间点
Am = m.Param(df['A_obs'].values)
Cm = m.Param(df['C_obs'].values)
Pm = m.Param(df['P_obs'].values)

# 直接定义动态变量，不用数组，给定初始值
A = m.Var(1.0, lb=0)
B = m.Var(0.0, lb=0)
C = m.Var(0.0, lb=0)
P = m.Var(0.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.Intermediate(k1 * A)
r2 = m.Intermediate(k2 * A * B)
r3 = m.Intermediate(k3 * C * B)
r4 = m.Intermediate(k4 * A)
r5 = m.Intermediate(k5 * A)
r6 = m.Intermediate(k6 * A * B)

# 微分方程直接写即可，GEKKO自动处理所有时间点
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)

# 目标函数取变量[1:]切片（去掉0点）和测量值匹配
m.Minimize((A[1:] - Am)**2)
m.Minimize((C[1:] - Cm)**2)
m.Minimize((P[1:] - Pm)**2)

m.options.IMODE = 5
m.options.SOLVER = 3
m.options.RTOL = 1E-6
m.options.OTOL = 1E-6
m.options.NODES = 6
m.solve()

# 输出拟合得到的参数
print(f'拟合参数结果：\nk1={k1.value[0]:.4f}\nk2={k2.value[0]:.4f}\nk3={k3.value[0]:.4f}\nk4={k4.value[0]:.4f}\nk5={k5.value[0]:.4f}\nk6={k6.value[0]:.4f}')

数组变量的正确用法示例 #

如果确实有多平行变量的需求，比如你有4个组分需要统一管理，可以按如下方式使用数组变量，每个元素仍是完整的动态变量：

# 定义4个组分的浓度数组
conc = m.Array(m.Var, 4, lb=0)
# 给定初始值
conc[0].value = 1.0 # A的初始值
conc[1].value = 0.0 # B的初始值
conc[2].value = 0.0 # C的初始值
conc[3].value = 0.0 # P的初始值

# 写方程时用索引访问即可
m.Equation(conc[0].dt() == - r1 - r2 - r4 - r5 - r6)

内容的提问来源于stack exchange，提问作者juanmvenegas