Python小白的数学建模课-B5.新冠疫情SEIR模型

Python⼩⽩的数学建模课-B5.新冠疫情SEIR模型

1. SEIR 模型

1.1 SEIR 模型的提出

建⽴传染病的数学模型来描述传染病的传播过程，要根据传染病的发病机理和传播规律，结合疫情数据进⾏拟合分析，可以认识传染病的发展趋势，预测疫情持续时间和规模，分析和模拟各种防控措施对疫情发展的影响程度，为传染病防控⼯作提供决策指导，具有重要的理论意义和现实意义。

SI 模型是最简单的传染病传播模型，把⼈分为易感者（S 类）和患病者（I 类）两类，通过 SI 模型可以预测传染病⾼潮的到来；提⾼卫⽣⽔平、强化防控⼿段，降低病⼈的⽇接触率，可以推迟传染病⾼潮的到来。在 SI 模型基础上发展的 SIS 模型考虑患病者可以治愈⽽变成易感者，SIS 模型表⾯传染期接触数 σσ 是传染病传播和防控的关键指标，决定了疫情终将清零或演变为地⽅病长期存在。在 SI 模型基础上考虑病愈免疫的康复者（R 类）就得到 SIR 模型，通过 SIR 模型也揭⽰传染期接触数 σσ 是传染病传播的阈值，满

⾜ s0>1/σs0>1/σ 才会发⽣传染病蔓延，由此可以分析各种防控措施，如：提⾼卫⽣⽔平来降低⽇接触率λλ、提⾼医疗⽔平来提⾼⽇治愈率 μμ，通过预防接种达到体免疫来降低 s0s0 等。

传染病⼤多具有潜伏期（incubation period），也叫隐蔽期，是指从被病原体侵⼊肌体到最早临床症状出现的⼀段时间。在潜伏期的后期⼀般具有传染性。不同的传染病的潜伏期长短不同，从短⾄数⼩时到长达数年，但同⼀种传染病有固定的（平均）潜伏期。例如，流感的潜伏期为 1～3天，冠状病毒感染的潜伏期为4~7天，新型冠状病毒肺炎传染病（COVID-19）的潜伏期为1-14天（* 来⾃：新型冠状病毒肺炎诊疗⽅案试⾏第⼋版，潜伏时间 1～14天，多为3～7天，在潜伏期具有传染性），的潜伏期从数周到数⼗年。

SEIR 模型考虑存在易感者（Susceptible）、暴露者（Exposed）、患病者（Infectious）和康复者（Recovered）四类⼈，适⽤于具有潜伏期、治愈后获得终⾝免疫的传染病。易感者（S 类）被感染后成为潜伏者（E类），随后发病成为患病者（I 类），治愈后成为康复者（R类）。这种情况更为复杂，也更为接近实际情况。

SEIR 模型的仓室结构⽰意图如下：

1.2 SEIR 模型假设

1. 考察地区的总⼈数 N 不变，即不考虑⽣死或迁移；

2. ⼈分为易感者（S 类）、暴露者（E 类）、患病者（I 类）和康复者（R 类）四类；

3. 易感者（S 类）与患病者（I 类）有效接触即变为暴露者（E 类），暴露者（E 类）经过平均潜伏期后成为患病者（I 类）；患病者（I

类）可被治愈，治愈后变为康复者（R 类）；康复者（R类）获得终⾝免疫不再易感；

4. 将第 t 天时 S 类、E 类、I 类、R 类⼈的占⽐记为 s(t)s(t)、e(t)e(t)、i(t)i(t)、r(t)r(t)，数量分别为 S(t)S(t)、E(t)E(t)、I(t)I(t)、

R(t)R(t)；初始⽇期 t=0t=0 时，各类⼈占⽐的初值为 s0s0、e0e0、i0i0、r0r0；

5. ⽇接触数 λλ，每个患病者每天有效接触的易感者的平均⼈数；

6. ⽇发病率 δδ，每天发病成为患病者的暴露者占暴露者总数的⽐例；

7. ⽇治愈率 μμ，每天被治愈的患病者⼈数占患病者总数的⽐例，即平均治愈天数为 1/μ1/μ；

