1500字范文 > 新型冠状病毒SIR预测模型 MATLAB代码

新型冠状病毒SIR预测模型 MATLAB代码

时间：2019-06-24 03:02:36

What's the point of TOC

SIR模型简介数据获得和预处理求解思路以及代码实现结果展示回复评论区

SIR模型简介

1、关于传染病的数学模型有哪些？ - 知乎用户的回答 - 知乎
/question/367466399/answer/982558966
2、关于传染病的数学模型有哪些？ - 酱紫君的回答 - 知乎
/question/367466399/answer/982597090
3、徐宝春. 基于SIR模型的SARS传染病研究[D]. .

模型及其动力学方程如下：

数据获得和预处理

数据获取方式为官方渠道。模型中的ItI_tIt和RtR_tRt的实际值由下列公式得到。

求解思路以及代码实现

易感人群初值S0S_{0}S0难以估计，为待定量。动力学方程组中的β\betaβ和γ\gammaγ同样为待定量。解决思路如下：

1、待定S0S_{0}S0，β\betaβ和γ\gammaγ。目标函数为误差函数（见4），使用遗传算法寻找其最小值。

2、输入初值条件：第0天时的易感人数S0S_0S0，感染人数I0I_0I0和恢复人数R0R_0R0。

3、求解SIR模型的控制方程，获得第1~n天的SIR预测值：其为常微分方程组（ODE），可使用龙格库塔4阶方法。

4、求出实际值与预测值的误差。这部分误差由两部分构成，感染人数ItI_tIt和恢复人数RtR_tRt。而且由于数据源的滞后性和消息来源的不同，需要自行调整权重值。

遗传算法的目标函数代码：

function ra=optSIR(X)% beta=0.00002;% gamma=1/20;beta=X(1)/1e6;gamma=X(2)/1e2;S0=X(3)*1e4;x0=[S0,25,2];ts=0:1:300;[t,x] = ode45(@(t,x) SIRModel(t,x,beta,gamma), ts, x0);I =[7621 5768 4696 3747 2857 2364 2042 1721 1458578526502440227169412827 25];R =[730 616 446 362 358 275 219 184 1321209270553129211714 02];T =[35 34 33 32 31 30 29 28 27 26252423171615 0];intensifiedPoint=15;iP=intensifiedPoint;omega=[iP iP iP iP 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];eps_I=0;eps_R=0;for i=1:length(T)eps_I=eps_I+omega(i)*(x(T(i)+1,2)-I(i))^2/I(i)^2;eps_R=eps_R+omega(i)*(x(T(i)+1,3)-R(i))^2/R(i)^2;%eps_I=eps_I+omega(i)*(x(T(i)+1,2)-I(i))^2;%eps_R=eps_R+omega(i)*(x(T(i)+1,3)-R(i))^2;endomegaI=1e0;ra=(omegaI*sqrt(eps_I)+sqrt(eps_R))/length(T);IsPlot=0;if IsPlot==1plot(t,x(:,1),'k',t,x(:,2),'r',t,x(:,3),'b');hold onxlabel('时间/天');ylabel('人数');legend('S','I','R');plot(T,I,'ro');hold onplot(T,R,'bo');hold offendendfunction y=SIRModel(t,x,beta,gamma)y=[-beta*x(1)*x(2),beta*x(1)*x(2)-gamma*x(2),gamma*x(2)]';end

附遗传代码工具箱使用设置简略过程

画图代码

clearIsPlot=1;ra=[24.777628430066810.0194108949530886750.7805747296079624];beta=ra(1)*1e-6;gamma=ra(2)*1e-2;S0=ra(3)*1e4;x0=[S0,25,2];ts=0:1:300;[t,x] = ode45(@(t,x) SIRModel(t,x,beta,gamma), ts, x0);I =[7621 5768 4696 3747 2857 2364 2042 1721 1458578526502440227169412827 25];R =[730 616 446 362 358 275 219 184 1321209270553129211714 02];I=I/S0;R=R/S0;T =[35 34 33 32 31 30 29 28 27 26252423171615 0];% The data is from Wuhan city only. The data source cannot be revealed because it is related to sensitive word (Wuhan Wei Jian Wei and Hubei WJW).if IsPlot==1figurex=x/S0;plot(t,x(:,1),'k',t,x(:,2),'r',t,x(:,3),'b');hold onxlabel('时间/天');ylabel('比例');plot(T,I,'ro');hold onplot(T,R,'bo');hold offlegend('易感人群','感染人群预测值','康复人群预测值','感染人群实际值','康复人群实际值');axis([0,300,0,1])endfunction y=SIRModel(t,x,beta,gamma)y=[-beta*x(1)*x(2),beta*x(1)*x(2)-gamma*x(2),gamma*x(2)]';end