《埃博拉病毒的传播预测与控制美赛.doc》由会员分享，可在线阅读，更多相关《埃博拉病毒的传播预测与控制美赛.doc（26页珍藏版）》请在三一文库上搜索。
1、Team # 41029 page 26 of 26埃博拉病毒的传播预测与控制 摘要2014年非洲爆发了历史上最为严重的病毒疫情-埃博拉。据科学研究报道，这个病毒一旦感染人体，将有着高达90%以上的死亡率，这是一种世上最厉害的感染病毒（生物安全等级为4级），如何消灭埃博拉成为当前的首要任务。当然，疾病的传播、患病人口的预测、药物的生产和运输，都是消灭埃博拉必须考虑的因素。根据病毒传播率、感染者人数的预测、药物的合理分配和隔离人数的比重等因素，本文运用随机微分方程、产销平衡和最优控制三种算法分别建立了随机微分方程模型、线性规划模型和最优隔离控制模型。这三个模型分别解决了埃博拉病毒的传播规律、感染
2、者人数的预测问题、药物的运输问题和以隔离控制为决定性作用因素的优化问题。针对模型一：将环境因素作为随机变量，结合病毒传播率，本文建立了随机微分方程模型，对以后10个月的患病人口总数进行了预测。利用数值解方法,对埃博拉病毒感染者人数进行预测,并通过仿真过程验证了疾病传播率的一个临界值，得出能使埃博拉传播速度降低直至消亡的一个条件。针对模型二：假设几内亚、利比里亚和塞拉利昂为需求地，美国、中国、日本、俄罗斯、法国以及瑞士为药物生产地。利用产销平衡原理，建立了时间优化模型，求得产地与需求地之间的最短运输时间为15.8小时。针对模型三：本模型基于传染病模型，利用极值原理给出了最优控制的设计方案，通过仿
3、真，验证了最优控制方案的优越性。同时，由协态方程得到当系统控制变量为0.50时，隔离效果最佳，也证明了隔离是控制疾病继续传播最有效的控制措施。本文三个模型均使用的官方数据，而且内容上层层优化，互相补充，使文章所述更为具体，更为实用，为埃博拉病毒问题的解决提供了一份可靠地，可行的，可依赖的数学模型。关键词：埃博拉病毒 预测 随机微分方程 优化问题 最优隔离控制1.问题重述不用翻译1995年5月14日，扎伊尔发现罕见传染病埃博拉。2014年，埃博拉病毒首次爆发就夺走了近300人的生命,2014年再度爆发,大约4000人命丧黄泉。现在世界医学组织已经宣布：他们研究的新药物可以阻止埃博拉病毒，并非晚期
4、病人。本文疾病的传播、所需药物数量、可行的运输系统、运送的地点、生产疫苗或药物的速度和其他起决定性作用的因素考虑，建立一个符合实际的实用模型，可以达到优化消灭埃博拉或减小当前压力的目的。除了为此次比赛建立的模型解决方法外，为世界医学组织准备一份1-2页非技术信函，以用于他们的宣告。2.问题分析本文关于埃博拉病毒的传播、患病人数的预测、所需药物数量、可行的运输系统、疫苗的预防和药物的治疗等几个方面展开讨论和研究。模型一主要解决疾病的传播和患病人口预测问题。由于人口密度、周围是否有患病人群、生活环境等因素的随机性，所以将其视为随机变量。然后本文将病毒传播率作为一个高斯白噪声过程带入常微分方程，得到
5、关于埃博拉病毒传播的随机微分方程。此时不考虑人口的出生率、死亡率和人口的出入境情况，本文根据官方数据，得到2014年3月22号至今的感染者人数，从而得到一个疾病的传播率，进而预测未来10个月的感染者的总数。模型二主要解决药物的运输时间与成本的问题。由于几内亚、利比里亚和塞拉利昂这三个国家患病人数最多，所以选择这三个国家作为需求地。现在具备疫苗或药物生产能力的国家：美国、中国、日本、俄国、法国和瑞士。本文选择这六个国家作为产地。本模型只考虑在生产地和需求地之间的药物运输。首先保证各国所使用的运输机为同款运输机，在运输过程中，速度均为同等速度。接下来，本文将产销平衡模型中的成本替换成运输所用时间，
6、这样成本最低变成时间最短。然后结合模型一中的患病人口预测结果，再加上每个病人对应药量的比例系数，则计算出任意时刻所需要的药物总量。在满足各需求地需求量的前提下，本文再利用线性规划模型得到最优调运方案，即时间优化模型。模型三在模型二的基础上，分析其他可以消灭埃博拉的决定性因素。本文使用最优隔离控制法，把易感染者、染病者、治愈者、隔离者以及总人口数作为初始值代入目标函数，则会存在一个最优控制因素，再将其对应的状态解代入协态方程，得到最优控制因素隔离的确切最优解，再通过数值仿真完成对本文模型的最后优化。3随机微分传播模型根据提供的官方数据得知，目前感染者人数已达1.3万人，集中分布在几内亚、利比里亚
7、和塞拉利昂三个国家。本文针对这三个国家的患病情况，建立埃博拉病毒的随机微分模型来描述病毒的传播过程，分析并预测未来感染人数的变化规律。3.1 符号说明符号符号说明感染者人数占总人口比例埃博拉传播过程中人与人之间的接触率由于得了患埃博拉所造成的死亡率增加值埃博拉病毒的传播率平均传播率环境干扰强度布朗运动3.2模型假设l 假设埃博拉病毒在研究过程中不会发生变异。l 埃博拉病毒的研究期是2014.3.22至2015.1.22，研究对象为几内亚、利比里亚以及塞拉利昂三个国家。本次疫情是埃博拉病毒发现以来，规模最大的一次暴发流行，且感染者集中在这三个国家，基于此假设条件的模型更具有实用性、有效性、针对性
8、。l 在病毒传播期内这三个国家的总人数不变；不考虑出生和死亡因素对传播的影响；感染者病愈后不会再感染。 3.3 模型的建立与求解为了解决环境因素对病毒传播过程随机干扰的问题，本文利用随机微分方程研究该过程。并结合实际数据，再利用确定性模型估计疾病传播率。最后本文预测出这三个国家未来10个月内埃博拉病毒感染者人数，并得出其概率分布。在模型中，将传播率设定为一常数，因此得出的解是一固定曲线。埃博拉病毒感染者人数占总人口的比例满足下面的常微分方程 （1） 由于受到环境因素的随机干扰，埃博拉病毒的传播率会跟着改变。本文将传播率设为一个高斯白噪声过程1代入（1），得到埃博拉传播过程的随机微分方程，即用来
9、代替（1）中的： （2） 其中，是零均值且方差为1的高斯白噪声，和为常数，分别代表埃博拉在传播过程中的平均传播率和环境干扰的强度。于是，所得解（即埃博拉病毒感染者人数占总人口的比例）就变为一随机过程；解曲线将会随着布朗运动变化，从而显示环境因素的干扰对解的影响。随机微分方程(2)得不到它的解析解，所以必须采用随机微分方程数值解方法对其进行分析。本文采用法： （3） 其中和分别代表时间节点和轨道，是时间节点的个数，是样本轨道数。对每一条固定的轨道，根据以上公式及已知条件，可求出不同时间节点上的毒感染者人数的比例值，从而，该过程的均值为方差为。首先，在假定各参数不变的情况下，本文通过确定性模型计算
10、在几内亚、利比里亚、塞拉利昂三个国家中病毒的传播速度。的统计显示：2014年埃博拉病毒最肆虐的一年，比如感染者人数从几十例突增至13000多例。本文以2014年3月为起点，利用上述随机微分方程模型对未来10个月感染者人数的比例变化趋势给出预测。根据的统计数据，2014年3月22日至2015年1月22日埃博拉病毒未至晚期的感染者人数为13282例5；2013年末在爆发埃博拉病毒之前三个国家总人口为22131341人6。故研究期内埃博拉病毒感染者人数占总人口的比例约为。一般情况下，埃博拉病毒的平均潜伏期为发作到死亡约为18天左右7，故取8。进一步，假设这些参数10个月内不变。由确定性情形下感染者比
11、例所满足的常微分方程（1）、这三个国家的感染者人数比例，通过仿真，得出疾病传播速率。然后，在式（3）中取,，对未来10个月感染者人数比例进行预测，如图1（a）所示： （a） （b）图1确定性和随机情况下2014.32015.12感染者人数比例（a）图表示确定性情况，（b）图表示随机情形，取，共15条轨道）图1（a）给出了确定情形下三个国家自2014年3月来感染者人数比例曲线。由图1（a）可看出，模型与实际情况一致。到2015年12月，三个国家的感染者人数占总人口比例将达到。考虑病毒传播率受到环境因素的干扰，本文采用随机模型（2）进行预测。图 1（b）给出了当环境干扰强度时，三个国家在这期间感染
12、者人数比例的波动状况。由图1（b）可看出，由于环境干扰的作用，病毒感染者比例在某些情况下会超过 。然后本文利用仿真来验证模型的结果：疾病传播率的一个临界值能使埃博拉病毒传播速度下降直至埃博拉病毒灭亡。埃博拉病毒的传播率对未来感染者人数的变化有着十分重要的作用。理论分析表明，降低病毒的传播率可减少未来的感染者人数。本文通过模拟发现，传播率有所增加，这表明现实情况己经超出了我们的预料，必须高度警惕起来。图2 自2014.32015.12疾病传播率的拟合结果通过仿真比较两种稳定性，可以考虑三种情形：（a） （b）图3几乎必然指数稳定、一阶矩指数稳定及二阶矩指数稳定的区域(a)图不变;(b)图中不变)
13、当不变时，得图3(a)，第l条线以下的区域代表保证几乎必然指数稳定，第2条线和第3条线下的相应区域分别保证一阶矩和二阶矩指数稳定9。由图3（a）看出，一阶矩指数稳定的条件比几乎必然指数稳定的条件要强；同样，二阶矩指数稳定的条件比一阶矩指数稳定的条件要强，这与理论结果相符。图4几乎必然指数稳定、一阶矩指数稳定及二阶矩指数稳定的区域(不变)比较图3和4，得到：不论在何种情况下，保证几乎必然指数稳定的区域大于一阶矩指数稳定的区域；保证一阶矩指数稳定的区域大于二阶矩指数稳定的区域，这个结论与理论结果是一致的。接下来讨论能使埃博拉病毒最终灭亡的条件。假定其他参数都不变，定义区间为，以月为单位，根据数值解
14、方法进行模拟，得图5(a)。由图5(a)可知，感染者比例会迅速增长并最终趋于1，也就是说感染者比例会在未来若干月后全面爆发，这个结果是非常可怕的。为了防止埃博拉病毒在未来爆发，必须采取相应的措施。由稳定性条件，在假定其他参数不变的前提下，可计算出的一个临界值：1.8596。然后验证这个临界值对病毒传播趋势的作用。假定其他参数均不变，时间区间，以月为单位。由于时间区间加大，为更加清晰地显示变化趋势，取节点数为100,得图5(b)。（a） （b）图5 未来感染者比例由图看出，未来感染者人数比例不会像图5(a)中那样迅速增长直至爆发。相比而言，感染者比例增长要缓慢得多，从而在一定程度上缓解了埃博拉病
15、毒爆发的风险，对现实具有重要的指导意义。但是，由于取的是临界值，并且从图中也无法看到感染者比例明显的下降趋势，所以为了达到更好的效果，进一步，在假定其它参数不变的情况下，再降低的值，来看感染者比例变化的情况。令，变化情况如下-新加的 （a） （b）由上面的分析可知，当将疾病的传染率控制指标降低到一定的临界值之下时(假定在其它参数不变的情况下)，埃博拉病毒的传播速度就能得到降低，并且在一定的时间区间内最终达到感染者比例趋于0，也就是灭亡了。仿真结果对得到的结论进行了验证，指出埃博拉病毒最终消亡的条件。同时，我们也看到了降低接触率的作用，当降传染率控制指标控制到临界值以下后，再进一步控制接触率使其
16、下降，在使感染者比例最终趋于零后，能够达到使感染者比例下降速度更快的效果。由此得出结论，一方面，我们必须采取相应措施降低疾病传播率到其临界值以下，这样可确保埃博拉病毒不在未来爆，,另一方面，也要采取措施降低相应的接触率，这样才能降低其传播速度，减小其传播直至爆发的风险。3.4结果分析该模型建立的随机微分方程，研究的是埃博拉病毒感染者人数占总人口比例的变化趋势，通过模型表征研究疾病的发展趋势、动态变化，便于对疫情进行监控，提出一定的策略来降低感染者人数比例。4时间优化模型埃博拉病毒是一种急性传染病，患者会在数天内死亡，而且传染率极高，同时死亡率高达，所以把时间控制在最短是消除埃博拉病毒的关键性因
17、素。4.1符号说明符号符号说明药物运达需求地点时的感染者人数培育一批药物或疫苗的时间培育药物或疫苗的速度从到运输药物的成本从到的运量4.2模型假设l 假设每位感染者的用药量均为一剂量。虽然目前已经研制出应对埃博拉病毒的疫苗或者药物，但是不同感染程度的患者所需实际的药剂量数据不易获得。l 假设疫苗与药物的生产地，开始培育的时间以及生产速度均相同，且培养药物和疫苗的周期均为，但两者的作用对象不同，疫苗作于与健康人群，而药物作用于感染者。l 假设各地生产药物的速度相同，每批生产的药量满足当前三个国家的需求量。l 假设运输问题中各个产地的产量相同。用来运送药物的飞机类型相同，且保持相同速度行进。4.3
18、模型建立与求解在药物或疫苗运输方面，我们通过建立线性规划模型，在满足各需求地需求量的前提下，制定相应调运方案，将这些物资运到各个需求地，使总运费最小，进而得出运输用时最短的调运方案。考虑到成本问题，我们仅研究供求相等的情况。已知有个生产地点可供应药物，其供应量分别为，由个需求地，其需求量分别为，从到运输药物的单位时间成本为，见表1： 表1 药物运输时间成本表产地需求地产量需求量建立如下数学模型：满足设，分别代表几内亚、利比里亚以及塞拉利昂；，分别代表能够生产药物的美国、中国、日本、俄罗斯、法国以及瑞士。根据模型一，我们预测出各国感染者人数，如图6。如果我们预期在2015年末所有非晚期感染者均被
19、注射药物，根据假设条件可确定各个国家药物需求量，即，那么。图6 三个国家预测感染者人数通过产地与需求地的距离，计算出时间成本，得到表2：表2 运输药物时间成本表产地需求地产量几内亚利比里亚塞拉利昂美国10.010.110.213512中国13.313.31313512日本15.415.815.813512俄国7.27.67.513512法国4.85.35.113512瑞士4.75.35.113512需求量21146234543647281072建立如下数学模型：满足在模型一中，未考虑药物治疗和疫苗预防，感染者的传播特性都是一样的。（1）当感染者接受药物治疗后，设传染率为治疗前的倍，则此时传染率
20、为，结合模型一中仿真系统的可得图7。从图7中可得感染者接受药物治疗后，感染者人数比例曲线趋于平缓，可见接受药物治疗对减小埃博拉病毒的压力有一定作用。（2）当非感染者注射疫苗后，就确保了没有新感染者，就切断了埃博拉病毒的传播途径，可得图8。（3）同时考虑药物治疗和疫苗预防的情况，综合（1）（2）情况，得出图9。4.4结果分析在满足各需求地需求量的前提下，利用Lingo软件确定出各个产地运往不同需求地的药物量，然后得到最优调运方案。根据表5.2.2中的运输时间筛选出本方案的最短运输时间为15.8小时。5最优隔离控制5.1符号说明符号符号说明被隔离患者人数对染病者实施隔离控制的比例未被隔离染病者死亡
21、率被隔离的病人治愈率易感者人数染病者人数 病愈者人数 出生率=死亡率 不易感病者输入比率 没有被隔离治疗者治愈率 环境中的总人口数5.2模型建立与求解控制传染病最好的方法是隔离控制,下图为加入被隔离染病者舱室的舱室图： 图3.3-1加入隔离后的舱室转移图通过模型的稳定性分析，我们知道要想消除疾病必须得隔断染病者输入，这也是当传染病暴发时，政府控制染病者流动的原因，所以我们假设在人口的输入输出中没有染病者。于是可以建立如下的微分代数系统：式中：表示易感者，染病者和治愈者的输入率。系统的初始值给定。5.3最优控制系统考虑到如何实施控制方案才能保证疾病大规模爆发的同时，我们付出的各种代价（包括政治代
22、价、政策代价、经济代价等）达到最小，我们采用这样的性能指标，其中，为末端时刻，此性能指标考虑了染病者人数和实际隔离控制所付出的代价。事实上，我们在整个时间段上采用最大强度的隔离控制必然会使染病者数量最小，然而任何控制都是有代价的，所以实施隔离控制的代价应该在性能指标中考虑到，性能指标中的和分别是相应的权重，表示对应的代价的重要程度。我们的目标是寻找最优控制使得其中，是关于的函数，控制约束集合，为指定常数。对于初值已给定的系统，给定目标函数：存在一个最优控制，使得5.4最优隔离控制极值原理给出了最优控制的必要条件，将最优控制问题转化为最小化的哈密顿函数问题，状态方程的不同以及终端条件的不同，都会
23、使得横截条件发生变化.下面应用极值原理求解最优控制律。将本文所讨论的最优控制问题转化为最小化哈密顿函数问题：设，其中是第个系统状态方程的右端，是协态变量，如果是上面系统的最优控制，是相应的状态解，那么存在协态变量使得：协态方程： 横截条件和最优控制5.5数值仿真数值仿真研究埃博拉爆发期间的最优检疫隔离控制策略。根据几内亚埃博拉期间的具体情况，取系统控制变量的最大值为0.50，表示平均至少需要两天才能够有效的隔离染病者。3.3-1 系统参数对应表参数值参数值0.230.030.970.030.0000340.03370.00790.00680.03863001121971579158637103
24、68922131341末端时刻630天（21个月）控制上界0.5控制下界0.05根据系统参数值运用MATLAB绘制易感者人数，染病者人数和治愈者人数对于时间的图像曲线，根据实际情况与数据显示，选出最优隔离控制的值。首先考虑未使用隔离和使用隔离的染病人数对比（见图3.3-2），得到结论：使用隔离会在一定程度上降低染病人数。接下来考虑使用隔离的情况。对于易感病人数对应时间的函数图像（见附录4），理想的图像应该显示出以下信息：首先在易感者初始值的基础上，曲线下降速率越快就越理想，当然最后不能保证所有人都不易感病，所以在易感病人口数这组图里，时是最接近理想曲线的，所以是这组数中的最优解。同理，在染病者
25、人数图像组里，首先峰值最小，而且大部分呈下降趋势，在这一组里，又是最优控制解。再看治愈者人数，治愈人数肯定为上升趋势，而且上升速度越快就越理想，同样还是为最优控制解（见图3.3-3）。(a) (b)Fig .3.3.2 (a)图为未实施隔离控制下染病者人数变化情况，(b)图为最优隔离控制下染病者人数变化图3.3-3（a） 在最优控制下，染病者的人员总数情况图3.3-3（b） 在最优控制下，易感者的人员总数情况图3.3-3（c） 在最优控制下，治愈者的人员总数情况 4.模型评价4.1模型优点模型一考虑了实际情况下环境因素的干扰，利用随机微分方程模型对埃博拉病毒的传播过程进行建模，结合实际数据，运
26、用随机微分方程数值解方法，对未来感染者人数占总人口的比例进行预测，在此基础上提出通过调节传染率控制指标来达到降低未来感染者比例的目的，为更进一步地对埃博拉病毒传播过程进行精确地建模提供了理论上的基础，为分析其流行趋势提供了更为有利的工具。l模型二主要从优化时间的角度研究药物对埃博拉病毒的传播过程的影响，将问题简化，使之通俗易懂。模型三是在模型一、二的基础上，提出一种控制疾病继续传播的有效措施。本文提出利用极值原理给出了最优的隔离策略，明确出以什么样的强度进行隔离能够在保证效果最好的同时，又能尽可能的节约成本。仿真结果验证了最优隔离控制的有效性，并且说明了在疾病大规模爆发时，釆用最大强度的隔离措
27、施对控制疾病是非常有意义的。4.2模型不足在药物运输过程中，没有全面考虑国与国之间的关系程度，所以药物的分配会有些差距。各国运输药物使用同款运输机，并且运行速度保持一致，实践起来，有一定难度。5.模型推广本文主要解决了传染病的传播、药物的运输、分配和优化问题。文中所涉及的方法和模型可以解决类似的病情，比如SARS、禽流感、疯牛病、艾滋病等。本文的模型还可以用来解决货物的运输或采购、企业绩效的预测、大自然的绿化优化等问题。6.参考文献1 Shanglai Guo. Stochastic ControlM. Beijing: Science Press, 1998.2 M.A.Kouritzin,
28、 L.Deli. One x Plieit solution to stoehastie differential equations. Stoehastl.eAnalrsisandAPPlicatl.ons,2000,18(4):571-580.3 Qiyuan Jiang, Jinxing Xie, Jun Ye. Mathematical modelM. Beijing: Higher Education Press, 2003.4 Min Xu. Stochastic model machine risk analysis of HIV transmissionD. Donghua U
29、niversity. 2007.12. 5 http:/zh.wikipedia.org/wiki/ 2015.2.96 http:/www.worldbank.org/2015.2.97 Heinz F, Thomas W G. Ebola haemorrhagic fever J. The Lancet, 2011,377: 849-862.8 Maosheng Liu, Yicang Zhou. Research of one type age structure of HIV transmission dynamics model. North China Institute of T
30、echnology, 2004, 25(2): 25-30.7.附录Appendix1.基于Lingo的产销平衡的数值解法：min 10.0x11 + 10.1x12 + 10.2x13 +13.3x21+15.4x31+15.8x32+15.8x33+7.2x41+7.6x42+7.5x43+4.8x51+5.3x52+5.1x53+4.7x61+5.3x62+5.1x63s.t.x11 + x12 + x13=13512x21 + x22 + x23 =13512x31 + x32 + x33=13512x41 + x42 + x43=13512x51 + x52 + x53=13512x
31、61 + x62 + x63=13512x11 + x21 + x31 +x41+x51+x61= 21146x12 + x22 + x32+x42+x52+x62 = 23454x13 + x23 + x33 +x43+x53+x63= 36472end模型的最优解见下图：2. 埃博拉病毒患病人数随时间的分布情况Tab1.The Cumulative MonthsPeriod of TimeCasesThe Late Patients122th May- 21th Apr21579222th Apr- 21th May27089322th May- 21th June599191422th
32、June- 21th July1093433522th July- 21th Aug26151188622th Aug- 21th Sept62633363722th Sept- 21th Oct99465087822th Oct- 21th Nov1535110246922th Nov- 21th Dec19497119091022th Dec- 25th Jan22092132823.基于Euler迭代法下的随机微分方程的%The eqaution that approximatively solving stochastic by Euler iterative.function mec
33、ca_zjM=500;N=1000;zMesh=repmat(0,M N);Max=50;%Parameters and innitial values ahout the model.sigma=1;beita=1/5;lamga=1/(4*(sigma2);X0=-20;-30;derta=0.1;dertat=lamga*(derta2);timef=40;stepMax=ceil(timef/dertat);for i=1:Maxstep=1;% Starting in the process of Euler iterative. X=X0;while(step=stepMax) %
34、 X=pointP(:,step);% temp=(plan(step*dertat)+affine(0;0,step*dertat,0)*X); temp=(plan(step*dertat)+affine(X,step*dertat,1); ekq=temp(1)/(2*(sigma2)*(1-hfunction(X,beita); nkq=temp(2)/(2*(sigma2)*(1-hfunction(X,beita); xkq=1/(lamga*(sigma2)*(1-hfunction(X,beita)-4; ckq=2*cosh(-derta*ekq)+2*cosh(-derta
35、*nkq)+xkq; pkl=exp(-derta*ekq)/ckq; pkr=exp(derta*ekq)/ckq; pkd=exp(-derta*nkq)/ckq; pku=exp(derta*nkq)/ckq; pko=xkq/ckq; seed=rand(1); if(seedpkl) X=X+-derta;0; elseif(seedpkl+pkr) X=X+derta;0; elseif(seedpkl+pkr+pkd) X=X+0;-derta; elseif(seedpkl+pkr+pkd+pku) X=X+0;derta; else X=X; end J=round(N*(X(1)+30)/90); I=round(M*(X(2)+50)/80); if(IM) I=M; end if(JN) J=N; end zMesh(I,J)=zMesh(I,J)+1; step=step+1;endiendtemp=max(max(zMesh);zMesh=255*zMesh/temp;zMesh=uint8(zMesh);imshow(zMesh);imwrite(zMesh,1.bmp,bmp);4.图表4.1感染病人总数关于时间的图组