开放式IT项目实战训练
题 目 学 院 姓 名 学 号 指导教师 完成时间
说明文档
基于AnyLogic的传染病传播仿真建模 计算机与通信学院 陶婷婷 10240541 赵宏 2013年5月11日
： ： ： ： ：
：
目 录
第一章 需求分析 .................................................................................................. 1 第二章 总体设计 .................................................................................................. 2
2.1系统框架图 .................................................................................................. 2 2.2 系统程序功能实现步骤的详细说明 ......................................................... 2 第三章 详细设计 .................................................................................................. 6 第四章 功能测试 .................................................................................................. 9 第五章 源码 ........................................................................................................ 12 总 结 .................................................................................................................... 15 致 谢 .................................................................................................................... 16 参 考 文 献 .......................................................................................................... 17
第一章 需求分析
传染病模型有着悠久的历史，一般认为始于1760年Daniel Bernoulli在他的一篇论文中对接种预防天花的研究。真正的确定性传染病数学模型研究的前进步伐早在20世纪初就开始了，Hamer, Ross等人在建立传染病数学模型的研究中做出了大量的工作.直到1927年Kermack与McKendrick在研究流行于伦敦的黑死病时提出了的SIR仓室模型，并于1932年继而建立了SIS模型，在对这些模型的研究基础上提出了传染病动力学中的阐值理论.Kermack与McKendrick的SIR模型是传染病模型中最经典、最基本的模型，为传染病动力学的研究做出了奠基性的贡献.该模型模拟了易感人群、感染人数和恢复人数三者之间数目的变化，并通过一定的图形界面展示出来。
SIR模型中S表示易感者，I表示感染者，R表示移出者。模型中把传染病流行范围内的人群分成三类：S类，易感者（Susceptible），指未得病者，但缺乏免疫能力，与感病者接触后容易受到感染；I类，感病者（Infective），指染上传染病的人，它可以传播给S类成员；R类，移出者（Removal），指被隔离，或因病愈而具有免疫力的人。
随着社会的发展和在民众物质，经济的提高，传染病的种类与传播途径也越来越多。由于如今各地方人流量越来越大、人群来源更复杂等的特性，使得传染病通过各种方式快速地波及全球的机会越来越大，传染病流行的全球化越发明显。单是2003年的SARS，疫情单通过航空网络传播到全球25个国家，共有8000人被感染，775人死亡，经济损失400亿美元。另外，当国际性的传染病爆发时，由于人口平均所占空间较少，极大的增加了潜在传播传染病的机率从而使得干预准备的时间减少，因此分析传染病模型，成为解决疾病的一种重要手段，对防止传染病的传入传出意义重大，同时人们也对传染病监控的模式、重点、反应时限等均提出了特殊的要求。 针对流感大流行，许多干预措施，如关闭学校、工厂、禁止大型集会、旅行限制、边界检疫、病例诊断和治疗、隔离患者、密接者追踪与管理、接种疫苗、风险沟通、加强公共卫生和个人卫生等，都被广泛采用。通过传染病传播动力学模型，可以研究这些干预措施对传染病传播的影响，进而评估各种干预措施，及围堵策略的效率。
本文就是基于这样的背景，旨在通过分析传染病传播扩散的动力学模型，探索如何在生活中科学有效地防控传染病的扩散。
1
第二章 总体设计
2.1系统框架图
仓室模型的主要思想是将人群划分为若干个类（仓室），分别代表处于不同疾病状态的人群，如易感者、感染者和移出者，然后采用数学手段建立这些变量的动力学方程，进而研究疾病的传播动力学过程。SIR 仓室模型这一简洁的建模思想获得了研究者的青睐，很快就被发展成各种具有更多复杂仓室结构的模型，在20 世纪此后的几十年中，一直占据着传染病动力学建模的主流地位,如下表2.1.1： 表2.1.1 仓室模型表 模型 条件 SI模型 患病后难以治愈 SIS模型 1. 建模简单，容易实现模SIR模患病者治愈后获利终身免疫力 拟 型 病人康复后只有暂时免疫力，2. 假设人群SIRS单位时间内将有部分康复者丧接触概率均模型 等，只能适用失免疫力而可能再次被感染 在三类人群中增加一类，感染而未发病于小规模人者(Exposed)，可在SIR或SIRS模型的群、 基础上得到更复杂的SEIR或SEIRS模型 患病后可以治愈，恢复者不具有免疫力 特点 SIR模型 在不考虑出生与死亡等种群动力学因素的情况下，传染病若无潜伏期 考虑传染病的潜伏期 模型主要描述的是以下易感人群，感染人数，恢复人数三者之间的关系，如下图2.1.2：
感染率 恢复率
图2.1.2 联系图
还可以设计简单的时间堆叠图和饼图来表示它们三者之间的变化关系。
2.2 系统程序功能实现步骤的详细说明
Anylogic可以实现多种的建模仿法，并且可以将多种方法融合在一起使用，
2
最近本人了解了一下Anylogic的主体建模方法，总体感觉简单易用，功能强大，与老牌的Swarm相比，具备很多非常显著的特点。% M\" N9 U* h9 h& w7 N% AAnylogic的主体仿真总体感觉，学习难度不大，功能却很强大，如果能够将系统动力学、连续系统仿真和离散事件仿真结合起来将威力无穷。目前，Anylogic引入国内的时间比较短，可供参考的资料比较少，学习起来有一定的难度，光参考案例和手册，很难掌握其精华，我学习了几天认为，sir模型实现功能有： 一﹑模型假设
在疾病传播期内所考察的地区范围不考虑人口的出生、死亡、流动等种群动力因素。总人口数N(t)不变，人口始终保持一个常数N。人群分为以下三类：易感染者(Susceptibles)，其数量比例记为s(t)，表示t时刻未染病但有可能被该类疾病传染的人数占总人数的比例；感染病者(Infectives)，其数量比例记为i(t)，表示t时刻已被感染成为病人而且具有传染力的人数占总人数的比例；恢复者(Recovered)，其数量比例记为r(t)，表示t时刻已从染病者中移出的人数（这部分人既非已感染者，也非感染病者，不具有传染性，也不会再次被感染，他们已退出该传染系统。）占总人数的比例。
病人的日接触率（每个病人每天有效接触的平均人数）为常数λ，日治愈率（每天被治愈的病人占总病人数的比例）为常数μ，显然平均传染期为1／μ，传染期接触数为σ=λ／μ。该模型的缺陷是结果常与实际有一定程度差距，这是因为模型中假设有效接触率传染力是不变的。 二﹑模型构成
在以上三个基本假设条件下，易感染者从患病到移出的过程框图2.2.1表示如下：
s i r λsi μi 图2.2.1 过程框图
在假设1中显然有：
s(t) ? i(t) ? r(t) ? 1 （1） 对于病愈免疫的移出者的数量应为：
d? N??Ni （2）
dts0（s0＞0），i0（i0＞不妨设初始时刻的易感染者，染病者，恢复者的比例分别为0），r0?0.
SIR基础模型用微分方程组表示如下：
3
?di???si??i?dt?ds????si （3） ?dt?d????i?dts(t) ，i(t)的求解极度困难，在此我们先做数值计算来预估计s(t) ， i(t)的一般变化规律。
二﹑数值计算
在方程（3）中设λ=1，μ=0.3，i(0)= 0.02，s(0)=0.98，用MATLAB软件编程： function y=ill(t,x)
a=1;b=0.3;
y=[a*x(1)*x(2)-b*x(1);-a*x(1)*x(2)]; ts=0:50;
x0=[0.20,0.98];
[t,x]=ode45('ill',ts,x0); plot(t,x(:,1),t,x(:,2)) pause
plot(x(:,2),x(:,1))
输出的简明计算结果列入表（2.2）。i(t) , s(t)的图形以下两个图形，i~s图形称为相轨线,初值i(0)=0.02,s(0)=0.98相当于图2中的P0点,随着t的增,(s,i)沿轨线自右向左运动.由表2.2.1、图2.2.1、图2.2.2可以看出,i(t)由初值增长至约t=7时达到最大值,然后减少,t→∞,i→0,s(t)则单调减少,t→∞,s→0.0398. 并分析i(t),s(t)的一般变化规律。
表2.2.1 i(t),s(t)的数值计算结果
t i(t) s(t) t i(t) s(t) 0 0.0200 0.9800 9 0.2863 0.1493 1 0.0390 0.9525 10 0.2418 0.1145 2 0.0732 0.9019 15 0.0787 0.0543 3 0.1285 0.8169 20 0.0223 0.0434 4 0.2033 0.6927 25 0.0061 0.0408 5 0.2795 0.5438 30 0.0017 0.0401 6 0.3312 0.3995 35 0.0005 0.0399 7 0.3444 0.2839 40 0.0001 0.0399 8 0.3247 0.2027 45 0 0.0398
4
图2.2.2 病人比例图
图2.2.3 i-s图形
5
第三章 详细设计
第一步：利用anylogic仿真软件创建模型 第二步：设置Main参数
1. 拖入四个Parameter，分别命名为ContactRate、Infectivity、TotalPopulation、AverageIllnessDuration，Type都为double，Default Value分别为5、0.05、1000和15，分别表示接触率、感染性、总人数和平均患病时间。
2. 拖入三个Stock Variable，分别命名为Susceptible、Infectious、Recovered，Initial Value分别为TotalPopulation-1、1和null。拖入两个Flow Aux Variable，分别命名为InfectionRate和RecoveryRate，值分别为Infectious * ContactRate * Susceptible / TotalPopulation * Infectivity和Infectious / AverageIllnessDuration。
3. 通过双击建立彼此之间的联系，建好后如下图3.1所示：
图 3.1 Main的参数
上面五个连接的变量分别表示的意义为易感人群、感染率、感染人数、恢复率和恢复人数。
1. 拖入三个Data Set，分别命名为dsSusceptible、dsInfectious、dsRecovered，勾中Use time as horizontal axis value，Vertical axis value填入Susceptible、Infectious和Recovered，都填入keep up to 101 latest samples。三个Data Set分别统计易感人数、感染人数和恢复人数。
第三步：设置Main界面
1. 拖入一个Pie chart，添加三个Data Set，Title分别为Susceptible、Infectious和Recovered，Value也分别填入Susceptible、Infectious和Recovered。选择Update automatically，Begin at time 0，Recurrence time：1。
2. 拖入一个Time Stack Chart，添加三个Data Set，Title分别为Susceptible、Infectious和Recovered，Data set分别填入dsSusceptible、dsInfectious和dsRecovered。Time Window为100，Vertical Scale为100%，选择Update automatically，Begin at time 0，Recurrence time：1。
3. 设置好的Main界面如下图3.2所示：
6
图3.2 Main类界面
第四步：设置Simulation
1.拖入三个Plain Variable，分别命名为ContactRate、Infectivity、AverageIllnessDuration，Type都为double，Initial value分别为5、0.05、15。
2.Simulation的设置界面下的General中的ContactRate、Infectivity和AverageIllnessDuration分别填入ContactRate、Infectivity和AverageIllnessDuration。
3.拖入三个Slider，Minimum Value都为0，Maximum Value分别为20、1、100，Default Value分别为ContactRate、Infectivity和AverageIllnessDuration，Action处分别填入ContactRate = value;、Infectivity = value;和AverageIllnessDuration = value;。旁边的Text的Dynamic中的Text处分别填入format( ContactRate )、format( Infectivity )和format( AverageIllnessDuration )。
4.在默认的按钮中的标签写入 Run the model and switch to Simulation view， 标签处写入run();
getEngine().getPresentation().setPresentable( getEngine().getRoot() ); 表示点击后运行仿真并跳转到 Main 界面。 5.拖入一个Time plot，Name为plot，Time Window为100，Vertical Scale选择Auto，选中Do not update automatically。
6.回到Main界面，拖入一个Button，Label为Save results and Finish，Action处填入如下代码。
//初始化Simulation中的TimePlot
TimePlot plot = ((Simulation)getEngine().getExperiment()).plot;
//点击该按钮后在Simulation中的Timeplot中添加一个表示感染人数的Data Set，以下为设置该Data Set的属性（具体参见addDataSet函数API）。
7
plot.addDataSet( dsInfectious,
\"Infectious when CR = \" + ContactRate + \" Inf = \" + Infectivity + \" AID = \" + AverageIllnessDuration,
spectrumColor( plot.getCount(), 5 ), true, Chart.INTERPOLATION_LINEAR, 1, Chart.POINT_NONE );
getEngine().getExperiment().stop();
在Model Time选项栏中的Time units选择days，Stop选择Stop at specified time，Start time为0，Stop time为100。该模型通过按钮中的函数控制，使得在一张图中显示出不同参数设置下的感染人数变化情况，形式比较新颖，对后续的参数分析也有很重要的意义。具体如下图3.3所示：
图3.3 Simulation的设置
8
第四章 功能测试
模型验证：
上世纪初在印度孟买发生的一次瘟疫中几乎所有病人都死亡了。死亡相当于移出传染系统，有关部门记录了每天移出者的人数，即有了Kermack等人用这组数据对SIR模型作了验证。
首先，由方程（2），（3）可以得到
dsd???si????si???? ，两边积分得dtdtd?的实际数据，dt?ss01ss??rds????rd ?lns???r??er0?0rs?0ss0所以：s?t??s0e??r?t? (12) ?dr??i???1?r?s???1?r?s0e??r (13) dt??当 r?
1
?时，取（13）式右端e??rTaylor展开式的前3项得：
?s0?2r2?dr ????1?r?s0??s0r?2?? dt??在初始值?0=0 下解高阶常微分方程得： r?t??1????t????s??1??th????? (14) 2?0s0???2??s0??12其中?2??s0??1??2s0i0?2，th??? 从而容易由（14）式得出：
d??dt?2s0?2ch?????2??2?2??? (15)
然后取定参数 s0, σ等，画出（15）式的图形，如图2.2.3中的曲线，实际数据在图中用圆点表示，可以看出，理论曲线与实际数据吻合得相当不错。 被传染比例的估计
在一次传染病的传播过程中，被传染人数的比例是健康者人数比例的初始值s0与sm之差，记作x,即
x?s0?sm (16) 当i0很小，s0接近于1时，由（9）式可得
9
x??x??ln?1??0 (17) ????s0?1取对数函数Taylor展开的前两项有
?1xx?1???s?2s2?00????0 (18) ??1?时（18）式给出
记 s0?1?? , ? 可视为该地区人口比例超过阈值
?的部分。当 ??1?
1??x?2s0??s0???2? (19)
???这个结果表明，被传染人数比例约为?的2倍。对一种传染病，当该地区的
1卫生和医疗水平不变，即?不变时，这个比例就不会改变。而当阈值提高时，??减小，于是这个比例就会降低。
用anylogic运行模型结果会出现如下界面返回最初界面时三者之间会有曲线的变化关系，则表明构造的模型合理，最后运行界面如下图4.1：
10
图4.1 运行界面
11
第五章 源码
先是定义了4个参数ContactRate、Infectivity、TotalPopulation、AverageIllnessDuration以及对它们设定了初始值部分代码如下： public double _TotalPopulation_DefaultValue_xjal() { final Main self = this; return 1000 ; }
public double _Infectivity_DefaultValue_xjal() { final Main self = this; return 0.05 ;
}
public double _ContactRate_DefaultValue_xjal() { final Main self = this; return 5 ; }
pblic double _AverageIllnessDuration_DefaultValue_xjal() { final Main self = this; return 15 ; }
然后定义了三个存量以及两个流量来控制它们三者之间的变化，而且对它们分别设置了初始值：
public void _assign_Susceptible_Formula_xjal() { Susceptible = TotalPopulation - 1 ; }
public void _assign_Infectious_Formula_xjal() { Infectious = 1 ; }
public void _assign_InfectionRate_Formula_xjal() { InfectionRate =
Infectious * ContactRate * Susceptible / TotalPopulation * Infectivity ; }
public void _assign_RecoveryRate_Formula_xjal() { RecoveryRate =
Infectious / AverageIllnessDuration ; }
画出时间堆叠图里面加入三个数据项其中Susceptible的代码如下： public void update( DataSet _d ) {
if ( time() == _lastUpdateX ) { return; }
12
double _a = Susceptible ;
_d.add( time(), _a ); _lastUpdateX = time(); }
接下来再画出饼图在界面显示出它们三者所占的比例，设置为自动更新且更新时间为1，代码如下：
private double _chart_DataItem0Value() { return Susceptible ; }
if ( _e == _chart_autoUpdateEvent_xjal ) return 1；
然后再添加一个按钮显示结果，如下：
TimePlot plot = ((Simulation)getEngine().getExperiment()).plot; plot.addDataSet( InfectiousDS,
\"Infectious when CR = \" + ContactRate + \" Inf = \" + Infectivity + \" AID = \" + AverageIllnessDuration,
spectrumColor( plot.getCount(), 5 ), true, Chart.INTERPOLATION_LINEAR, 1, Chart.POINT_NONE );
getEngine().getExperiment().stop();
最后设计Simulation界面先拖入三个滑块通过链接以控制三个参数值，如下： public double getShapeControlDefaultValueDouble( int _shape, int index ) { switch(_shape) {
case slider: return ContactRate ; case slider1: return Infectivity ;
case slider2: return AverageIllnessDuration ; }
添加一个时间窗绘图，代码如下： plot = new TimePlot(
Simulation.this, true, 380, 110, 400, 440,
null, null, 40, 20,
340, 190, white, black, black, 210, Chart.SOUTH,100 , null, Chart.SCALE_AUTO,
0, 0, Chart.GRID_DEFAULT, Chart.GRID_DEFAULT, darkGray, darkGray, _items, _titles, _appearances ); }
13
添加一个控制按钮回到最初的界面，并且改变最初的界面，如下： public boolean isShapeControlEnabled( int _shape, int index ) { switch( _shape ) { case button: return
getEngine().getState() == Engine.IDLE ;
default: return super.isShapeControlEnabled( _shape, index ); } } run();
getEngine().getPresentation().setPresentable( getEngine().getRoot() );
最后插入一张显示anylogic标志的图片。
14
总 结
该模型采用了数值计算,图形观察与理论分析相结合的方法,先有感性认识 (表2.1.1,图2.2.1),再进行数值验证和估算,可以看作计算机技术与建模方法的巧妙配合。可取之处在于它们比较全面地达到了建模的目的,即描述传播过程、分析感染人数的变化规律,预测传染病高潮到来时刻,度量传染病蔓延的程度并探索制止蔓延的手段和措施。主要成果有：
(1) 综述了国内外传染病研究情况、仿真技术研究及其应用情况，并介绍了仿真软件AnyLogic中SIR模型。
(2) 研究并改进了SIR模型，使之更加适应于各种传染病传播场所和途径。 (3) 提出相关传染病防控措施，同时定量分析了各个防控措施的效果。
15
致 谢
一个月的模型设计,有过山穷水尽的困惑；有过柳暗花明的惊喜，在不知道从哪下手的时候我战胜了自己，勇敢面对人生的信念是我战胜自己的根本。回忆当初刚开始学习anylogic时候我对它怀着满腔的热情，同学帮忙安装。现在我已经基本能够自己做成模型，这里我要感谢我的老师，是赵老师的悉心教导才使我打好了基础，在我最最困惑的时候又是学长给与我很大的帮助，使我走出困惑。
感谢赵老师的辛勤指导，为我解决设计过程中遇到的困难，不仅使我掌握到了一些关于仿真模型设计的基本知识，同时让我学习到了更多课外实用性的知识，大大加强了我对学习的兴趣，加深了对仿真的了解，并且顺利完成了本次的创新课程设计，真诚的感谢你。
16
参 考 文 献
[1] 高景德,王祥珩,李发海.交流电机及其系统的分析[M]. 北京:清华大学出版社,1993,8.120-125.
[2] 李永东.PWM供电的异步电机电压定向矢量控制[J]. 电气传动,1991,4(1):52-56. [3] 张强.浅析电力行业企业的整改措施[N]. 中国电力报,2004,10-11. [4] 马知恩等.传染病动力学的数学建模与研究. 北京:科学出版社，2004. [5] 于飞.机场生产流程仿真研究[D]. 南京:南京航空航天大学，2005.
[6] 倪桂明,杨东援.机场系统计算机仿真研究的应用与发展[J]. 系统工程学报，2002. [7] 曹燕. 机场值机柜台的优化管理市场周刊[J]. 2004.
[8] 李洪涛等.国内民航市场构成及消费行为分析. 中国民航航空,2010.
[9] Kermack, W.O. and McKendrick, A.G. Contributions to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London, Series A, 1927, (115): 700-721.
[10] Hethcote, H. W. (2000). \"The mathematics of infectious diseases.\" Siam Review 42(4):599-653.
17
因篇幅问题不能全部显示，请点此查看更多更全内容