numpy 如何用插值法拟合常微分方程组？

ctzwtxfj 于 11个月前发布在其他

关注(0)|答案(3)|浏览(109)

给定数据x0、x1和x2，假设我想在以下常微分方程系统中拟合参数k0、k1、k2、p1和p2
下面的代码实现了我想要的功能

pip install lmfit

个字符

的数据
现在，我不想用一个ODE拟合x0，我想拟合它的一个插值，所以我有更多的自由。理想情况下，我想定义一个（或多个）插值参数，这也将进入优化问题。关于最好的方法有什么想法吗？

**我的想法：**类似于this answer in Mathematica，我们可以考虑定制一个插值函数，它将为连续数据点之间的每个间隔取一个参数，这些参数将与原始系统中的其余参数k1，k2，p1和p2进行位拟合。我不确定这是否理想（或快速），但如果我们忽略x0的第一个方程并使用插值，它将给予更好的近似。例如，使用以下插值的参数化版本

from scipy.interpolate import UnivariateSpline
spline = UnivariateSpline(t_measured, x0_measured)
x0_spline = spline(t_interpolated)
plt.figure(figsize=(8, 4))
plt.scatter(t_measured, x0_measured, color='blue', label='Measured x0', zorder=3)
plt.plot(t_interpolated, x0_spline, color='green', linestyle='-', label='Spline Interpolation of x0')
plt.title("Spline Interpolation of x0")
plt.xlabel("Time")
plt.ylabel("x0")
plt.legend()
plt.grid(True)
plt.show()

型

的

numpy

来源：https://stackoverflow.com/questions/77587951/how-to-fit-a-system-of-odes-with-interpolation

3条答案

按热度按时间

bbuxkriu1#

如果你的数据集可以用你的模型拟合，它可能会错过观测值（curse of dimensionality）。
你的过程可能是正确的，但由于缺少点而无法正确收敛。对于这个具有5个未知参数的3D ODE系统，3x 6点似乎不足以形成正确的错误景观。
这里有一个执行优化的过程，不需要插值（委托给solve_ivp）来交叉检查结果：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

def system(t, x, k0, k1, k2, p1, p2):
    return np.array([
        -k0 * x[0],
        p1 * x[0] - k1 * x[1],
        p2 * x[1] - k2 * x[2]
    ])

def solver(parameters, t=np.linspace(0, 1, 10), x0=np.ones(3)):
    solution = integrate.solve_ivp(system, [t.min(), t.max()], x0, args=parameters, t_eval=t)
    return solution.y

def residuals_factory(t, x):
    def wrapped(parameters):
        return 0.5 * np.sum(np.power(x - solver(parameters, t=t, x0=x[:, 0]), 2))
    return wrapped

# Synthetic dataset:
texp = np.linspace(0, 35, 15)
p0 = np.array([ 0.03693555,  0.38054633, -0.06252069,  1.41453107, -0.11159681])
x0 = np.array([0.24, 0.20, 0.89])
xexp = solver(p0, t=texp, x0=x0)

# Optimization
residuals = residuals_factory(texp, xexp)
solution = optimize.minimize(residuals, x0=[1, 1, 1, 1, 1])  # Initial guess

# Regression:
tlin = np.linspace(texp.min(), texp.max(), 200)
xhat = solver(solution.x, t=tlin, x0=x0)

fig, axe = plt.subplots()
for i in range(xexp.shape[0]):
    axe.scatter(texp, xexp[i, :])
    axe.plot(tlin, xhat[i, :])
axe.grid()

字符串
x1c 0d1x的数据
该方法既能对初始参数进行回归，又能保证拟合函数通过初始点，当维数小于10点时，收敛困难。
误差景观似乎有多个最小值，这可能需要调整初始参数以找到搜索到的最优值。还要注意，误差函数（残差）需要足够多的带有适当噪声的项，以便在最优值周围有足够的曲率和陡度来驱动梯度下降。

更新

我们可以修改回归过程以放松通过初始点的约束。新的签名看起来像：

def solver(parameters, t=np.linspace(0, 1, 10)):
    x0 = parameters[-3:]
    parameters = parameters[:-3]
    solution = integrate.solve_ivp(system, [t.min(), t.max()], x0, args=parameters, t_eval=t)
    return solution.y

def residuals_factory(t, x):
    def wrapped(parameters):
        return 0.5 * np.sum(np.power(x - solver(parameters, t=t), 2))
    return wrapped

型
那么这个问题有8个参数，因为初始条件也有影响。我们也可以将梯度下降限制为仅为正常数。

residuals = residuals_factory(texp, xexp)
solution = optimize.minimize(
    residuals, x0=[1, 1, 1, 1, 1, 1, 1, 1, 1],
    bounds=[(0, np.inf)]*8
)
# array([0.05157739, 0.32032917, 0.8680015 , 0.49081592, 0.90863933, 0.68856984, 0.17679386, 0.8890076 ])

型
下图显示了数据集的结果。

的
它没有太多的适应度，但它确认约束已经放松，并允许捕获比第一个签名更多的方差。
另一方面，拟合过程仍然适用于合成数据集（已知服从模型的数据集）上具有合理噪声的8个参数。

的

观察结果

您的数据显示出一些无法用您的模型拟合的振荡（例如：指数衰减无法拟合x0动态）;
如果振荡是由噪声引起的，那么它的幅度很高，并且与缺乏数据（非常少的点）相结合，这使得拟合困难
结果表明，只要数据点足够多，噪声合理，程序就能回归出参数

你有几个选择：

