用动力学方法模拟小区之间疫情传播

用动力学方法模拟小区之间疫情传播

这一篇本来是想做成一个研究的，但是由于找不到继续深化的点，而且之后要搞其他课题，所以把已经完成的部分记录下来。
这个研究的目的就是建立一个反映城市内部各个小区之间疫情扩散的动力学仿真模型。该模型基于经典SEIR传染病模型，可以研究和预测疫情在城市内部的传播过程。

系统动力学方法

在VENSIM软件中建立两个小区之间疫情传播的系统动力学流图，小区1的新增潜伏、感染个体与小区1到小区2的流动人口有着十分密切的关系。具体的关系可以在系统动力学流图中看出来。

设置步长为1天，总共仿真的时长为100天，假设疫情从小区1开始爆发，可以模拟出100天内小区1和小区2各自的易感、感染、潜伏、移除人数。

多个小区之间

由于系统动力学流图只能构建两个小区之间的传播关系(当多于3个小区时，小区两两之间都会产生人口的流动，流图将会变得十分复杂)，因此对多个小区之间的疫情传播也进行了仿真。

import pandas as pd
import numpy as np

小区1人口 = 8000000
小区2人口 = 6000000
小区3人口 = 5000000
od = np.array([[0,200000,400000],[200000,0,200000],[400000,300000,0]])
x = locals()
for i in range(3):
    x['小区' + str(i+1)+'人口'] = 0
r i in range(3):
    x['小区' + str(i+1)] = np.zeros((100,4))
小区1[0][0] = 2400000
小区1[0][1] = 9999
小区2[0][0] = 1800000
小区3[0][0] = 1500000
for i in range(3):
    x['小区' + str(i+1)+'流量'] = np.zeros((100,6))
三种个体 = np.zeros((3,3,3))

