有关SARS传染病的数学预测模型

1、有关SARS传染病的数学预测模型摘 要本文针对问题一，首先从附件1所给模型参数选取的合理性和科学性入手，分析了K和L的价值作用，并结合模型的实际预测结果，对模型的实用性和合理性进行了评价。同时，根据SARS的传播特点，指出了该模型的不足之处。针对问题二，在克服前模型不足的前提下，把人群划分为五大类，建立了SARS传染病动力学预测方程，并用遗传算法对所给参数进行估计，最后利用龙格库塔数值积分方法分别做出了这五类人群变化的趋势线，与实际情况的变化相吻合，并根据题意做出了评述。针对问题三，通过1997年到2003年8月北京海外游客的数据，就非典对旅游业产生的影响进行了分析和预测。首先不考虑非典的影响

2、，即不考虑2003年4-8月份的数据的情况下，利用神经网络和GM（1，1）模型法分别进行预测，再结合标准差法确定组合权重实现组合预测，得出4-12月份的结果为下：28.9204、30.3630、30.1892、28.7201、31.5473、30.4872、31.2696、29.5585、25.7050。其次在有非典影响的情况下，引入心理影响因子收缩因子，将非典对旅游业的负面影响用收缩因子进行描述，根据4-8月份的数据用最小二乘法估计收缩因子的参数，从而得到9-12份的预测因子，最后结合在不考虑非典影响情况下得到的预测数据，便得到了9-12份受非典影响后的预测数据，结果为22.7059，24.

3、0840，23.3815，20.7795。 最后，根据传染病模型的特点和作用，提出了建立数学模型对疫情分析、预测控制方面的重要意义。关键词： 龙格库塔 神经网络 GM（1，1）模型 组合预测模型 传染病动力学模型 遗传算法一、问题的提出SARS（Severe Acute Respiratory Syndrome，严重急性呼吸道综合症, 俗称：非典型肺炎）是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，我们从中得到了许多重要的经验和教训，认识到定量地研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。就SARS 的传播建立数学模型

4、，具体要求如下：（1）对附件1所提供的一个早期的模型，评价其合理性和实用性。（2）建立模型，说明为什么优于附件1中的模型；特别要说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，这样做的困难在哪里？对于卫生部门所采取的措施做出评论，如：提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计。（3）收集SARS对经济某个方面影响的数据，建立相应的数学模型并进行预测。（4）给当地报刊写一篇通俗短文，说明建立传染病数学模型的重要性。 二、模型一的评价附件所给的模型能很好的反映SARS早期的传染规律，主要体现在以下几点：第一、 数L的选取具有科学背景，据中科院的报

5、告，SARS的潜伏期一般为23周，因此将L选为20天作为一个周期将有效的避免数据的波动性，并且在L天内的数据能很好的吻合指数规律。第二、 数K的选取使模型有较强的适应性，对不同地区不同时段的数据，都能通过参数K来调整模型，并且参数K很好的反映了某种社会环境下，一个病人传染其他人的平均概率，并充分显示了全社会的警觉程度，政府和公众的措施以及医疗卫生条件的好坏。第三、 对实际情况的分析具有较高的准确性，这表现在以下两个方面：从香港和广州K值回落速度的快慢，反映了香港的医疗卫生条件比广州好，这与实际情况相符；从香港、广东短期内K值调整的幅度，反映了政府和社会对非典的高度关注，这和实际情况也是统一的。

6、综上所述，该模型有很好的合理性和实用性。但它的缺点也是显而易见的，过分依赖对数据的统计，未考虑将人群分类和对病人的隔离，使它对后期的K值难于确定，导致预测结果过于粗糙。三、模型二的分析与建立1. 问题分析与假设：SARS不同于一般传染病模型，它有着自身的传播特点。在SARS传播的初期，政府和社会对它的传播速度和危害程度认识不够，没有采取足够的控制措施，使得SARS病毒传播的速度相当快。当非典感染者的数量不断增加，政府开始采取各种措施对病毒进行控制，使得非典感染者的数量逐渐减少。所以SARS病毒的传播可以分为三个阶段，分别为：控制前期, 没有采取任何控制措施的阶段。过渡时期，从SARS被人们重视

