Python小白的數學建模課-B4. 新冠疫情 SIR模型

時間  2021-07-23
標籤 html python 編程 數組 ide 函數 工具 spa code orm 欄目 Python 简体版
原文   原文鏈接

Python小白的數學建模課-B4. 新冠疫情 SIR模型


傳染病的數學模型是數學建模中的典型問題，常見的傳染病模型有 SI、SIR、SIRS、SEIR 模型。html

SIR 模型將人羣分爲易感者（S類）、患病者（I類）和康復者（R 類），考慮了患病者治癒後的免疫能力。python

本文詳細給出了 SIR 模型微分方程、相空間分析的建模、例程、結果和分析，讓小白都能懂。編程

『Python小白的數學建模課 @ Youcans』帶你從數模小白成爲國賽達人。數組


1. 疫情傳播 SIR 模型

傳染病的傳播特性不可能經過真實的試驗開展研究，所以須要針對不一樣的傳染病傳播方式和流行特色創建相應的數學模型，並對模型進行理論研究和數值模擬。經過研究發現傳染病傳播的特徵閾值，就能夠爲預防和控制傳染病提供數據支撐和防控策略。ide

1927年，W. Kermack 在論文 「Contributions to the mathcmatical theory of epidemics」 中研究了倫敦黑死病和孟買瘟疫的流行過程，創造性地提出了 SIR 模型。函數

SI 模型和 SIS 模型將人羣分爲感染者（S類）和患病者（I類）兩類人羣，但大多數傳染病的患者在治癒後就有很強的免疫力，終身或在一段時期內再也不會被感染而變成病人。這類人羣稱爲病癒免疫的康復者（R 類）。康復者已經退出傳染系統，對於致死性疾病的死亡者也能夠用該類別描述其傳播特性。工具

SIR 模型適用於具備易感者、患病者和康復者三類人羣，能夠治癒，且治癒後終身免疫再也不復發的疾病，例如天花、肝炎、麻疹等免疫力很強的傳染病。spa

SIR 模型假設：code

  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 模型的微分方程：
由orm

\[\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 小白，又涉及到看着就心煩的微分方程組，因此咱們寧願把程序寫得累贅一些，便於讀者將程序與前面的微分方程組逐項對應。

  1. 定義初值 \(y_0\) 和 \(y\) 的定義區間 \([t_0,\ t]\)，注意初值爲數組向量 \(y_0=[i_0, s_0]\)。
  2. 調用 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模型

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程