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

传染病的数学模型是数学建模中的典型问题，常见的传染病模型有 SI、SIR、SIRS、SEIR 模型。

SIR 模型将人群分为易感者（S类）、患病者（I类）和康复者（R 类），考虑了患病者治愈后的免疫能力。

本文详细给出了 SIR 模型微分方程、相空间分析的建模、例程、结果和分析，让小白都能懂。

『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。

欢迎关注『Python小白的数学建模课 @ Youcans』系列，每周持续更新
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python小白的数学建模课-10.微分方程边值问题

1. 疫情传播 SIR 模型

传染病的传播特性不可能通过真实的试验开展研究，因此需要针对不同的传染病传播方式和流行特点建立相应的数学模型，并对模型进行理论研究和数值模拟。通过研究发现传染病传播的特征阈值，就可以为预防和控制传染病提供数据支撑和防控策略。

1927年，W. Kermack 在论文 “Contributions to the mathcmatical theory of epidemics” 中研究了伦敦黑死病和孟买瘟疫的流行过程，创造性地提出了 SIR 模型。

SI 模型和 SIS 模型将人群分为感染者（S类）和患病者（I类）两类人群，但大多数传染病的患者在治愈后就有很强的免疫力，终身或在一段时期内不再会被感染而变成病人。这类人群称为病愈免疫的康复者（R 类）。康复者已经退出传染系统，对于致死性疾病的死亡者也可以用该类别描述其传播特性。

SIR 模型适用于具有易感者、患病者和康复者三类人群，可以治愈，且治愈后终身免疫不再复发的疾病，例如天花、肝炎、麻疹等免疫力很强的传染病。

在这里插入图片描述

SIR 模型假设：

考察地区的总人数 N 不变，即不考虑生死或迁移；
人群分为易感者（S类）、患病者（I类）和康复者（R 类）三类；
易感者（S类）与患病者（I类）有效接触即被感染，变为患病者（I类）；患病者（I类）可被治愈，治愈后变为康复者；康复者（R类）获得终身免疫不再易感；无潜伏期；
将第 t 天时 S类、I 类、R 类人群的占比记为 $s (t)$ 、 $i (t)$ 、 $r (t)$ ，数量为 $S (t)$ 、 $I (t)$ 、 $R (t)$ ；初始日期 $t = 0$ 时， S类、I 类、R 类人群占比的初值为 $s_0$ 、 $i_0$ 、 $r_0$ 。
日接触数 $λ\lambda$ ，每个患病者每天有效接触的易感者的平均人数；
日治愈率 $μ\mu$ ，每天被治愈的患病者人数占患病者总数的比例，即平均治愈天数为 $1/μ1/\mu$ ；
传染期接触数 $σ=λ/μ\sigma = \lambda / \mu$ ，即每个患病者在整个传染期内有效接触的易感者人数。

SIR 模型的微分方程：
由
${Ndsdt=−NλsiNdidt=Nλsi−μis(t)+i(t)+r(t)=1\begin{cases} & N \frac{ds}{dt} = -N \lambda s i\\ & N\frac{di}{dt} = N\lambda s i - \mu i\\ & s(t) + i(t) + r(t) = 1 \end{cases}$

得：
${dsdt=−λsi,s(0)=s0didt=λsi−μi,i(0)=i0\begin{cases} & \frac{ds}{dt} = -\lambda s i , &s(0)=s_0\\ & \frac{di}{dt} = \lambda s i - \mu i , &i(0)=i_0\\ \end{cases}$

SIR 模型不能求出解析解，只能通过数值计算方法求解。

2. SIR 模型的 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)$ ，即 y 在 t 处的导数，以函数的形式表示
y0: array：　　初始条件 $y_0$ ，注意 SIR模型是二元常微分方程组，初始条件为数组向量 $y_0=[i_0, s_0]$
t: array：　　求解函数值对应的时间点的序列。序列的第一个元素是与初始条件 $y_0$ 对应的初始时间 $t_0$ ；时间序列必须是单调递增或单调递减的，允许重复值。
args: 向导数函数 func 传递参数。当导数函数 $f (y, t, p 1, p 2, . .)$ 包括可变参数 p1,p2… 时，通过 args =(p1,p2,…) 可以将参数p1,p2… 传递给导数函数 func。

odeint() 的返回值：

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

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

导入 scipy、numpy、matplotlib 包。
定义导数函数 $f (y, t)$ 。注意对于常微分方程（例如 SIS模型）和常微分方程组（SIR模型），y 分别表示标量和向量，函数定义略有不同，以下给出两种情况的例程以供对比。

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

def dySIS(y, t, lamda, mu):  # SIS 模型，导数函数dy_dt = lamda*y*(1-y) - mu*y  # di/dt = lamda*i*(1-i)-mu*ireturn dy_dt

常微分方程组的导数定义（SIR模型）

def dySIR(y, t, lamda, mu):  # SIR 模型，Y=[i,s] 点的导数dy/dti, s = ydi_dt = lamda*s*i - mu*i  # di/dt = lamda*s*i-mu*ids_dt = -lamda*s*i  # ds/dt = -lamda*s*ireturn np.array([di_dt,ds_dt])

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

定义初值 $y_0$ 和 $y$ 的定义区间 $t_0,\ t]$ ，注意初值为数组向量 $y_0=[i_0, s_0]$ 。
调用 odeint() 求 $y$ 在定义区间 $t_0,\ t]$ 的数值解。