7、到政府采取控强有效制措施前的一段时期。控制后期,在政府开始控制后的时期。我们所要做的是对SARS在北京传播的情况进行预测。北京是在SARS刚刚大肆传播就采取很强有力的措施，因此，北京的过渡期可以包括控后期；我们将北京的SARS传播规律用“控制前”和“控制后”两个时期来模拟。由分析问题我们可以知道，控制前和控制后，SARS的传播源有很大的差别。控制前，每个病人都将成为传染源，而控制后只有游离的带菌者成为有效传染载体。因此，我们将控制前和控制后进行分段处理，分别建立模型。由此，我们可以对问题做如下假设：基本假设：1. 只考虑患SARS的病人，患其它病的病人归为健康者。2. 不考虑隐性患者，即只要感

8、染上SARS病毒最后都会表现症状。3. 由于SARS的传播时间相对较短，所以不考虑这段时间内的人口出生率和自然死亡率。4. 卫生部公布的数据正确可靠。5. SARS患者康复后具有免疫力，退出传染系统。控制前的传播模型的相关假设：1. 在SARS传播期内北京地区的总人数N视为常数，即不考虑人口的流入流出（流入流出人口与总人口相比可忽略不计）2. 将SARS所有可能的传染途径视为与传染源的直接接触。3. 将人群分为四类：健康者：用表示健康者的总人数；病人 ：用表示病人的总人数；移出者：包括“被治愈者”和“死亡者”两种情况，用表示移出者的总人数；处于潜伏期者：这些人还没发病，但他们最终将发病，用Q表

9、示处于潜伏期者的总人数。控制后的传播模型的相关假设：1. 不考虑被隔离但未感染病毒的情况，因为这部分人由于被隔离，对病毒的传播不产生影响。2. 被隔离人群完全断绝与外界的接触，不具有传染性。3. 将人群分为五类：健康者：用表示健康者的总人数；病人 ：用表示病人的总人数；移出者 ：用表示移出者的总人数；疑似病人 ：包括已出现有关症状但未确诊的被隔离者和还未出现症状但已疑为带菌者而被隔离观察的人群，用X表示疑似病人的总人数；游离带菌者： 没有被隔离的病毒携带者，用Z表示游离带菌者的总人数。2. 模型的建立：1) 控制前模型的建立：符号说明：每个病人单位时间内有效接触每个健康者的概率。：退出率，为S

10、ARS患者的日死亡率和日治愈率之和。：处于潜伏期的病人的日发病率。控制前模型方程的建立分析各个变量间的关系，结合一般传染病的传播规律，我们可以建立如下的动力学模型： (1)由于SARS控制前的相关数据无法查找，所以我们只列出模型的方程，而不做过多的分析，我们着重分析控制后的模型。2) 控制后模型的建立符号说明：疑似病人中每日被排除的人数占疑似人数的比例：疑似病人中每日确诊为病人的人数占疑似人数的比例：每个游离带菌者转化为病人的日转化率：每个游离带菌者发病后被收治前平均每天有效感染每个健康者的概率；：被游离带菌者有效感染的人中被隔离的人的比率控制后模型方程的建立同样，我们可以仿照控制前模型的建立

11、方法，列出控制后模型的动力学模型： (2)四、模型二的求解及说明在模型二的建立过程中，由于政府在控制前期没有采取有力的措施对疫情进行控制，所以相关的数据无法查找，无法对控制前的模型做很有意义的解析分析。因此，我们未对控制前模型进行求解。下面我们来对控制后的模型进行求解。很明显，从我们建的模型中无法求出精确的解析解，因此，我们采用龙格库塔方法来求模型的近似解。我们对参数进行估计，并根据实际情况对方程组进行合理的简化。考虑到健康者S对于病人数I来说是一个很大的数，的变化量很小，所以健康者S在一段时期内可以视为一个恒定的数量，故我们忽略（2）中的第一个方程，然后对简化后的模型进行求解。方程组的参数辩