8. 传染期接触数 σ=λ/μσ=λ/μ，即每个患病者在整个传染期内有效接触的易感者⼈数。

1.3 SEIR 模型的微分⽅程

由

Ndsdt=−NλsiNdedt=Nλsi−NδeNdidt=Nδe−NμiNdrdt=Nμi(1)(2)(3)(4)

(1)Ndsdt=−Nλsi(2)Ndedt=Nλsi−Nδe(3)Ndidt=Nδe−Nμi(4)Ndrdt=Nμi

得：

SEIR 模型不能求出解析解，可以通过数值计算⽅法求解。

2. SEIR 模型的 Python 编程

霜叶红于二月花是什么植物2.1 Scipy ⼯具包求解微分⽅程组

SIS 模型是常微分⽅程初值问题，可以使⽤ Scipy ⼯具包的 scipy.integrate.odeint() 函数求数值解。

scipy.integrate.odeint(func, y0, t, args=())

**scipy.integrate.odeint() **是求解微分⽅程的具体⽅法，通过数值积分来求解常微分⽅程组。

odeint() 的主要参数：

func: callable(y, t, …) 　导数函数 f(y,t)f(y,t) ，即 y 在 t 处的导数，以函数的形式表⽰

y0: array：初始条件 y0y0，注意 SEIR模型是⼆元常微分⽅程组，初始条件为数组向量 y0=[i0,s0]y0=[i0,s0]

t: array：求解函数值对应的时间点的序列。序列的第⼀个元素是与初始条件 y0y0 对应的初始时间 t0t0；时间序列必须是单调递增或单调递减的，允许重复值。

args: 向导数函数 func 传递参数。当导数函数 f(y,t,p1,p2,..)f(y,t,p1,p2,..) 包括可变参数 p1,p2.. 时，通过 args =(p1,p2,..) 可以将参数p1,p2.. 传递给导数函数 func。

华大影院odeint() 的返回值：

y: array 　数组，形状为 (len(t),len(y0)，给出时间序列 t 中每个时刻的 y 值。

2.2 odeint() 求解 SEIR 模型的编程步骤

1. 导⼊ scipy、numpy、matplotlib 包。

2. 定义导数函数 f(y,t)f(y,t)。注意对于常微分⽅程（例如 SI模型）和常微分⽅程组（SEIR模型），y 分别表⽰标量和向量，函数定义略

有不同，以下给出两种情况的例程以供对⽐。

常微分⽅程的导数定义（SIS模型）

def dySIS(y, t, lamda, mu): # SIS 模型，导数函数

dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i

return dy_dt

常微分⽅程组的导数定义（SEIR模型）

def dySEIR(y, t, lamda, delta, mu): # SEIR 模型，导数函数

s, e, i = y

ds_dt = - lamda*s*i # ds/dt = -lamda*s*i

de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e

di_dt = delta*e - mu*i # di/dt = delta*e - mu*i

return np.array([ds_dt,de_dt,di_dt])

Python 可以直接对向量、向量函数进⾏定义和赋值，使程序更为简洁。但考虑读者主要是 Python ⼩⽩，⼜涉及到看着就⼼烦的微分⽅程组，所以我们宁愿把程序写得累赘⼀些，便于读者将程序与前⾯的微分⽅程组逐项对应。

1. 定义初值 y0y0 和 yy 的定义区间 [t0, t][t0, t]，注意初值为数组向量 y0=[s0,e0,i0]y0=[s0,e0,i0]。

2. 调⽤ odeint() 求 yy 在定义区间 [t0, t][t0, t] 的数值解。

SEIR 模型是三元常微分⽅程，返回值 y 是 len(t)*3 的⼆维数组。

2.3 Python例程：SEIR 模型

# modelCovid4_v1.py

# Demo01 of mathematical modeling for Covid2019

# SIR model for epidemic diseases.

# Crated：2021-06-13

# Python⼩⽩的数学建模课 @ Youcans

# 1. SEIR 模型，常微分⽅程组

from scipy.integrate import odeint # 导⼊ scipy.integrate 模块

import numpy as np # 导⼊ numpy包

import matplotlib.pyplot as plt # 导⼊ matplotlib包

def dySIS(y, t, lamda, mu): # SI/SIS 模型，导数函数

dy_dt = lamda*y*(1-y) - mu*y # di/dt = lamda*i*(1-i)-mu*i

return dy_dt

def dySIR(y, t, lamda, mu): # SIR 模型，导数函数

s, i = y # youcans

ds_dt = -lamda*s*i # ds/dt = -lamda*s*i

di_dt = lamda*s*i - mu*i # di/dt = lamda*s*i-mu*i

return np.array([ds_dt,di_dt])

def dySEIR(y, t, lamda, delta, mu): # SEIR 模型，导数函数

金属学报s, e, i = y # youcans

ds_dt = -lamda*s*i # ds/dt = -lamda*s*i

de_dt = lamda*s*i - delta*e # de/dt = lamda*s*i - delta*e

di_dt = delta*e - mu*i # di/dt = delta*e - mu*i

return np.array([ds_dt,de_dt,di_dt])

# 设置模型参数

number = 1e5 # 总⼈数

lamda = 0.3 # ⽇接触率, 患病者每天有效接触的易感者的平均⼈数

delta = 0.03 # ⽇发病率，每天发病成为患病者的潜伏者占潜伏者总数的⽐例

mu = 0.06 # ⽇治愈率, 每天治愈的患病者⼈数占患病者总数的⽐例

sigma = lamda / mu # 传染期接触数

fsig = 1-1/sigma

tEnd = 300 # 预测⽇期长度

t = np.arange(0.0,tEnd,1) # (start,stop,step)

i0 = 1e-3 # 患病者⽐例的初值

e0 = 1e-3 # 潜伏者⽐例的初值

s0 = 1-i0 # 易感者⽐例的初值

Y0 = (s0, e0, i0) # 微分⽅程组的初值

# odeint 数值解，求解微分⽅程初值问题

ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI 模型

ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS 模型

ySIR = odeint(dySIR, (s0,i0), t, args=(lamda,mu)) # SIR 模型

ySEIR = odeint(dySEIR, Y0, t, args=(lamda,delta,mu)) # SEIR 模型基础教育参考

# 输出绘图

print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig))