SIR 模型是二元常微分方程，返回值 y 是 len(t)*2 的二维数组。

2.3 Python例程：SI 模型、SIS 模型与SIR 模型的比较

# modelCovid3_v1.py
# Demo01 of mathematical modeling for Covid2019
# SIR model for epidemic diseases.
# Copyright 2021 Youcans, XUPT
# Crated：2021-06-12
# Python小白的数学建模课 @ Youcans# 1. SIR 模型，常微分方程组
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*ireturn dy_dtdef dySIR(y, t, lamda, mu):  # SIR 模型，导数函数i, s = ydi_dt = lamda*s*i - mu*i  # di/dt = lamda*s*i-mu*ids_dt = -lamda*s*i  # ds/dt = -lamda*s*ireturn np.array([di_dt,ds_dt])# 设置模型参数
number = 1e5  # 总人数
lamda = 0.2  # 日接触率, 患病者每天有效接触的易感者的平均人数
sigma = 2.5  # 传染期接触数
mu = lamda/sigma  # 日治愈率, 每天被治愈的患病者人数占患病者总数的比例
fsig = 1-1/sigma
tEnd = 200  # 预测日期长度
t = np.arange(0.0,tEnd,1)  # (start,stop,step)
i0 = 1e-4  # 患病者比例的初值
s0 = 1-i0  # 易感者比例的初值
Y0 = (i0, s0)  # 微分方程组的初值print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda,mu,sigma,fsig))# odeint 数值解，求解微分方程初值问题
ySI = odeint(dySIS, i0, t, args=(lamda,0))  # SI 模型
ySIS = odeint(dySIS, i0, t, args=(lamda,mu))  # SIS 模型
ySIR = odeint(dySIR, Y0, t, args=(lamda,mu))  # SIR 模型# 绘图
plt.title("Comparison among SI, SIS and SIR models")
plt.xlabel('t-youcans')
plt.axis([0, tEnd, -0.1, 1.1])
plt.axhline(y=0,ls="--",c='c')  # 添加水平直线
plt.plot(t, ySI, ':g', label='i(t)-SI')
plt.plot(t, ySIS, '--g', label='i(t)-SIS')
plt.plot(t, ySIR[:,0], '-r', label='i(t)-SIR')
plt.plot(t, ySIR[:,1], '-b', label='s(t)-SIR')
plt.plot(t, 1-ySIR[:,0]-ySIR[:,1], '-m', label='r(t)-SIR')
plt.legend(loc='best')  # youcans
plt.show()

2.4 SI 模型、SIS 模型与SIR 模型的比较

在这里插入图片描述

本图为例程 2.3 的运行结果，参数和初值为： $λ=0.2，μ=0.08，(i0,s0,r0)=(0.0001,0.9999,0)\lambda =0.2，\mu=0.08，(i_0,s_0,r_0)=(0.0001,0.9999,0)$ 。

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

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

曲线 i(t)-SIR、s(t)-SIR、r(t)-SIR 分别是 SIR 模型的易感者（S类）、患病者（I 类）、康复者（R 类）人群的占比。图中易感者比例 s(t) 单调递减并收敛到非零的稳态值 $s∞s_\infty$ ，康复者比例 r(i) 单调递增并收敛到非零的稳态值 $r∞r_\infty$ ，患病者比例 i(t) 先上升达到峰值，然后再逐渐减小趋近于 0 。

3. SIR 模型参数的影响