为每条曲线收集更多的点，以更好的精度（或拒绝已知为离群值的点）。这将允许您接受或拒绝当前模型;
如果你只有很少的数据点，那么改变你的ODE模型来捕捉动态，而不增加太多的复杂性（那么我们可能需要对生成这些数据的底层过程有更多的了解）;
重新定义您的问题，您是否正在寻找：
*regression：需要从模型中获取常数;
或者插值：你想通过点（无论底层模型如何）来预测其他地方的点？

更新2

带有RBF核的高斯过程可以通过你的系列，具有明显的适应性：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process.kernels import RBF
from sklearn.gaussian_process import GaussianProcessRegressor

texp = np.array([0, 5, 9, 18, 28, 38]).reshape(-1, 1)
xexp = np.array([
    [0.24, 0.71, 0.95, 0.26, 0.05, 0.22],
    [0.2, 0.62, 0.95, 0.51, 0.13, 0.05], 
    [0.89, 0.66, 0.95, 0.49, 0.28, 0.05]
])

tlin = np.linspace(0, 40, 200).reshape(-1, 1)
fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(9,3))

for i in range(3):

    kernel = 1 * RBF(length_scale=1.0)
    model = GaussianProcessRegressor(kernel=kernel, alpha=0.01**2)
    
    model.fit(texp, xexp[i, :])
    xpred, xstd = model.predict(tlin, return_std=True)

    axe = axes[i]
    axe.scatter(texp, xexp[i, :])
    axe.plot(tlin, xpred)
    axe.fill_between(tlin.ravel(), xpred - 1.96 * xstd, xpred + 1.96 * xstd, alpha=0.5)
    axe.grid()

型

的
在这种情况下，您可以使用精度提示在已知数据之间预测/插值点，但您没有要提取的回归参数。
您可以调整高斯过程alpha参数以平滑曲线（低于alpha=0.1**2）：
x1c4d 1x的

赞(0）回复(0）举报 11个月前

qncylg1j2#

这本身不是一个答案，但下面的代码使我更接近我想要的。现在我只需要构建一个自定义插值（like this, eg），这样我就可以更改它以最小化整体误差。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from lmfit import minimize, Parameters, report_fit
from scipy.interpolate import interp1d

# Measured data
t_measured = np.array([0, 5, 9, 18, 28, 38])
x0_measured = np.array([0.24, 0.71, 0.95, 0.26, 0.05, 0.22])
x1_measured = np.array([0.2, 0.62, 0.95, 0.51, 0.13, 0.05])
x2_measured = np.array([0.89, 0.66, 0.95, 0.49, 0.28, 0.05])

# Initial conditions
x00 = x0_measured[0]
x10 = x1_measured[0]
x20 = x2_measured[0]
y0 = [x00, x10, x20]

# Data dictionary
data = {'x0': x0_measured, 'x1': x1_measured, 'x2': x2_measured}

# Create an interpolation function for x0
x0_interp = interp1d(t_measured, x0_measured, kind='cubic', fill_value="extrapolate")

def f(y, t, paras):
    """
    System of differential equations using interpolated x0
    """
    x1 = y[0]
    x2 = y[1]

    # Use the interpolated x0 value
    x0 = x0_interp(t)

    try:
        k1 = paras['k1'].value
        k2 = paras['k2'].value
        p1 = paras['p1'].value
        p2 = paras['p2'].value
    except KeyError:
        k1, k2, p1, p2 = paras

    # Model equations, no equation for x0 as it's interpolated
    f1 = p1 * x0 - k1 * x1
    f2 = p2 * x1 - k2 * x2
    return [f1, f2]

def g(t, x10, x20, paras):
    """
    Solution to the ODE with initial conditions for x1 and x2
    """
    x = odeint(f, [x10, x20], t, args=(paras,))
    return x

def residual(paras, t, data):
    x10 = paras['x10'].value
    x20 = paras['x20'].value
    model = g(t, x10, x20, paras)

    x1_model = model[:, 0]
    x2_model = model[:, 1]

    # Compute residuals only for x1 and x2
    return np.concatenate([(x1_model - data['x1']).ravel(),
                           (x2_model - data['x2']).ravel()])

# Set parameters including bounds
params = Parameters()
params.add('x00', value=x00, vary=True)
params.add('x10', value=x10, vary=True)
params.add('x20', value=x20, vary=True)
params.add('k1', value=0.1, min=0.0001)
params.add('k2', value=0.1, min=0.0001)
params.add('p1', value=0.1, min=0.0001)
params.add('p2', value=0.1, min=0.0001)

# Fit model
result = minimize(residual, params, args=(t_measured, data), method='leastsq')

# Display fit report
report_fit(result)

# Generate fitted data
t_fine = np.linspace(0., 38., 100)
data_fitted = g(t_fine, x10, x20, result.params)
x0_fitted = x0_interp(t_fine)
x1_fitted = data_fitted[:, 0]  # x1 data from the ODE solution
x2_fitted = data_fitted[:, 1]  # x2 data from the ODE solution

# Plots
plt.figure(figsize=(12, 4))

# Plot for x0
plt.subplot(1, 3, 1)
plt.scatter(t_measured, x0_measured, marker='o', color='b', label='measured x0', s=75)
plt.plot(t_fine, x0_fitted, '-', linewidth=2, color='red', label='interpolated x0')
plt.legend()

# Plot for x1
plt.subplot(1, 3, 2)
plt.scatter(t_measured, x1_measured, marker='o', color='b', label='measured x1', s=75)
plt.plot(t_fine, x1_fitted, '-', linewidth=2, color='red', label='fitted x1')
plt.legend()

# Plot for x2
plt.subplot(1, 3, 3)
plt.scatter(t_measured, x2_measured, marker='o', color='b', label='measured x2', s=75)
plt.plot(t_fine, x2_fitted, '-', linewidth=2, color='red', label='fitted x2')
plt.legend()

plt.show()

字符串

的数据