plt.title("Comparison among SI, SIS, SIR and SEIR models")

plt.xlabel('t-youcans')

plt.axis([0, tEnd, -0.1, 1.1])

plt.plot(t, ySI, 'cadetblue', label='i(t)-SI')

plt.plot(t, ySIS, 'steelblue', label='i(t)-SIS')

plt.plot(t, ySIR[:,1], 'cornflowerblue', label='i(t)-SIR')

# plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], 'cornflowerblue', label='r(t)-SIR')

plt.plot(t, ySEIR[:,0], '--', color='darkviolet', label='s(t)-SIR')

plt.plot(t, ySEIR[:,1], '-.', color='orchid', label='e(t)-SIR')

plt.plot(t, ySEIR[:,2], '-', color='m', label='i(t)-SIR')

plt.plot(t, 1-ySEIR[:,0]-ySEIR[:,1]-ySEIR[:,2], ':', color='palevioletred', label='r(t)-SIR') plt.legend(loc='right') # youcans

plt.legend(loc='right') # youcans

plt.show()

2.4 SI /SIS/SIR 模型与SEIR 模型的⽐较

例程 2.3 的参数和初值为：λ=0.3，δ=0.03，μ=0.06，(s0,e0,i0)=(0.001,0.001,0.998)λ=0.3，δ=0.03，μ=0.06，

(s0,e0,i0)=(0.001,0.001,0.998)，上图为例程的运⾏结果。

曲线 i(t)-SI 是 SI 模型的结果，患病者⽐例急剧增长到 1.0，所有⼈都被传染⽽变成患病者。

曲线 i(t)-SIS 是 SIS 模型的结果，患病者⽐例快速增长并收敛到某个常数，即稳态特征值 i∞=1−μ/λ=0.8i∞=1−μ/λ=0.8，表明疫情稳定，并将长期保持⼀定的患病率，称为地⽅病平衡点。

