跨年龄组SEPIAR模型

#4.1 代码
#导入包
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import pandas as pd
from lmfit import Model, Parameters

# 定义带季节性参数的传播模型
def f(t, y, br, dr, ny, no, betayy_base, betaoy_base, betaoo_base, betayo_base, rho, omega1, omega2, omega3, gamma, gamma1, k1, k2):
    Sy, Ey, Ay, Py, Iy, Ry, So, Eo, Ao, Po, Io, Ro = y
    
    # 计算季节性调整的 beta 值，以12个月为一个周期，t - 20用于调节高峰位置，用sin函数模拟季节性波动，其中该函数“+1”是为保证该调整因子始终大于0
    betayy = betayy_base * (1 + np.sin((t - 20) * 2 * np.pi / 12)) 
    betaoy = betaoy_base * (1 + np.sin((t - 20) * 2 * np.pi / 12))
    betaoo = betaoo_base * (1 + np.sin((t - 20) * 2 * np.pi / 12))
    betayo = betayo_base * (1 + np.sin((t - 20) * 2 * np.pi / 12))

    
    # 计算传播动力学模型
    dSy_dt = ny * br - Sy * betayy * (Iy + k1 * Ay + k2 * Py) / ny - Sy * betaoy * (Io + k1 * Ao + k2 * Po) / ny - dr * Sy
    dEy_dt = Sy * betayy * (Iy + k1 * Ay + k2 * Py) / ny + Sy * betaoy * (Io + k1 * Ao + k2 * Po) / ny - rho * omega1 * Ey - (1 - rho) * omega2 * Ey - dr * Ey
    dAy_dt = rho * omega1 * Ey - gamma * Ay - dr * Ay
    dPy_dt = (1 - rho) * omega2 * Ey - omega3 * Py - dr * Py
    dIy_dt = omega3 * Py - gamma1 * Iy - dr * Iy
    dRy_dt = gamma * Ay + gamma1 * Iy - dr * Ry
    dSo_dt = So * betaoo * (Io + k1 * Ao + k2 * Po) / no - So * betayo * (Iy + k1 * Ay + k2 * Py) / no - dr * So
    dEo_dt = So * betaoo * (Io + k1 * Ao + k2 * Po) / no + So * betayo * (Iy + k1 * Ay + k2 * Py) / no - rho * omega1 * Eo - (1 - rho) * omega2 * Eo - dr * Eo
    dAo_dt = rho * omega1 * Eo - gamma * Ao - dr * Ao
    dPo_dt = (1 - rho) * omega2 * Eo - omega3 * Po - dr * Po
    dIo_dt = omega3 * Po - gamma1 * Io - dr * Io
    dRo_dt = gamma * Ao + gamma1 * Io - dr * Ro
    
    # 返回微分方程的结果
    return [dSy_dt, dEy_dt, dAy_dt, dPy_dt, dIy_dt, dRy_dt, dSo_dt, dEo_dt, dAo_dt, dPo_dt, dIo_dt, dRo_dt]


# 拟合函数，计算季节性参数影响后的结果
def fit_func(t, betayy_base, betaoy_base, betaoo_base, betayo_base):
    # 初始条件假设，可根据实际情况调整
    Ey0 = 100
    Ay0 = 5
    Py0 = 10
    Ry0 = 0  
    Eo0 = 20
    Ao0 = 1
    Po0 = 2
    Ro0 = 0  
    Sy0 = ny - Ey0 - Iy0 - Ay0 - Py0 - Ry0
    So0 = no - Eo0 - Io0 - Ao0 - Po0 - Ro0
    y0 = [Sy0, Ey0, Ay0, Py0, Iy0, Ry0, So0, Eo0, Ao0, Po0, Io0, Ro0]
    
    # 求解微分方程组，使用四阶龙格-库塔法
    result = solve_ivp(f, [t[0], t[-1]], y0, args=(
        br, dr, ny, no, betayy_base, betaoy_base, betaoo_base, betayo_base, rho, omega1, omega2, omega3, gamma, gamma1, k1, k2), t_eval=t)
    return [result.y[4], result.y[10]] 


