如何用插值拟合 ODE 系统？

Question

给定数据

x0

、

x1

和

x2

，假设我想在以下 ODE 系统中拟合参数

k0

、

k1

、

k2

、

p1

和

p2

以下代码实现了我想要的功能

pip install lmfit

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

# Function definitions
def f(y, t, paras):
    """
    System of differential equations
    """

    x0 = y[0]
    x1 = y[1]
    x2 = y[2]

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

    # Updated model equations
    f0 = -k0 * x0
    f1 = p1 * x0 - k1 * x1
    f2 = p2 * x1 - k2 * x2
    return [f0, f1, f2]

def g(t, x000, paras):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x000
    """
    x = odeint(f, x000, t, args=(paras,))
    return x

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

    # Assume you have data for x1 and x3 as well
    x0_model = model[:, 0]
    x1_model = model[:, 1]
    x2_model = model[:, 2]

    # Compute residuals for all variables
    return np.concatenate([(x0_model - data['x0']).ravel(),
                           (x1_model - data['x1']).ravel(),
                           (x2_model - data['x2']).ravel()])


# 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}

# Set parameters including bounds
params = Parameters()
params.add('x00', value=x00, vary=False)
params.add('x10', value=x10, vary=False)
params.add('x20', value=x20, vary=False)
params.add('k0', value=0.2, min=0.0001, max=2.)
params.add('k1', value=0.3, min=0.0001, max=2.)
params.add('k2', value=0.1, min=0.0001, max=2.)  # Assuming k2 is needed
params.add('p1', value=0.1, min=0.0001, max=2.)  # Assuming p1 is needed
params.add('p2', value=0.1, min=0.0001, max=2.)  # Assuming p2 is needed

# Fit model
result = minimize(residual, params, args=(t_measured, data), method='leastsq')  # leastsq nelder
# Check results of the fit
data_fitted = g(np.linspace(0., 38., 100), y0, result.params)

# Updated plotting to include x0, x1, and x2
plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
plt.scatter(t_measured, x0_measured, marker='o', color='b', label='measured x0', s=75)
plt.plot(np.linspace(0., 38., 100), data_fitted[:, 0], '-', linewidth=2, color='red', label='fitted x0')
plt.legend()

plt.subplot(1, 3, 2)
plt.scatter(t_measured, x1_measured, marker='o', color='b', label='measured x1', s=75)
plt.plot(np.linspace(0., 38., 100), data_fitted[:, 1], '-', linewidth=2, color='red', label='fitted x1')
plt.legend()

plt.subplot(1, 3, 3)
plt.scatter(t_measured, x2_measured, marker='o', color='b', label='measured x2', s=75)
plt.plot(np.linspace(0., 38., 100), data_fitted[:, 2], '-', linewidth=2, color='red', label='fitted x2')
plt.legend()

plt.show()

现在，我不想用 ODE 来拟合

x0

，而是想对其进行插值，这样我就有更多的自由。理想情况下，我想定义一个（或多个）插值参数，该参数也将进入优化问题。关于最好的方法有什么想法吗？

我的想法：类似于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()

Answer 1

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

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

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

def f(y, t, paras):
    """
    Modified 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):
    """
    Modified 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()])

# 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}

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