12、识由国家公布数据，从4月21日算起， 由题中附录1可知每个游历带菌者平均每天能感染的人数为为讨论方便起见，引入如下记号：此时系统（2）等价于以下初值问题 其中设某地区第天的实际累计病例数为，实际出院人数为，实际死亡人数为，则该地区第天的实际染病人数为，移出人数为，为使方程组（2）更好地描述SARS传播规律，则需选择，使第天的理论值与实测值 误差尽可能小，这可表示为如下的参数辨识问题： （3）我们考虑用遗传算法来求解该问题，其因变量为，并且易知,我们采用实数编码，取遗传种群大小为n=5。取适应函数为目标函数J。遗传算法的流程图如下： Step1:随机产生十组初始染色体 Step2:代入约束条件中

13、用龙格库塔数值积分，得出数组I ,T Step3:计算适应度即目标函数J Step3:用轮盘赌方法繁殖新种群 Step4:对种群进行杂交、变异、自然选择 Step5:如果迭代次数小于规定次数，返回Step2，否则继续下一步 Step6:输出种群，取其最优者搜索结果：代入式（2）式用龙格库塔数值积分算法得下图。（具体程序见附录）图1模型说明与附件的模型相比，我们的模型考虑的因素更多，分析较为全面，能更准确的反映非典的传播规律，同时还能够反映病情传播的全过程，使预测的效果更好。另外，我们的模型对政府采取控制措施有明确的指导作用，要真正预测疫情的发展，必须有可靠的初始数据，并且对当地的人口分布，人口

14、流动有清楚的了解，以及对当地医疗机构水平，政府和公众的重视程度等要有明确的认识，从而使确定的参数更符合实际，但是这些条件无法量化，甚至有些是随机的，另外有些是保密的，所以这使得参数的确定非常困难。卫生部门首先要作好宣传工作，提高人们对疫情的认识。其次，对公共场所，要做好消毒工作；对流动人口及本地人口，进行严格监控，发现异常人口，立即确诊，并作好隔离工作。如果隔离时间延迟，从模型可知，控前期时间延长，必定使患病人数在疫情达到高潮时数量很大，从而在采取有效的措施后，患病人数降到某一数量将需要更长的时间。五、模型三的建立与求解SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，特别在旅游业

15、方面影响最大，所以在模型三中我们采用附件3所提供的数据，来分析预测2003年后四个月北京市接待海外旅游人数，同没有爆发SARS得到预测数据相比较，得出SARS对旅游业影响的程度。首先我们考虑用神经网络中线性网络模型和GM（1，1）模型来分别单项预测在没有爆发SARS的情况下，2003年后四个月北京市接待海外旅游人数，然后采用标准差法进行权重分配，建立组合预测模型得出最后的预测结果。单项预测模型的建立与求解：1. 神经网络模型我们将从1997年1月到2003年3月即SARS爆发以前的数据作为输入变量，在实现时，利用solvelin函数设计网络，可以得到网络的权值和阀值，然后利用仿真函数simul

16、in得到网络运算的结果，最后输出模型的预测结果。主要过程如下： time=1:75;t=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6.15.4 17.1 23.5;w,b=solvelin(p,t)a=simulin(p,w,b)（具体的实现程序见附录）模型的预测结果为：表1 神经网络预测结果4月5月6月7月8月9月10月11月12月25.198727.494228.649828.352827.552127.823430.455628.809224.1310 图2图2为神经网络线形网络模型中原始数据和预测结果的比较图。分析

17、比较可以发现神经网络的预测结果和原始数据拟合的较好。2. GM（1，1）模型 考虑分为横向和纵向两种情况进行预测，所谓横向预测就是将1997年1月到2003年3月即SARS爆发以前的数据作为一个序列来预测后9个月的情况；纵向分别利用1997-2002年各个月份的情况来预测2003年的情况。（1）建立横向预测模型，其主要过程如下：原始序列 采用1997年1月到2003年3月即SARS爆发以前的数据作为原始序列代入数据有 累加生成序列 对原始序列做累加生成，得到新的序列，其元素为 所以 建立GM（1，1）模型的微分方程： 令，以累加生成序列，采用最小二乘法求得参数的估计值，即。所以由此建立的GM（

18、1，1）模型白化方程为： 求解上述微分方程，得GM（1，1）预测模型为： 模型的结果表2 横向GM（1，1）预测结果4月5月6月7月8月9月10月11月12月31.992133.950432.396329.115835.419133.751333.291230.730925.0645图3模型检验由计算可得（检验的内容和计算方法见附录），模型的关联度为0.3247，方差比为0.3640，小误差概率等于1。查附表1有，该模型的方差比大于0.35，可知该模型的预测较好。（2）建立纵向预测模型的主要过程和横向预测相同，在这里我们就不多加描述，只将模型得到的重要数据和预测结果列出。表3 参数a. u的估