# 读取数据
data = pd.read_csv('./SEPIAR_data.csv', encoding='gbk')
t_data = np.arange(0, 48) 
Iy_data = data.iloc[:, 1]
Io_data = data.iloc[:, 2]
Iy0 = Iy_data[0] 
Io0 = Io_data[0]

# 模型参数设置
ny = 3350788 #0-14岁人口数
no = 23918526 #14岁以上人口数
br = 0.00789/12 #出生率 
dr = 0.00532/12  #死亡率
k1 = 0.3 #无症状个体与有症状个体的相对传播率
k2 = 0.3  #无症状前个体与有症状个体的相对传播率
rho = 0.2 #无症状个体比例
omega1 = 1/3  #E至A的速率
omega2 = 1/3  #E至P的速率
omega3 = 1   #P至I的速率
gamma = 1/2  #A至R的速率
gamma1 = 1/2  #I至R的速率

params = Parameters()
params.add('betayy_base',min = 0.55, max = 1)
params.add('betayo_base', min = 0.1, max = 0.5)
params.add('betaoo_base', min = 0.05, max = 0.1)
params.add('betaoy_base', min = 0.1, max = 0.5)

#拟合
fit_model = Model(fit_func)
result1 = fit_model.fit([Iy_data, Io_data], params, t=t_data)

#输出拟合结果
result1

#绘图
plt.figure(figsize=(10, 6))
plt.plot(t_data, Iy_data, label='Iy Data')  
plt.plot(t_data, Io_data, label='Io Data')  
plt.plot(t_data, result1.best_fit[0], 'r-', label='Best Iy Fit')  
plt.plot(t_data, result1.best_fit[1], label='Best Io Fit')  
plt.xlabel('Time (days)')
plt.ylabel('Infected')
plt.title('Model Fit to Data')
plt.legend()
plt.grid(True)
# plt.savefig("C:/Users/86187/Desktop/图.svg")

plt.show()

# 绘制舱室的变化曲线
plt.figure(figsize=(10, 6)) 
# 使用拟合函数计算所有状态变量
all_states = fit_func(t_data, result1.params['betayy_base'].value, result1.params['betayo_base'].value, result1.params['betaoo_base'].value, result1.params['betaoy_base'].value)  
Ey0 = 80
Ay0 = 10
Py0 = 10
Ry0 = 0  
Eo0 = 10
Ao0 = 2
Po0 = 2
Ro0 = 0  
Sy0 = ny - Ey0 - Iy0 - Ay0 - Py0 - Ry0
So0 = no - Eo0 - Io0 - Ao0 - Po0 - Ro0
y0 = [Sy0, Ey0, Ay0, Py0, Iy0, Ry0, So0, Eo0, Ao0, Po0, Io0, Ro0]

# 求解 ODE 以获取所有状态变量
result = solve_ivp(f, [t_data[0], t_data[-1]], y0, args=(br, dr, ny, no, result1.params['betayy_base'].value, result1.params['betayo_base'].value, result1.params['betaoo_base'].value, result1.params['betaoy_base'].value, rho, omega1, omega2, omega3, gamma, gamma1, k1, k2), t_eval=t_data)
Sy_fit, Ey_fit, Py_fit, Iy_fit, Ay_fit, Ry_fit, So_fit, Eo_fit, Po_fit, Io_fit, Ao_fit, Ro_fit = result.y  

#plt.plot(t_data, Sy_fit, label='Sy')
plt.plot(t_data, Ey_fit, label='Ey')
plt.plot(t_data, Py_fit, label='Py')
plt.plot(t_data, Iy_fit, label='Iy')
plt.plot(t_data, Ay_fit, label='Ay')
plt.plot(t_data, Ry_fit, label='Ry')
#plt.plot(t_data, So_fit, label='So')
plt.plot(t_data, Eo_fit, label='Eo')
plt.plot(t_data, Po_fit, label='Po')
plt.plot(t_data, Io_fit, label='Io')
plt.plot(t_data, Ao_fit, label='Ao')
plt.plot(t_data, Ro_fit, label='Ro')
plt.xlabel('Time(days)')
plt.ylabel('Individuals')
plt.legend()
plt.title('SEPIAR Model Dynamics')
plt.tight_layout()
#plt.savefig("SEPIAR dynamic")
plt.show()