SIR 模型中有日接触率 $λ\lambda$ 与日治愈率 $μ\mu$ 两个参数，还有 $i_0、s_0$ 两个初始条件，共有 4 个可以调整的参数条件都会影响微分方程的解，也就是会影响患病者、易感者比例的时间变化曲线。其中的各种组合无穷无尽，如果没有恰当的研究方法、不能把握内在的规律，即使在几十、几百组参数条件下进行模拟，仍然只是盲人摸象、管中窥豹。

3.1 初值条件 $i_0, s_0$ 的影响

下面考察初值条件 $i_0, s_0$ 的影响。固定参数 $λ=0.2,μ=0.02\lambda=0.2, \mu=0.02$ 不变，不同初值条件下 $i (t), s (t)$ 的变化曲线如下图所示。
在这里插入图片描述

通过对于该参数条件下不同初值条件的单因素分析，可以看到患病者比例、易感者比例的初值条件对疫情发生、达峰、结束的时间早晚具有直接影响，但对疫情曲线的形态和特征影响不大。不同初值条件下的疫情曲线，几乎是沿着时间指标平移的。

这说明如果不进行治疗防控等人为干预，疫情传播过程与初始患病率无关，该来的总会来。

图中患病率达到高峰后逐步降低，直至趋近于 0；易感率在疫情爆发后迅速下降，直至趋近于 0。但这一现象是基于具体的参数条件 $λ=0.2,μ=0.02\lambda=0.2, \mu=0.02$ 的观察，无法确定其是否普遍规律。

3.2 日接触率 $λ\lambda$ 的影响

首先考察日接触率 $λ\lambda$ 的影响。固定参数 $μ=0.25,i0=0.002,s0=1−i0\mu=0.25, i_0=0.002,s_0=1-i_0$ 不变，$\lambda = [0.2, 0.25, 0.5, 1.0, 2.0] $ 时 $i (t), s (t)$ 的变化曲线如下图所示。
在这里插入图片描述

通过对于该条件下日接触率的单因素分析，可以看到随着日接触率 $λ\lambda$ 的增大，患病率比例 $i (t)$ 出现的峰值更早、更强，而易感者比例 $s (t)$ 从几乎不变到迅速降低，但最终都趋于稳定。

对本图我们好像感觉到存在一些规律，但又似乎说不清，很难总结出来。即便总结出某些特征，也只能说是在该固定参数条件下的特征，不能说是 SIR 模型的共有特征。即便反复地改变固定参数的取值或日接触率的范围，情况也差不多。

3.3 日治愈率 $μ\mu$ 的影响

下面考察日治愈率 $μ\mu$ 的影响。固定参数 $λ=0.2,i0=0.002,s0=1−i0\lambda=0.2, i_0=0.002,s_0=1-i_0$ 不变， $μ=[0.4,0.2,0.1,0.05,0.025]\mu = [0.4, 0.2, 0.1, 0.05, 0.025]$ 时 $i (t), s (t)$ 的变化曲线如下图所示。
在这里插入图片描述

通过对于该条件下日治愈率的单因素分析，可以看到随着日治愈率 $μ\mu$ 的减小，患病率比例 $i (t)$ 出现的峰值更强、也稍早，而易感者比例 $s (t)$ 从几乎不变到迅速降低，但最终都趋于稳定。

对于本图的观察和分析情况与上图是差不多的，看起来内容更丰富，似乎也有规律可循，但很难说的清，只能做一些简单的描述。即便进行更多的模拟，情况也差不多。

这是因为，对于SIR 模型这类微分方程，对结果具有决定性影响的特征参数，往往不是模型中的某个参数，而是多个参数特定关系的组合，因此仅从单因素实验很难充分反映模型中的内在特征。

不过，对于数学建模，通过几组单因素实验获得一系列曲线、图表，再从各个角度对结果进行一些描述和解读，就已经足够了。

4. SIR 模型的相空间分析

4.1 SIR 模型的相轨迹方程

SIR 模型不能求出解析解，可以通过相空间方法来研究解的周期性、稳定性。

由于患病者比例 $i (t)$ 和易感者比例 $s (t)$ 都是时间 t 的函数，因此当 t 取任意值时都对应着 $i - s$ 平面上的一个点，当 t 连续变化时对应着 $i - s$ 平面上的一条轨迹，称为相轨迹。通过相轨迹图可以分析微分方程的性质。