for t in range(1,100):
    if t <=40:
        for k in range(3):
            for i in range(3):
                for j in range(3):
                    三种个体[i][j]=[od[i][j]*(x['小区' + str(i+1)][t-1][0]/sum(x['小区' + str(i+1)][t-1])),od[i][j]*(x['小区' + str(i+1)][t-1][1]/sum(x['小区' + str(i+1)][t-1])),od[i][j]*(x['小区' + str(i+1)][t-1][2]/sum(x['小区' + str(i+1)][t-1]))]
            a1 = 0
            a2 = 0
            for m in range(3):
                a1 += 0.10*5*x['小区' + str(m+1)][t-1][1]*三种个体[k][m][0]/x['小区' + str(m+1)+'人口']
                a2 += 0.10*1*x['小区' + str(m+1)][t-1][2]*三种个体[k][m][0]/x['小区' + str(m+1)+'人口']
            x['小区' + str(k+1)+'流量'][t-1][3] = a1 + a2
            x['小区' + str(k+1)+'流量'][t-1][0] = x['小区' + str(k+1)+'流量'][t-1][3]+0.10*5*(x['小区' + str(k+1)][t-1][1]+sum(三种个体[:,k,1]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k,:,0]))/x['小区' + str(k+1)+'人口']+0.10*1*(x['小区' + str(k+1)][t-1][2]+sum(三种个体[:,k,2]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k][:][0]))/x['小区' + str(k+1)+'人口']
            x['小区' + str(k+1)+'流量'][t-1][4] = 1/7*sum(三种个体[k,:,1])
            x['小区' + str(k+1)+'流量'][t-1][1] = x['小区' + str(k+1)+'流量'][t-1][4]+1/7*(x['小区' + str(k+1)][t-1][1]-sum(三种个体[k,:,1]))
            x['小区' + str(k+1)+'流量'][t-1][5] = 0.2 *  sum(三种个体[k,:,2]) 
            x['小区' + str(k+1)+'流量'][t-1][2] = x['小区' + str(k+1)+'流量'][t-1][5]+0.2*(x['小区' + str(k+1)][t-1][2]-sum(三种个体[k,:,2]))
            x['小区' + str(k+1)][t][0] = x['小区' + str(k+1)][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][0]
            x['小区' + str(k+1)][t][1] = x['小区' + str(k+1)][t-1][1]+x['小区' + str(k+1)+'流量'][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][1]
            x['小区' + str(k+1)][t][2] = x['小区' + str(k+1)][t-1][2]+x['小区' + str(k+1)+'流量'][t-1][1]-x['小区' + str(k+1)+'流量'][t-1][2]    
            x['小区' + str(k+1)][t][3] = x['小区' + str(k+1)][t-1][3]+x['小区' + str(k+1)+'流量'][t-1][2]                                                               
    if t >40 and t<=60:
         for k in range(3):
            for i in range(3):
                for j in range(3):
                    三种个体[i][j]=0
            a1 = 0
            a2 = 0
            for m in range(3):
                a1 += 0
                a2 += 0
            x['小区' + str(k+1)+'流量'][t-1][3] = 0
            x['小区' + str(k+1)+'流量'][t-1][0] = 0+0.10*0*(x['小区' + str(k+1)][t-1][1]+sum(三种个体[:,k,1]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k,:,0]))/x['小区' + str(k+1)+'人口']+0.10*0*(x['小区' + str(k+1)][t-1][2]+sum(三种个体[:,k,2]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k][:][0]))/x['小区' + str(k+1)+'人口']
            x['小区' + str(k+1)+'流量'][t-1][4] = 1/7*sum(三种个体[k,:,1])
            x['小区' + str(k+1)+'流量'][t-1][1] = x['小区' + str(k+1)+'流量'][t-1][4]+1/7*(x['小区' + str(k+1)][t-1][1]-sum(三种个体[k,:,1]))
            x['小区' + str(k+1)+'流量'][t-1][5] = 0.2 *  sum(三种个体[k,:,2]) 
            x['小区' + str(k+1)+'流量'][t-1][2] = x['小区' + str(k+1)+'流量'][t-1][5]+0.2*(x['小区' + str(k+1)][t-1][2]-sum(三种个体[k,:,2]))
            x['小区' + str(k+1)][t][0] = x['小区' + str(k+1)][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][0]
            x['小区' + str(k+1)][t][1] = x['小区' + str(k+1)][t-1][1]+x['小区' + str(k+1)+'流量'][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][1]
            x['小区' + str(k+1)][t][2] = x['小区' + str(k+1)][t-1][2]+x['小区' + str(k+1)+'流量'][t-1][1]-x['小区' + str(k+1)+'流量'][t-1][2]    
            x['小区' + str(k+1)][t][3] = x['小区' + str(k+1)][t-1][3]+x['小区' + str(k+1)+'流量'][t-1][2]
    if t>60:
        for k in range(3):
            for i in range(3):
                for j in range(3):
                    三种个体[i][j]=[od[i][j]*(x['小区' + str(i+1)][t-1][0]/sum(x['小区' + str(i+1)][t-1])),od[i][j]*(x['小区' + str(i+1)][t-1][1]/sum(x['小区' + str(i+1)][t-1])),od[i][j]*(x['小区' + str(i+1)][t-1][2]/sum(x['小区' + str(i+1)][t-1]))]
            a1 = 0
            a2 = 0
            for m in range(3):
                a1 += 0.10*5*x['小区' + str(m+1)][t-1][1]*三种个体[k][m][0]/x['小区' + str(m+1)+'人口']
                a2 += 0.10*1*x['小区' + str(m+1)][t-1][2]*三种个体[k][m][0]/x['小区' + str(m+1)+'人口']
            x['小区' + str(k+1)+'流量'][t-1][3] = a1 + a2
            x['小区' + str(k+1)+'流量'][t-1][0] = x['小区' + str(k+1)+'流量'][t-1][3]+0.10*5*(x['小区' + str(k+1)][t-1][1]+sum(三种个体[:,k,1]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k,:,0]))/x['小区' + str(k+1)+'人口']+0.10*1*(x['小区' + str(k+1)][t-1][2]+sum(三种个体[:,k,2]))*(x['小区' + str(k+1)][t-1][0]-sum(三种个体[k][:][0]))/x['小区' + str(k+1)+'人口']
            x['小区' + str(k+1)+'流量'][t-1][4] = 1/7*sum(三种个体[k,:,1])
            x['小区' + str(k+1)+'流量'][t-1][1] = x['小区' + str(k+1)+'流量'][t-1][4]+1/7*(x['小区' + str(k+1)][t-1][1]-sum(三种个体[k,:,1]))
            x['小区' + str(k+1)+'流量'][t-1][5] = 0.2 *  sum(三种个体[k,:,2]) 
            x['小区' + str(k+1)+'流量'][t-1][2] = x['小区' + str(k+1)+'流量'][t-1][5]+0.2*(x['小区' + str(k+1)][t-1][2]-sum(三种个体[k,:,2]))
            x['小区' + str(k+1)][t][0] = x['小区' + str(k+1)][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][0]
            x['小区' + str(k+1)][t][1] = x['小区' + str(k+1)][t-1][1]+x['小区' + str(k+1)+'流量'][t-1][0]-x['小区' + str(k+1)+'流量'][t-1][1]
            x['小区' + str(k+1)][t][2] = x['小区' + str(k+1)][t-1][2]+x['小区' + str(k+1)+'流量'][t-1][1]-x['小区' + str(k+1)+'流量'][t-1][2]    
            x['小区' + str(k+1)][t][3] = x['小区' + str(k+1)][t-1][3]+x['小区' + str(k+1)+'流量'][t-1][2]

from matplotlib import pyplot as plt 
plt.plot(小区1[:,2])
plt.plot(小区2[:,2])
plt.plot(小区3[:,2])

得到的结果如上图所示。

结果分析

对一些措施及其力度进行了调整，得到以下的疫情传播趋势。