基于两组观测数据求解常微分方程中参数wf的方法咨询

可行性说明 #

完全可以通过两组观测数据反求参数wf。你仅需要求解1个未知参数，却有2个观测点约束，求解过程是适定的，不存在欠定无法收敛的问题。

实现思路 #

• 将原模型改造为接收wf作为入参、返回指定时刻（t=0、t=tr）的C(t)预测值的函数
• 定义损失函数：计算预测值与实际观测值的残差平方和
• 调用数值优化工具，在wf的合理取值范围（0 < wf < 1）内最小化损失函数，得到最优wf取值

完整实现代码 #

from scipy.integrate import odeint
import numpy as np
from scipy.optimize import minimize

# 固定参数（不需要优化的参数）
S = 0.0001
M = 30.03
K = 113.6561
Vr = 58
R = 8.3145
T = 298.15
Q = 0.000133
Vp = 0.000022
Mr = 36
Pvap = 1400
tr = 1200
mass = 40000

# ---------------------- 此处替换为你的实际观测值 ----------------------
obs_C0 = 0  # t=0时刻观测到的C取值
obs_Ctr = 0.057  # t=tr时刻观测到的C取值，此处示例值为wf=0.001时的模拟结果
# -------------------------------------------------------------------

# 定义ODE模型，将wf作为外部参数传入
def model(x, t, wf):
    C = x[0]  # Cair(t)
    c = x[1]  # Cv(t)
    a = Q + (K * S / Vr)
    b = (K * S * M) / (Vr * R * T)
    s = (K * S * M) / (Vp * R * T)
    w = (1 - wf) * 1000
    Peq = (c * Pvap) / (c + w * c * M / Mr)
    Pair = (C * R * T) / M
    dcdt = -s * (Peq - Pair)
    if t <= tr:
        dCdt = -a * C + b * Peq
    else:
        dCdt = -a * C
    return [dCdt, dcdt]

# 定义预测函数：给定wf，返回t=0和t=tr时刻的C值
def predict_C(wf):
    # 计算初始状态
    Cv0 = (mass / Vp) * wf
    Cair0 = obs_C0  # 用t=0的观测值作为初始条件
    # 求解0到tr时刻的ODE
    t_eval = [0, tr]
    x = odeint(model, [Cair0, Cv0], t_eval, args=(wf,))
    return x[:, 0]  # 返回两个时刻的C值

# 定义损失函数：计算预测值与观测值的残差平方和
def loss_func(wf):
    # wf是数组格式，取标量
    wf_val = wf[0]
    # 限制wf范围在0~1之间，超出范围给大损失
    if wf_val <= 0 or wf_val >=1:
        return 1e10
    pred_C0, pred_Ctr = predict_C(wf_val)
    # 计算残差平方和
    loss = (pred_C0 - obs_C0)**2 + (pred_Ctr - obs_Ctr)**2
    return loss

# 优化求解
init_guess = [0.001]  # 初始猜测值，可以根据常识调整
bounds = [(1e-6, 0.999)]  # wf的合理范围
result = minimize(loss_func, x0=init_guess, bounds=bounds, method='L-BFGS-B')

# 输出结果
if result.success:
    optimized_wf = result.x[0]
    print(f"反求得到的最优wf值为：{optimized_wf:.6f}")
    print(f"拟合损失：{result.fun:.8f}")
else:
    print("优化失败，原因：", result.message)

注意事项 #

• 代码中的obs_C0和obs_Ctr需要替换为你实际的观测数据
• 如果观测的C(0)和原代码初始假设的0不符，代码中已经自动用观测值作为初始条件，无需额外调整
• 如果优化结果误差过大，可以尝试调整初始猜测值init_guess，或者更换优化方法（比如将method参数改为Nelder-Mead）

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