19、计结果a-0.0940-0.1029-0.1099-0.0842-0.0809-0.0684-0.0504-0.0624-0.0977u17.207917.184715.622416.554520.671921.469724.006820.587612.8222模型的结果表4 纵向GM（1，1）预测结果4月5月6月7月8月9月10月11月12月27.947128.111128.276128.442128.609128.777028.945929.115829.2867 3、组合预测模型的建立与求解：采用标准差法确定组合权重。设神经网络线形网络模型、横向GM（1，1）模型、纵向GM（1，1）模型

20、的预测误差的标准差分别为且，取；m为模型的个数。 计算各单项模型的权重分别为，根据这个组合权重，建立组合预测模型如下：其中，y为组合预测值，为神经网络预测值，为横向GM（1，1）预测值，为纵向GM（1，1）预测值。 运行附录中的程序，可以得到各单项模型4-12月份的权重见下表：表5 单项模型权重4月5月6月7月8月9月10月11月12月0.31220.32570.31940.31160.30160.33100.32500.35810.30000.23510.25420.24530.23440.22020.26170.25320.29980.21800.45270.42010.43540.454

21、00.47830.40730.42190.34210.4819 将各单项模型的权重代入组合预测模型，我们就可以得到组合预测的结果，见下表：表6 组合预测结果4月5月6月7月8月9月10月11月12月28.920430.363030.189228.720131.547330.487231.269629.558525.7050图4现在我们来分析非典对经济造成的影响。根据非典传播的规律，我们可以将非典对经济产生的影响分为三个阶段。第一阶段为非典防治阶段，在这个阶段人们特别恐惧，失去了消费信心；第二阶段与非典过渡时期，在这一阶段，虽然隔离制度取消了，非典病人也没有了，但人们还是后怕，心有余悸；第三阶段

22、为正常阶段。 由此可知4-5月应属于第一阶段，6-8月份应属于第二阶段。根据国家统计局估计，到年底非典的影响也将不能完全消除，所以我们认为9-12月份也属于第二阶段。 非典将使经济发展的速度减慢，因此我们定义Y为非典因子或收缩因子，显然是关于时间t的函数，即。 设为某天未受非典影响旅游的人数值； 为某天受非典影响旅游的人数值，则。 以下为几个引理：引理1 在防治阶段收缩因子随时间的增加而减小；引理2 在第二阶段收缩因子随时间的增加而增加；引理3 第二阶段收缩因子将随总长度的增长而减小。所以，我们可以定义如下函数： 对于4月份来说，我们无法准确得出非典爆发的初始时间，所以设4月份非典已传播的时间

23、T0，我们需要确定的参数有。 结合非典时期北京接待海外旅游实际人数和组合预测得到的结果，可以列出如下方程： 由上式易得 对于(1)、(2)、(3)、(4)中的后三个方程我们采用最小二乘方法拟和求解所以构造函数如下+ 调用fmincon函数求解，得利用参数可以求得得到最终预测值为表7 最终预测值9月10月11月12月22.705924.084023.381520.7795 图5 图5为有非典影响和无非典影响两种情况预测结果的对比图六、模型检验稳定性检验：由于模型2的方程是微分方程，因此须检验稳定性。故在最优解的邻域作微小的随机扰动，取几组新值 并对其数值积分，由图可得模型比较稳定。 （图见附录）

24、预测模型准确度的测试为了检验预测模型的准确度，我们考虑弃掉2002年4月到12月的实际数据，而用我们已建立的模型进行预测，得到一个理论值，然后用理论值与实际值的差异来描绘其准确度。理论值：27.24 28.27 27.82 24.5527 26.65 16.4501 27.3750 24.8229实际值：29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9易知除开16.4501点严重偏离实测值以外其他点均符合的非常好。因此，该模型具有很好的可信度。七、模型的推广与评价模型二不仅可以适用非典的分析和预测，还可以适用于其它具有相似特征的传染病的情况。在模型中，我们应用遗传