曲线 i(t)-SIR 是 SIR 模型的结果，患病者⽐例 i(t) 先上升达到峰值，然后再逐渐减⼩趋近于常数。

曲线 s(t)-SEIR、e(t)-SEIR、i(t)-SEIR、r(t)-SEIR 分别表⽰ SEIR模型中易感者（S类）、潜伏者（E类）、患病者（I类）和康复者（R 类）⼈的占⽐。

图中易感者⽐例 s(t) 单调递减并收敛到接近于 0 的稳定值。潜伏者⽐例 e(t) 曲线存在波峰，先逐渐上升⽽达到峰值，然后再逐渐减⼩，最终趋于 0。患病者⽐例 i(t) 曲线与潜伏者⽐例曲线类似，上升达到峰值后逐渐减⼩，最终趋于 0；但患病者⽐例曲线发展、达峰的时间⽐潜伏者曲线要晚⼀些，峰值强度也较低。康复者⽐例 r(i) 单调递增并收敛到⾮零的稳态值。以上分析只是对本图进⾏的讨论，并⾮普遍结论，取决于具体参数条件。

⽐较相同参数条件下 SIR 和 SEIR 模型的结果，SIR 模型中患病者⽐例 i(t) 的波形起点、峰值和终点到来的时间都显著早于 SEIR 模型，峰值强度也⾼于 SEIR 模型。这表明具有潜伏期的传染病，疫情发⽣和峰值的到来要晚于没有潜伏期的传染病，⽽且持续时间更长。

可贵的沉默教学设计3. SEIR 模型参数的影响

SEIR 模型中有⽇接触率 λλ 、⽇发病率 δδ 和⽇治愈率 μμ 三个参数，还有 i0、e0、s0i0、e0、s0 等初始条件，我们先⽤单因素分析的⽅法来观察参数条件对于疫情传播的影响。

3.1 初值条件 i0、e0、s0i0、e0、s0 初始条件的影响

SEIR 模型中有 i0、e0、s0i0、e0、s0 等 3个初始条件，组合众多⽆法穷尽。考虑实际情况中，疫情初始阶段尚⽆康复者，⽽潜伏者⽐例往往⾼于确诊的发病者，我们假定 e0/i0=2，r0=0e0/i0=2，r0=0，考察不同 i0i0 时的疫情传播情况。

通过对于该参数下不同的患病者、潜伏者⽐例初值条件的模拟，可以看到患病者、潜伏者⽐例的初值条件对疫情发⽣、达峰、结束的时间早晚具有直接影响，但对疫情曲线的形态和特征影响不⼤。不同初值条件下的疫情曲线，⼏乎是沿着时间指标平移的。

这说明如果不进⾏防控等⼈为⼲预，疫情传播过程与初始患病者、潜伏者⽐例关系并不⼤，该来的总会来。

图中患病率达到⾼峰后逐步降低，直⾄趋近于 0；易感率在疫情爆发后迅速下降，直⾄趋近于 0。但这⼀现象是基于具体的参数条件的观察，仅由该图并不能确定其是否普遍规律。

3.2 ⽇接触率λλ的影响

⾸先考察⽇接触率 λλ 的影响。

保持参数 δ=0.1，μ=0.06,(i0=0.001,e0=0.002,s0=0.997)δ=0.1，μ=0.06,(i0=0.001,e0=0.002,s0=0.997) 不变，λ=

[0.12,0.25,0.5,1.0,2.0]λ=[0.12,0.25,0.5,1.0,2.0] 时 i(t),s(t)i(t),s(t) 的变化曲线如下图所⽰。

血淋巴

通过对于该条件下⽇接触率的单因素分析，可以看到随着⽇接触率 λλ 的增⼤，患病率⽐例 i(t)i(t) 出现的峰值更早、更强，⽽易感者⽐例 s(t)s(t) 逐渐降低，但最终都趋于稳定。

3.3 ⽇发病率δδ的影响

本文发布于:2024-09-23 05:22:16，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/369867.html

上一篇：2020年秋季学期西南大学[0349]《数学建模》作业答案

下一篇：捕食者—被捕食者、竞争、共生三种模型的参数估计问题

标签：模型疫情患病者传染病

留言与评论（共有 0 条评论）