使用 octave/matlab ODE 求解器时拟合 ODE 的参数

Question

我在 Octave 中使用 OdePkg 来解决一个僵硬 ODE 系统，例如通过 ode5r：

function yprime = myODEs(t,Y,param)
    yprime = [
        - param(1) * Y(1);                      # ODE for Y(1)
        param(1) * Y(1) - param(2) Y(2) * Y(3); # ODE for Y(2)
        param(2) Y(2) * Y(3)                    # ODE for Y(3)
                                                # etc.
];

time_span = [1, 24]         # time span of interest
Y0        = [1.0, 1.1, 1.3] # initial values for dependent variables Y
param     = [7.2, 8.6, 9.5] # parameters, to be optimized

[t, Y] = ode5r(@myODEs, time_span, Y0, ..., param);

求解器将因变量 Y 存储在一个关于时间 t（向量）的矩阵中：

t     Y(1)  Y(2)  Y(3)
0.0   1.0   1.1   1.3
0.1   ...   ...   ...
0.5   ...   ...   ...
0.9   ...   ...   ...
...   ...   ...   ...
4.0   ...   ...   ...
...   ...   ...   ...
24.0  ...   ...   ...

我想在 param 中拟合参数，以便生成的变量 Y 最适合我的参考值，例如：

t         Y(1)  Y(2)  Y(3)
0.5       1.1   N/A   N/A
1.0       1.9   N/A   N/A
4.0       2.3   2.7   2.1
5.0       N/A   2.6   2.2
24.0      0.9   1.5   2.0

哪个 Octave/Matlab（欢迎使用其他语言）例程可以执行多参数（最小二乘/样条）拟合？ 如何在拟合中组合不同初始值的参数集 Y0？ 如果你能给我一些提示和可能性，我会很高兴。

最好的问候， 西蒙

Answer 1

这对于 scipy 应该相对简单。 scipy.optimize.leastsq() 采用一个函数，该函数应返回给定参数向量的残差数组。它将最小化残差的平方和。要处理具有不同初始值的多个数据集，您只需为每个数据集运行一次 ODE，计算每个数据集/运行对的残差，然后将残差向量连接在一起。这是一个粗略的草图：

import numpy
from scipy import integrate, optimize

# The initial guess.
p0 = numpy.array([7.2, 8.6, 9.5])

# The collected datasets.
# A list of (t, y0, y) tuples.
# The y's are all (len(y0), len(t))-shaped arrays. The output of
# integrate.odeint is also in this format.
datasets = [...]

def odes(y, t, params):
    dydt = [
        -params[0] * y[0],
        params[0]*y[0] - params[1]*y[1]*y[2],
        params[1]*y[1]*y[2],
    ]
    return np.array(dydt)

def residuals(params, datasets):
    res = []
    for t, y0, y in datasets:
        res.append(integrate.odeint(odes, y0, t, args=(params,)) - y)

    # Stack them horizontally and flatten the array into the expected vector.
    # You're on your own for handling missing data. Look into the numpy.ma
    # module.
    all_residuals = numpy.hstack(res).ravel()
    return all_residuals

opt_params, err = optimize.leastsq(residuals, p0, args=(datasets,))

Answer 2

你的意思是每个函数 y(t) 都需要拟合吗？在这种情况下，每组 Yi 与时间的租赁平方或样条拟合将正常工作。如果不查看您的数据，就无法判断哪个是最好的。

如果您的意思是要在给定时间点的所有 Yi 值上拟合一条曲线，然后观察该曲线随时间的变化，您将不得不提出另一个自变量。

更新：最小二乘拟合就是这样 - 我没有特别推荐的例程。 SciPy has one，我确定。很抱歉我没有更好的推荐。我现在才学 Python。

我不知道您所说的“健康指标”是什么意思。最小二乘拟合计算最小化拟合和每个点的数据之间的误差的均方的系数。

将多个数据集组合成一个拟合的一种方法：合并它们并重新运行计算。

Answer 3

我开发了一个全面的 Matlab 工具箱，用于将 ODE 的参数和初始值拟合到多个实验数据集。它可以根据每个实验处理不同的初始值，可在www.potterswheel.de 获得。