25、算法的模型进行对参数进行估计，得到的预测结果与实际情况吻合较好。在模型三中我们采用标准差法确定组合权重，使组合预测模型能够应用于其它行业，具有较强的实用性。当然，模型也有一定的不足，如在模型三中的收缩因子具有一定的随意性，这可能导致模型预测的结果产生较大的误差。模型二主要针对北京地区的情况，将模型用于其它地方可能产生一定的误差。八、短文SARS与数学模型记得在非典高峰期，全国上下人心惶惶。有人说非典要持续一年两年，或者说非典不可战胜，也有人说药到自然除，不需要隔离等等，甚至有人为了避免被隔离而殴打警察等等。这时，为了让人们认识到SARS的传染规律，积极配合国家的非典预防政策。科学家们充分利用传

26、染病规律和已拥有的病人数据资料，建立了SARS的传染病数学模型。该模型明确指出SARS并不可怕，高潮期过后病人数会逐渐减少直到为零，并根据过去的数据预测了现有每天新增病例人数，由于理论预测与实际吻合很好，使得人们逐渐相信该模型的预测能力，相信SARS很快过去。另外，该模型也指出了隔离病人和疑似病人能有效地减少传染源，控制SARS的疫情扩张，这使得对政府采取的各种预防控制措施具有积极意义。同时，数学模型不仅能预测传染病病人人数的涨落趋势，还能反映非典对经济各方面的影响。面对非典带来的负面影响，需要较为精确的经济发展态势，以正确地引导政府的宏观调控和相关行业的公司部门经营业务的调整，来尽最大努力减

27、少非典带来的经济损失。总之，从各个方面来说，在SARS这场战争中，人们应用数学模型这一法宝解决了许多实际问题。参考文献：1 杨建刚，人工神经网络，浙江大学出版社，20002 MATLAB6.5 辅助神经网络分析与设计，电子工业出版社，20033卢奇等，组合预测模型在我国能源消费系统中的建构及应用，系统工程理论与实践，第23卷，24-28页，2003附录：一、GM（1，1）的关联度分析：关联度说明了原始序列和生成序列之间的关性程度，用来描述模型模拟值序列对原始序列值的拟合的程度。 其中：的关联系数 是分辨率，一般取50%，即。当时，关联度大于0.6，模型的拟合精度就能达到比较满意的程度。后验差分

28、析：计算均方差比：计算小误差概率：附表1 GM（1，1）模型精度等级标准PC模型等级标准1级（好）2级（合格）3级（勉强合格）4级（不合格）二、模型检验中得到的图象附图1 附图2附图3三、模型中的主要原程序：神经网络模型：clf;figure(gcf)%setfsize(600,250);echo onclcpauseclctime=1:63;t=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6. 9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9. 10.1

29、12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5. 11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5. 11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7. 13.7 29.7 23.1; m1=12;p=delaysig(t,1,m1)pauseclcplot(time,t);alabel(Time,Target Signal,Signal to be Predicted) pauseclcw,b

30、=solvelin(p,t)pauseclca=simulin(p,w,b)plot(time,a,time,t,red)alabel(Time,Outut_Target +,Output and Target Signals)pausee=t-a;plot(time,e);hold onplot(min(time) max(time),0 0,:r)hold offalabel(Time,Error,Error Signal)echo offfor i=1:9 for j=1:m1 q(j,1)=t(51+i+j-1); end a(75+i)=simulin(q,w,b); t(75+i)

31、=simulin(q,w,b);endplot(a)hold onplot(t,red)s=std(e)a遗传算法：function =qp2(K)%编码tica=0,0,0,0,0,0;%下界b=1,1,1,1,1,0;%上界n=15;%种群数l,m=size(a);%始种群选for i=1:n for j=1:m r(i,j)=unifrnd(a(j),b(j),1,1); endendglobal R;for i=1:n R=r(i,:); ts=1:38; x0=514 402 666 51; t,x=ode45(jifen,ts,x0); XX(i)=X;endclear Rclea

