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

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

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

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

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

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

---

1. 疫情传播 SIR 模型

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

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

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

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

SIR 模型假设：

1. 考察地区的总人数 N 不变，即不考虑生死或迁移；
2. 人群分为易感者（S类）、患病者（I类）和康复者（R 类）三类；
3. 易感者（S类）与患病者（I类）有效接触即被感染，变为患病者（I类）；患病者（I类）可被治愈，治愈后变为康复者；康复者（R类）得到终身免疫再也不易感；无潜伏期；
4. 将第 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\)。
5. 日接触数 \(\lambda\)，每一个患病者天天有效接触的易感者的平均人数；
6. 日治愈率 \(\mu\)，天天被治愈的患病者人数占患病者总数的比例，即平均治愈天数为 \(1/\mu\)；
7. 传染期接触数 \(\sigma = \lambda / \mu\)，即每一个患病者在整个传染期内有效接触的易感者人数。

SIR 模型的微分方程：
由

\[\begin{cases} \begin{align} & 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{align} \end{cases} \]

得：

\[\begin{cases} \begin{align*} & \frac{ds}{dt} = -\lambda s i , &s(0)=s_0\\ & \frac{di}{dt} = \lambda s i - \mu i , &i(0)=i_0\\ \end{align*} \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,p1,p2,..)\) 包括可变参数 p1,p2.. 时，经过 args =(p1,p2,..) 能够将参数p1,p2.. 传递给导数函数 func。

odeint() 的返回值：

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

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

1. 导入 scipy、numpy、matplotlib 包。
2. 定义导数函数 \(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*i
    return dy_dt

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

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

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

3. 定义初值 \(y_0\) 和 \(y\) 的定义区间 \([t_0,\ t]\)，注意初值为数组向量 \(y_0=[i_0, s_0]\)。
4. 调用 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*i
    return dy_dt

def dySIR(y, t, lamda, mu):  # SIR 模型，导数函数
    i, s = y
    di_dt = lamda*s*i - mu*i  # di/dt = lamda*s*i-mu*i
    ds_dt = -lamda*s*i  # ds/dt = -lamda*s*i
    return 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 的运行结果，参数和初值为：\(\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_\infty=1-\mu/\lambda = 0.6\)，代表疫情稳定，并将长期保持必定的患病率，称为地方病平衡点。

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

---

3. SIR 模型参数的影响

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

3.1 初值条件 \(i_0, s_0\) 的影响

下面考察初值条件 \(i_0, s_0\) 的影响。固定参数 \(\lambda=0.2, \mu=0.02\)不变，不一样初值条件下 \(i(t), s(t)\) 的变化曲线以下图所示。

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

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

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

3.2 日接触率 \(\lambda\) 的影响

首先考察日接触率 \(\lambda\) 的影响。固定参数 \(\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\) 的影响。固定参数 \(\lambda=0.2, i_0=0.002,s_0=1-i_0\)不变， \(\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 能够获得：

\[\frac{di}{ds} = \frac{1}{\sigma s} - 1 ,\; i(s=s_0)=i_0 \]

该微分方程的解为：

\[i = (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 = y
    di_dt = lamda*s*i - mu*i  # di/dt = lamda*s*i-mu*i
    ds_dt = -lamda*s*i  # ds/dt = -lamda*s*i
    return 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 的运行结果（\(\lambda =0.2，\mu=0.08，1/\sigma=0.4\)），是 SIR 模型的相轨迹图。

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

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

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

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

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

例如，患病者比例在 \(s=1/\sigma\) 时达到峰值，即便把不一样 \(\sigma\) 值下的患病者比例的时间变化曲线放在一张图中也没法观察到这一特征。进一步地，既然在 \(s=1/\sigma\) 时达到峰值，那么 \(s_0\) 与 \(1/\sigma\) 的关系天然就成为重要的分界线，并在图中能够观察到分界线两侧具备明显不一样的特征。

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

看到这里，小白同窗可能会对相轨迹研究的意义有所认识，但仍是对这种分析方法难以理解望而却步。不要紧，还记得咱们在”09 微分方程模型“中说的吗：

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

---

5. SIR 模型结果讨论

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

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

【本节完】

版权声明：

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

原创做品，转载必须标注原文连接：(https://www.cnblogs.com/youcans/p/14978514.html)。

Copyright 2021 Youcans, XUPT

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小白的数学建模课-B2.新冠疫情 SI模型
Python小白的数学建模课-B3.新冠疫情 SIS模型