使用Scipy Minimize估计模型参数：非时序观测数据适配问题

我想用Scipy的minimize方法估计模型中的beta和gamma参数，观测要求是模型达到平衡态时：

• 患病率I的平衡值为0.4
• 发病率J_diff的平衡值为0.3

但原代码执行minimize时直接返回初始设定的x0值，无法完成参数估计，需要修正代码。

• 参数被硬编码覆盖：在peak_infections函数中，刚从输入x取出beta和gamma，立刻又被赋值为15, 2/5，导致优化过程中x的参数根本没被用到，残差始终不变，minimize直接返回初始值。
• 仅考虑I的残差：没有把J_diff（对应代码中的cInc）的平衡目标纳入残差计算。
• 模拟时间不足：只模拟到20年，可能系统还没达到平衡态，无法获取真实的平衡值。

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

def get_equilibrium_vals(x):
    # Total population, N.
    N = 1
    # Initial conditions
    I0, R0 = 0.001, 0
    U0 = N - I0 - R0
    Lf0, Ls0, J0 = 0, 0, I0
    
    # 用输入的x参数，不再硬编码覆盖
    beta = x[0]
    gamma = x[1]
    mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
    u, v, w = 0.083, 0.88, 0.0006

    # 延长模拟时间，确保系统达到平衡（模拟到1000年，取最后一个时间点的值）
    times = np.arange(0, 1001, 10)

    def deriv(y, t, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
        U, Lf, Ls, I, R, cInc = y
        b = (mu * (U + Lf + Ls + R)) + (muTB * I)
        lamda = beta * I
        clamda = 0.2 * lamda
        dU = b - ((lamda + mu) * U)
        dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
        dLs = (u * Lf) - ((w + clamda + mu) * Ls)
        dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
        dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
        cI = w*Ls + v*Lf + (rho * R)
        return dU, dLf, dLs, dI, dR, cI

    # 求解ODE
    solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), times, 
                   args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
    U, Lf, Ls, I, R, cInc = solve.T
    
    # 返回最后一个时间点的平衡态值（确保系统已稳定）
    return I[-1], cInc[-1]

def residual(x):
    # 目标平衡值
    target_I = 0.4
    target_J_diff = 0.3
    
    # 获取当前参数下的平衡态值
    sim_I, sim_J_diff = get_equilibrium_vals(x)
    
    # 计算两个目标的残差平方和
    return (sim_I - target_I)**2 + (sim_J_diff - target_J_diff)**2

# 初始参数猜测
x0 = [12, 0.4] 
# 执行优化
res = minimize(residual, x0, method="Nelder-Mead", options={'fatol':1e-06})
print("估计的beta和gamma参数：", res.x)

• 移除参数硬编码：删除了原代码中beta, gamma = 15, 2/5的赋值，确保使用输入的x参数进行优化。
• 纳入双目标残差：残差函数同时计算I和J_diff与目标值的平方和，满足两个平衡态要求。
• 延长模拟时间：将模拟时间延长到1000年，取最后一个时间点的值作为平衡态估计，避免因模拟时间不足导致的非平衡值干扰。
• 优化函数重命名：将peak_infections改为get_equilibrium_vals，更贴合功能（获取平衡态值而非峰值）。

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