32、r xclear t%for i=1:K for J=1:m fit=fitness(r,n);%求适应值 r0=r; fit0=fit; r=fanzhi(r,fit,n,J);%繁殖 r=jiaocha(r,n,J);%交叉 r=bianyi(r,n,J,a,b);%变异 fit=fitness(r,n); r=xuanze(r,fit,r0,fit0,n,J,a,b);%自然选择 endendfit=fitness(r,n);disp(最大值);fit,l=max(fit);f=fun(r(l,:);disp(f);disp(最优解);disp(r(l,:);toc%产生适应函数 fun

33、ction fit=fitness(r,n)for i=1:n fit(i)=fun(r(i,:);endmax1=max(fit);min1=min(fit);for i=1:n fit(i)=(fit(i)-min1)/(max1-min1);end %繁殖种群function r=fanzhi(r,fit,n,J)p=zeros(1,n);s=0;for i=1:n s=s+fit(i);endif s=0 fit=fit/s;endfor i=1:n if i=1 p(i)=p(i); else p(i)=p(i)+p(i-1); end end for j=1:n R=rand(1,

34、1); for i=1:n if i=1 if R=0 r(j,J)=r(1,J); end end if i1 if R=p(i-1) r(j,J)=r(i,J); end end end end%杂交function r=jiaocha(r,n,J)for i=1:n c1=unidrnd(n); c2=unidrnd(n); t=r(c1,:); r(c1,:)=r(c2,:); r(c2,:)=t;endfor i=1:2:n-1 a=rand(1,1); x=a*r(i,J)+(1-a)*r(i+1,J); y=a*r(i+1,J)+(1-a)*r(i,J); r(i,J)=x; r

35、(i+1,J)=y;end%变异function r=bianyi(r,n,J,a,b) for j=1:n p=0.01; t=rand(1,1); if tp r(j,J)=a(J)+(b(J)-a(J)*rand(1,1); end end %自然选择 function r=xuanze(r,fit,r0,fit0,n,J,a,b)c,d=sort(fit0);l,n1=size(d);a1,b1=sort(fit);%b记录位置r(b1(1),:)=r0(d(n1),:);%删出操作以防局部收敛s=3;for i=1:n x=0; for k=1:n if r(i,:)=r(k,:)

36、x=x+1; end end if x=n-s for j=1:m r(i,J)=unifrnd(a(J),b(J),1,1); end endendfor i=1:2*n N1=unidrnd(n); N2=unidrnd(n); Q=r(N1,:); r(N1,:)=r(N2,:); r(N2,:)=Q;endfunction f=fun(XX)N=339 482 588 693 774 877 988 1114 1199 1347 1440 1553 1636 1741 1803 1897 1960 2049 2136 2177 2227 2265 2304 2347 2370 2388

37、 2405 2420 2434 2437 2444 2444 2456 2465 2490 2499 2504 2512;%实际累计病例数D=28 35 39 42 48 56 59 66 75 82 91 96 100 103 107 110 112 114 116 120 129 134 139 140 141 145 147 150 154 156 158 160 163 167 168 172 175 181;%实际死亡人数C=46 55 64 73 76 78 78 83 96 109 115 118 121 134 41 152 168 175 186 203 244 252 25

38、7 273 307 332 349 395 447 528 582 667 704 747 828866 928 1006 1087;%实际出院人数x11=N-C-D;X21=C+D;s=0;for j=1:38 t,x=il32(aa,bb,c,y1,y2,q); x1=XX(:,1); x2=XX(:,2); s=s+(x1(j)-x11(j)2+(x2(j)-x21(j)2);endf=s;function y=jifen(t,x)A=R(1);B=R(2);c=R(3);y1=R(4);y2=R(5);q=R(6);%a=0.12;b=0.1393;c=0.7;y1=0.02;y2=0

39、.02222;q=0.023431;y=A*x(3)+x(2)*y2-q*x(1),-x(2)*y1-x(2)*y2+B*x(3)*c,B*x(3)*(1-c)-A*x(3),x(1)*q;GM（1，1）模型：A=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6. 9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9. 10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5. 11.4 26.0 19

40、.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5. 11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7. 13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9. 15.4 17.1 23.5 ;Aa=9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6 9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20

41、.1 15.9 10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5 11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5 11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7 13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9;x=zeros(1,74);%生成时间序列for i=1:74 for j=1:i x(i)=x(i)+A(j); endendfor i=1:73%B B(i,1)=-1/2*(x(i)+x(i+1); B(i