对于 SIR 模型，消去 dt 可以得到：
$dids=1σs−1,i(s=s0)=i0\frac{di}{ds} = \frac{1}{\sigma s} - 1 ,\; i(s=s_0)=i_0$
该微分方程的解为：
$(s_0 + i_0) - s +\frac{1}{\sigma} ln \frac{s}{s_0}$

4.2 Python例程：SIR 模型的相轨迹

# modelCovid3_v1.py
# Demo01 of mathematical modeling for Covid2019
# SIR model for epidemic diseases.
# Copyright 2021 Youcans, XUPT
# Crated：2021-06-12
# Python小白的数学建模课 @ Youcans# 2. SIR 模型，常微分方程组 相空间分析
from scipy.integrate import odeint  # 导入 scipy.integrate 模块
import numpy as np  # 导入 numpy包
import matplotlib.pyplot as plt  # 导入 matplotlib包def dySIR(y, t, lamda, mu):  # SIR 模型，导数函数i, s = ydi_dt = lamda*s*i - mu*i  # di/dt = lamda*s*i-mu*ids_dt = -lamda*s*i  # ds/dt = -lamda*s*ireturn np.array([di_dt,ds_dt])# 设置模型参数
number = 1e5  # 总人数
lamda = 0.2  # 日接触率, 患病者每天有效接触的易感者的平均人数
sigma = 2.5  # 传染期接触数
mu = lamda/sigma  # 日治愈率, 每天被治愈的患病者人数占患病者总数的比例
fsig = 1-1/sigma
print("lamda={}\tmu={}\tsigma={}\t(1-1/sig)={}".format(lamda, mu, sigma, fsig))# odeint 数值解，求解微分方程初值问题
tEnd = 200  # 预测日期长度
t = np.arange(0.0,tEnd,1)  # (start,stop,step)
s0List = np.arange(0.01,0.91,0.1)  # (start,stop,step)
for s0 in s0List:  # s0, 易感者比例的初值i0 = 1 - s0    # i0, 患病者比例的初值Y0 = (i0, s0)  # 微分方程组的初值ySIR = odeint(dySIR, Y0, t, args=(lamda,mu))  # SIR 模型plt.plot(ySIR[:,1], ySIR[:,0])# 绘图
plt.title("Phase trajectory of SIR models")
plt.axis([0, 1, 0, 1])
plt.plot([0,1],[1,0],'b-')
plt.plot([1/sigma,1/sigma],[0,1-1/sigma],'b--')
plt.xlabel('s(t)-youcans')
plt.ylabel('i(t)-xupt')
plt.text(0.8,0.9,r"$1/\sigma$ = {}".format(1/sigma),color='b')
plt.show()

4.3 SIR 模型的相轨迹分析

上图为例程 3.2 的运行结果（ $λ=0.2，μ=0.08，1/σ=0.4\lambda =0.2，\mu=0.08，1/\sigma=0.4$ ），是 SIR 模型的相轨迹图。

图中每一条 i-s 曲线，从直线 i(t)+s(t)=1 上的某一初值点出发，最终收敛于 s轴上的某一点，对应着某一个初值条件下的患病者与易感者比例随时间的变化关系。

利用相轨迹图可以分析和讨论 SIR 模型的性质：

任一条 i-s 曲线都收敛于s轴上的一点，即 $i∞=0i_\infty=0$ ，表明不论初始条件如何，患病者终将清零。
患病者比例在 $s=1/σs=1/\sigma$ 时达到峰值。若易感者比例的初值 $s0>1/σs_0>1/\sigma$ ，患病者比例先增长，在 $s=1/σs=1/\sigma$ 时达到峰值 $imax=(s0+i0)−(1+ln(σs0))/σi_{max} = (s_0 + i_0) - (1+ln (\sigma s_0))/\sigma$ ，然后下降，终将清零；若易感者比例的初值 $s0<1/σs_0<1/\sigma$ ，患病者比例单调递减，终将清零。
易感者比例单调递减，易感者的最终比例是相轨迹与 s轴在 (0,1/sigma) 内交点的横坐标。易感者最终比例虽然与初值有关，但集聚于靠近 i轴的区域，表明不论初始条件如何，大部分人都会感染疫情并康复。

对于小白来说，比较容易理解 2.4 节图中变量随时间的变化曲线，而对于本节相轨迹方法的思想、方法和图形都会觉得不容易理解，甚至感到困惑。虽然相轨迹的每一条线也对应着 t 从 t0 到 tend 的过程，但为什么要这么画，为什么轨迹这么怪怪的呢？

