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 模型假設：
- 考察地區的總人數 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/\mu\)；
- 傳染期接觸數 \(\sigma = \lambda / \mu\)，即每個患病者在整個傳染期內有效接觸的易感者人數。
SIR 模型的微分方程：
由
得：
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 模型的編程步驟
- 導入 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*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 小白，又涉及到看着就心煩的微分方程組，所以我們寧願把程序寫得累贅一些，便於讀者將程序與前面的微分方程組逐項對應。
- 定義初值 \(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*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 個可以調整的參數條件都會影響微分方程的解，也就是會影響患病者、易感者比例的時間變化曲線。其中的各種組合無窮無盡，如果沒有恰當的研究方法、不能把握內在的規律，即使在幾十、幾百組參數條件下進行模擬，仍然只是盲人摸象、管中窺豹。
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 可以得到：
該微分方程的解為：
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 模型的性質：
- 任一條 i-s 曲線都收斂於s軸上的一點，即 \(i_\infty=0\)，表明不論初始條件如何，患病者終將清零。
- 患病者比例在 \(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\)，患病者比例單調遞減，終將清零。
- 易感者比例單調遞減，易感者的最終比例是相軌跡與 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 個可以設置的參數都會影響微分方程的解，也就是會影響患病者、易感者比例的時間變化曲線。其中的各種組合無窮無盡，如果沒有恰當的研究方法、不能把握內在的規律，即使在幾十、幾百組參數條件下進行模擬，仍然只是盲人摸象、管中窺豹。
看到這里，小白同學可能會對相軌跡研究的意義有所認識，但還是對這種分析方法難以理解望而卻步。沒關系，還記得我們在”09 微分方程模型“中說的嗎：
不會從問題建立微分方程模型怎么辦，不會展開參數對穩定性、靈敏度的影響進行討論怎么辦？誰讓你自己做呢，當然是先去找相關專業的教材、論文，從中選擇比較接近、比較簡單的理論和模型，然后通過各種假設強行將題目簡化為模型中的條件，這就可以照貓畫虎了。
5. SIR 模型結果討論
最后，我們簡單總結一下 SIR 模型的特點：
- SIR 模型是一個單向模型，易感者（S）不斷轉變為患病者（I），患病者（I）不斷轉變為康復者（R），因此易感者比例 s(t) 單調遞減，康復者比例 r(i) 單調遞增。
- 若 \(s_0>1/\sigma\)，患病者比例 i(t) 先增長，當 \(s_0=1/\sigma\) 時達到峰值，然后下降，最終為 0；若 \(s_0<1/\sigma\)，患病者比例 i(t) 單調遞減，最終為 0。
- 不論初始條件如何，患病者數量最終都會清零。
- \(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模型