相轨迹图看上去比较怪，也不容易理解，是因为忽略时间轴而着重关注两个变量之间的关系，这种视角与我们日常观察问题和思考问题的习惯完全不同。也正是因为这个原因，相轨迹图能反映出时间变化曲线图中难以表达的一些重要特征。

例如，患病者比例在 $s=1/σs=1/\sigma$ 时达到峰值，即使把不同 $σ\sigma$ 值下的患病者比例的时间变化曲线放在一张图中也无法观察到这一特征。进一步地，既然在 $s=1/σs=1/\sigma$ 时达到峰值，那么 $s_0$ 与 $1/σ1/\sigma$ 的关系自然就成为重要的分界线，并在图中可以观察到分界线两侧具有明显不同的特征。

有了对这些特征的认识和把握，才能选择不同的参数条件，在时间变化曲线图上进行比较系统的比较。要知道 SIR 模型中有 $λ、μ\lambda、\mu$ 两个参数，还有 $i_0、s_0$ 两个初始条件，共有 4 个可以设置的参数都会影响微分方程的解，也就是会影响患病者、易感者比例的时间变化曲线。其中的各种组合无穷无尽，如果没有恰当的研究方法、不能把握内在的规律，即使在几十、几百组参数条件下进行模拟，仍然只是盲人摸象、管中窥豹。

看到这里，小白同学可能会对相轨迹研究的意义有所认识，但还是对这种分析方法难以理解望而却步。

没关系，还记得我们在”09 微分方程模型“中说的吗：

不会从问题建立微分方程模型怎么办，不会展开参数对稳定性、灵敏度的影响进行讨论怎么办？谁让你自己做呢，当然是先去找相关专业的教材、论文，从中选择比较接近、比较简单的理论和模型，然后通过各种假设强行将题目简化为模型中的条件，这就可以照猫画虎了。

5. SIR 模型讨论

最后，我们简单总结一下 SIR 模型的特点：

SIR 模型是一个单向模型，易感者（S）不断转变为患病者（I），患病者（I）不断转变为康复者（R），因此易感者比例 s(t) 单调递减，康复者比例 r(i) 单调递增。
若 $s0>1/σs_0>1/\sigma$ ，患病者比例 i(t) 先增长，当 $s0=1/σs_0=1/\sigma$ 时达到峰值，然后下降，最终为 0；若 $s0<1/σs_0<1/\sigma$ ，患病者比例 i(t) 单调递减，最终为 0。
不论初始条件如何，患病者数量最终都会清零。
$1/σ1/\sigma$ 是传染病蔓延的阈值，满足 $s0>1/σs_0>1/\sigma$ 才会发生传染病蔓延。因此，为了控制传染病的蔓延：
- 一方面要提高阈值 $1/σ1/\sigma$ ，这可以通过提高卫生水平来降低日接触率 $λ\lambda$ 、提高医疗水平来提高日治愈率 $μ\mu$ ；
- 另一方面要降低 $s_0$ ，这可以通过预防接种达到群体免疫来实现。

【本节完】

欢迎关注『Python小白的数学建模课 @ Youcans』原创作品

原创作品，转载必须标注原文链接：https://blog.csdn.net/youcans/article/details/117843875。

Crated：2021-06-12

欢迎关注『Python小白的数学建模课 @ Youcans』系列，持续更新
Python小白的数学建模课-01.新手必读
Python小白的数学建模课-02.数据导入
Python小白的数学建模课-03.线性规划
Python小白的数学建模课-04.整数规划
Python小白的数学建模课-05.0-1规划
Python小白的数学建模课-06.固定费用问题
Python小白的数学建模课-07.选址问题
Python小白的数学建模课-09.微分方程模型
Python小白的数学建模课-10.微分方程边值问题
Python小白的数学建模课-A1.国赛赛题类型分析
Python小白的数学建模课-A2.2021年数维杯C题探讨
Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
Python小白的数学建模课-B2. 新冠疫情 SI模型
Python小白的数学建模课-B3. 新冠疫情 SIS模型
Python小白的数学建模课-B4. 新冠疫情 SIR模型
Python小白的数学建模课-B5. 新冠疫情 SEIR模型
Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型
Python数模笔记-PuLP库
Python数模笔记-StatsModels统计回归
Python数模笔记-Sklearn
Python数模笔记-NetworkX
Python数模笔记-模拟退火算法