数学建模算法大全现代优化算法简介

1、第二十三章现代优化算法简介§ 1现代优化算法简介现代优化算法是80年代初兴起的启发式算法。这些算法包括禁忌搜索(tabusearch),模拟退火(simulated annealing),遗传算法(genetic algorithms ),人工神经网 络(neural networks )o它们主要用于解决大量的实际应用问题。目前，这些算法在理论 和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标一求NP-hard组合优化问题的全局最优解。虽然有这些目标， 但NP-hard理论限制它们只能以启发式的算法去求解实际问题。启发式算法包含的算法很多，例如解决复杂

2、优化问题的蚁群算法(Ant ColonyAlgorithms ) o有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限 制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。现代优化算法解决组合优化问题，如TSP (Traveling Salesman Problem)问题，QAP(Quadratic Assignment Problem )问题，JSP ( Job-shop Scheduling Problem )问题等效 果很好。本章我们只介绍模拟退火算法，初步介绍一下蚁群算法，其它优化算法可以参看 相关的参考资料。§ 2模拟退火算法2.1算法简介模拟退火算法得

3、益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温(这个过 程被称为退火)，粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形 成处于低能状态的晶体。如果用粒子的能量定义材料的状态，Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态 i之下的能量为E(i),那么材料在温度 T时从状态i进入状 态j就遵循如下规律：(1)如果E( j) E (i),接受该状态被转换。(2)如果E(j) E(i),则状态转换以如下概率

4、被接受：E(i) E(j)e KT其中K是物理学中的波尔兹曼常数，T是材料温度。在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处 于状态i的概率满足波尔兹曼分布：E(i)e KTPT(x i)一红e KTj SS表示状态空间集合。其中x表示材料当前状态的随机变量, 显然E(i)，.一 e 41lim KTT 0管 1se其中| S |表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，E(i) EmineKTlim E ET 0E(j)Emine KTj SE(i) Emine KTTim0E(j) Emine KTj SminlimT 0 e

5、j SminE(i) Emine KTE(j) EminKTE(j) Emine KTj Smin|Smin0其它Smini |E(i) Emin。其中 Emin minE(j)且 Smin j S上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思 想应用于优化问题就可以得到模拟退火寻优方法。考虑这样一个组合优化问题：优化函数为F : x R ,其中x S,它表示优化问题的一个可行解，R y | y R,y 0 , S表示函数的定义域。N(x) S表示x 的一个邻域集合。首先给定一个初始温度 To和该优化问题的一个

6、初始解x(0),并由x(0)生成下一个解x' N(x(0),是否接受x'作为一个新解x(1)依赖于下面概率：1若 f(x')f(x(0)P(x(0)x') f(x)")e T0其它换句话说，如果生成的解x'的函数值比前一个解的函数值更小，则接受 x(1)x'作为f(x') f(x(0)一个新解。否则以概率 eT0接受x'作为一个新解。泛泛地说，对于某一个温度 和该优化问题的一个解 x(k),可以生成x'o接受x'作为下一个新解x(k 1)的概率为：1若 f(x')f(x(k)(1)P(x(k) x

7、')f(x') f(x(k)eT0其它在温度Ti下，经过很多次的转移之后，降低温度 工，得到Ti 11。在T 1下重复上述 过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。我们注意到，在每个 Ti下，所得到的一个新状态x(k 1)完全依赖于前一个状态x(k),可以和前面的状态 x(0), ,x(k 1)无关，因此这是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态 x(k)生成x'的概率，在N(x(k)中是均匀分布的，且新状态x'被接受的概率满足式(1),那么经过有限次的转换，在温

8、度 Ti下的平衡态xi的分布由下式给出：f (Xi)当温度T降为0时,e TP0)一叵 e Tij SXi的分布为：*Pi若XiSmin其它并且 * P 1 Xi Smin这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率 1被找到。因此可以说模拟退火算法可以找 到全局最优解。在模拟退火算法中应注意以下问题：(1)理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算 机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢， 相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局

9、最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。(2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值Te,或连续几个温度下转换过程没有使状态发生变化算法就结束。(3)选择初始温度和确定某个可行解的邻域的方法也要恰当。2.2应用举例例 已知敌方100个目标的经度、纬度如下：经度纬度经度纬度经度纬度经度纬度53.712115.304651.17580.032246.325328.275330.33136.9348 156.543221.418810.819816.

10、252922.789123.104510.158412.481920.105015.45621.94510.205726.495122.122131.48478.964026.241818.176044.035613.540128.983625.987938.472220.173128.269429.001132.19105.869936.486329.72840.971828.1477 18.958624.663516.561823.614310.559715.117850.211110.29448.15199.532522.107518.55690.121518.872648.207716

11、.888931.949917.63090.77320.465647.413423.778341.86713.5667 143.54743.906153.352426.725630.816513.459527.71335.070623.92227.630651.961222.851112.793815.73074.95688.366921.505124.090915.254827.21116.20705.144249.243016.7044 117.116820.035434.168822.75719.44023.920011.581214.567752.11810.40889.555911.4

12、21924.45096.563426.721328.566737.584816.847435.66199.933324.46543.16440.77756.957614.470313.636819.866015.12243.16164.242818.524514.3598 58.684927.148539.516816.937156.508913.709052.521115.795738.43008.464851.818123.01598.998323.644050.115623.781613.79091.951034.057423.396023.06248.431919.98575.7902

13、40.880114.297858.828914.522918.66356.743652.842327.288039.949429.511447.509924.066410.112127.266228.781227.66598.083127.67059.155614.130453.79890.219933.64900.39801.349616.835949.98166.082819.363517.662236.954523.026515.732019.569711.511817.388444.039816.263539.713928.42036.990923.180438.339219.9950

14、24.654319.605736.998024.39924.15913.185340.140020.303023.98769.403041.108427.7149我方有一个基地，经度和纬度为(70,40)。假设我方飞机的速度为 1000公里/小时。 我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。 在敌方每一目标点的侦察时间不计， 求该架飞机所花费的时间 (假设我方飞机巡航时间可以充分长).这是一个旅行商问题。我们依次给基地编号为1,敌方目标依次编号为 2, 3,101,最后我方基地再重复编号为 102(这样便于程序中计算)。距离矩阵D (% )102102： 其中dij表示表

15、示i, j两点的距离，i, j 1,2, ,102 ,这里D为实对称矩阵。则问题是 求一个从点1出发，走遍所有中间点，到达点 102的一个最短路径。上面问题中给定的是地理坐标(经度和纬度)，我们必须求两点间的实际距离。设A,B两点的地理坐标分别为 (X1, y1)，(X2,y2),过A, B两点的大圆的劣弧长即为两点 的实际距离。以地心为坐标原点 O,以赤道平面为XOY平面，以0度经线圈所在的平 面为XOZ平面建立三维直角坐标系。则A, B两点的直角坐标分别为：A( Rc0sxi cos y1, Rsin x1 cos y1, Rsin y1)B(Rcosx2 cosy2, RsinX2 co

16、sy2, Rsin y2)其中R 6370为地球半径。A, B两点的实际距离ICOA OBd Rarccos；r ,OA OB化简得dRarccoscos(x1 x2) cosy1cos y2 siny1siny2。求解的模拟退火算法描述如下：(1)解空间解空间S可表为1,2,101,102的所有固定起点和终点的循环排列集合，即S ( 1, 102)| 11,( 2, 101)为2,3,101的循环排列，102102其中每一个循环排列表示侦察 100个目标的一个回路， i j表示在第i次侦察j点, 初始解可选为(1,2, ,102),本文中我们使用 MonteCarlo方法求得一个较好的初始解

17、。(2)目标函数此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求101min f( 1, 2, 102) d i i11 1而一次迭代由下列三步构成：(3)新解的产生2变换法任选序号 u,v ( u v)交换u与v之间的顺序，此时的新路径为：v11023变换法任选序号 u,v和w ,将u和v之间的路径插到w之后，对应的新路径为(设u v w)1 u 1 v 1 w u v w 1102(4)代价函数差d)v v 1 '对于2变换法，路径差可表示为 f (d d ) (d' u 1 vu v 1 7' u 1 I(5)接受准则1f 0Pexp( f /T) f

18、 0如果f0,则接受新的路径。否则，以概率exp(f/T)接受新的路径，即若exp( f /T)大于0到1之间的随机数则接受。(6)降温利用选定的降温系数进行降温即：T T ,得到新的温度，这里我们取0.9999。(7)结束条件用选定的终止温度 e 10 3°,判断退火过程是否结束。若 T e,算法结束，输出 当前状态。我们编写如下的 matlab程序如下：clc,clearload sj.txt%加载敌方100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中x=sj(:,1:2:8);x=x(:);y=sj(:,2:2:8);y=y(:);sj=x y;d1=70,

19、40;sj=d1;sj;d1;sj=sj*pi/180;%距离矩阵dd=zeros(102);for i=1:101for j=i+1:102temp=cos(sj(i,1)-sj(j,1)*cos(sj(i,2)*cos(sj(j,2)+sin(sj(i,2)*sin(sj(j,2);d(i,j)=6370*acos(temp); endendd=d+d'S0=;Sum=inf;rand('state',sum(clock);for j=1:1000S=1 1+randperm(100),102;temp=0;for i=1:101temp=temp+d(S(i),S

20、(i+1);endif temp<SumS0=S;Sum=temp;endende=0.1A30;L=20000;at=0.999;T=1;%退火过程for k=1:L%产生新解c=2+floor(100*rand(1,2);c=sort(c);c1=c(1);c2=c(2);%计算代价函数值df=d(S0(c1-1),S0(c2)+d(S0(c1),S0(c2+1)-d(S0(c1-1),S0(c1)-d(S0(c2),S0(c2+1);%接受准则if df<0S0=S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102);Sum=Sum+df;elseif ex

21、p(-df/T)>rand(1)S0=S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102);Sum=Sum+df;endT=T*at;if T<ebreak;endend%输出巡航路径及路径长度S0,Sum计算结果为44小时左右。其中的一个巡航路径如下图所示：§ 3蚁群算法3.1 蚁群算法简介蚁群是自然界中常见的一种生物，人们对蚂蚁的关注大都是因为“蚁群搬家，天 要下雨”之类的民谚。然而随着近代仿生学的发展，这种似乎微不足道的小东西越来越 多地受到学者们地关注。1991年意大利学者 M. Dorigo等人首先提出了蚁群算法，人们开始了对蚁群的研究： 相

22、对弱小，功能并不强大的个体是如何完成复杂的工作的（如寻找到食物的最佳路径并返回等)。在此基础上一种很好的优化算法逐步发展起来。蚁群算法的特点是模拟自然界中蚂蚁的群体行为。科学家发现，蚁群总是能够发现从蚁巢到食物源的最短路径。经研究发现，蚂蚁在行走过的路上留下一种挥发性的激素，蚂蚁就是通过这种激素进行信息交流。蚂蚁趋向于走激素积累较多的路径。找到最短路径的蚂蚁总是最早返回巢穴，从而在路上留下了较多的激素。由于最短路径上积累了较多的激素，选择这条路径的蚂蚁就会越来越多，到最后所有的蚂蚁都会趋向于选择这条最短路径。基于蚂蚁这种行为而提出的蚁群算法具有群体合作，正反馈选择，并行计算等三大特点，并且可以

23、根据需要为人工蚁加入前瞻、回溯等自然蚁所没有的特点。在使用蚁群算法求解现实问题时，先生成具有一定数量蚂蚁的蚁群，让每一只蚂蚁建立一个解或解的一部分，每只人工蚁从问题的初始状态出发，根据“激素”浓度来选择下一个要转移到的状态，直到建立起一个解，每只蚂蚁根据所找到的解的好坏程度在所经过的状态上释放与解的质量成正比例的“激素”。之后，每只蚂蚁又开始新的求解过程，直到寻找到满意解。为避免停滞现象，引入了激素更新机制。3.2 解决TSP问题的蚁群算法描述现以TSP问题的求解为例说明蚁群系统模型。首先引进如下记号：n为城市的个数；m为蚁群中蚂蚁的数量；dj为两城市i和j之间距离；b (t)为t时刻位于城市

24、i的n蚂蚁的个数，mbi(t)； j(t)为t时刻边弧(i, j)的轨迹强度(即ij连线上残留的i 1信息量)，且设j (0) c (c为常数)，i,j 1,2, ,n,i j； j(t)为t时刻边弧(i, j) 的能见度，反映由城市i转移到城市j的期望程度。根据上述原理，蚂蚁 k(k 1,2,m)在运动过程中根据各条路径上的信息量决定转移方向。与真实蚁群系统不同，人工蚁群系统具有一定的记忆功能。随着时间的推移，以前留下的信息逐渐消逝，经n个时刻，蚂蚁完成一次循环，各路径上信息量要作调整。由此得到下述的人工蚁群系统模型：1)设人工蚁群在并行地搜索 TSP的解，并通过一种信息素做媒介相互通信，在

25、每 个结点上且和该结点相连的边上以信息素量做搜索下一结点的试探依据，直到找到一个TSP问题的可行解。2)在时刻t人工蚁k由位置 j (t) j (t)Pik(t)iv(t) iv(t),v S0,其中参数为轨迹的相对重要性i转移至位置j的转移概率为j Sj S(0);为能见度的相对重要性(3)0); S为可行点集，即蚂蚁 k下一步允许选择的城市。分别反映了蚂蚁在运动过程中所积累的信息及启发式因子在蚂蚁选择路径中所起的不同作用。3)当m个人工蚁按(3)式找到了可行解，则将各边的信息量用下式修改。即调 整信息量的轨迹强度更新方程为(4)j(t 1) j(t)/(0,1)m kijijk 1其中 j

26、为第k只蚂蚁在本次循环中留在路径(i, j)上的信息量；j为本次循环中路径（i, j）上的信息量的增量；参数为轨迹的持久性；1为轨迹衰减度，表示信息消逝程度。对上述系统模型，采用人工蚁群方法求解的算法步骤可归结为：step 1： NC 0 （ NC为迭代步数或搜索次数）；各j和 j的初始化；将m个 蚂蚁置于n个顶点上。step 2：将各蚂蚁的初始出发点置于当前解集中；对每个蚂蚁k （k 1,2, ,m）按概率pj转移至下一顶点j ;将顶点j置于当前解集。step 3：计算各蚂蚁的目标函数值Zk（k 1, ,m）,记录当前的最好解。step 4：按更新方程修改轨迹强度。step 5： NC NC

27、 1,若NC 预定的迭代次数且无退化行为（即找到的都是 相同解），则转step 2。若为了简化计算，增加处理较大规模的TSP问题的能力，则可将（4）式修改为：ij(t 1) ij(t) (1) j ,(0,1)其中dij(i,j)BEQ 其它此处BE为本次最优路线上的边集。3.3人工蚁群算法性能的讨论人工蚁群算法是一种基于种群的进化算法。作为一个新兴的研究领域，虽它还远 未像GA、SA等算法那样形成系统的分析方法和坚实的数学基础，但目前已有一些基 本结果。在M. Dorigo三种不同的模型中，循环路径（i, j）上信息量的增量 j不同。1） Ant-quantity system 模型中，L,

28、若第k只蚂蚁在时刻t和t 1之间经过ijkdijdij0, 其它2）在 Ant-density system 模型中，k Q,若第k只蚂蚁在时刻t和t 1之间经过ijij 0, 其它3)在 Ant-cycle system 模型中，k ,若第k只蚂蚁在本次循环中经 过ijikLk0, 其它其中Q是反映蚂蚁所留轨迹数量的常数，Lk表示第k只蚂蚁在本次循环中所走路径的kk长度；且t 0时，ij （0） c, ij 0。算法中模型1）、2）利用的是局部信息，模 型3）利用的是整体信息。人工蚁群算法中，、等参数对算法性能也有很大的影响。值的大小表明留在每个结点上的信息量受重视的程度，值越大，蚂蚁选择以

29、前选过的点的可能性越大，但过大会使搜索过早陷于局部极小点；的大小表明启发式信息受重视的程度；Q值会影响算法的收敛速度，Q过大会使算法收敛于局部极小值，过小又会影响算法的收敛速度，随问题规模的增大Q的值也需要随之变化；蚂蚁的数目越多，算法的全局搜索能力越强，但数目加大将使算法的收敛速度减慢。第二十四章时间序列模型时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序 列的方法构成数据分析的一个重要领域，即时间序列分析。时间序列根据所研究的依据不同，可有不同的分类。1 .按所研究的对象的多少分，有一元时间序列和多元时间序列。2 .按时间的连续性可将时间序列分为离散时间序列和连续时间

30、序列两种。3 .按序列的统计特性分，有平稳时间序列和非平稳时间序列。如果一个时间序列 的概率分布与时间t无关，则称该序列为严格的(狭义的)平稳时间序列。如果序列的 一、二阶矩存在，而且对任意时刻t满足：(1)均值为常数(2)协方差为时间间隔的函数。则称该序列为宽平稳时间序列，也叫广义平稳时间序列。我们以后所研究的时间序列主要是宽平稳时间序列。4 .按时间序列的分布规律来分，有高斯型时间序列和非高斯型时间序列。§ 1确定性时间序列分析方法概述时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势 的。一个时间序列往往是以下几类变化形式的叠加或耦合。(1)长期趋势变动。它

31、是指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，它反映了客观事物的主要变化趋势。(2)季节变动。(3)循环变动。通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相 似的波动。(4)不规则变动。通常它分为突然变动和随机变动。通常用Tt表示长期趋势项，St表示季节变动趋势项， Ct表示循环变动趋势项，Rt表示随机干扰项。常见的确定性时间序列模型有以下几种类型：(1)加法模型yt Tt St Ct Rt(2)乘法模型yt Tt St Ct Rt(3)混合模型yt Tt St Rtyt St Tt Ct Rt其中yt是观测目标的观测记录，E(Rt) 0, E(Rt2)2。如果

32、在预测时间范围以内，无突然变动且随机变动的方差2较小，并且有理由认为过去和现在的演变趋势将继续发展到未来时，可用一些经验方法进行预测， 具体方法如下：1.1移动平均法设观测序列为y1， ，丫1，取移动平均的项数 N T 。一次移动平均值计算公式为：一 I /、M t(yt yt iyt n i)N1/ (yt 1yt N )(yt yt N ) M t 1(yt yt N )NNN二次移动平均值计算公式为：(2) _1(i)(i)(2) _1(i)(1) M t 77 (M tM t n i) M t i TT (M t M t n )NN当预测目标的基本趋势是在某一水平上下波动时，可用一次移

33、动平均方法建立预 测模型：(1) 1 .?ti M t 77(其？tNi),t N, N 1,N其预测标准误差为：T c2(夕 yt)2s ' tN -，T N最近N期序列值的平均值作为未来各期的预测结果。一般N取值范围：5 N 200。当历史序列的基本趋势变化不大且序列中随机变动成分较多时，N的取值应较大一些。否则 N的取值应小一些。在有确定的季节变动周期的资料中，移动 平均的项数应取周期长度。选择最佳N值的一个有效方法是，比较若干模型的预测误差。均方预测误差最小者为好。当预测目标的基本趋势与某一线性模型相吻合时，常用二次移动平均法，但序列 同时存在线性趋势与周期波动时，可用趋势移动

34、平均法建立预测模型：?T m aT bTm, m 1,2,其中 aT 2M： MT2), bT 2(M M，)。N 1例1某企业1月11月份的销售收入时间序列如下表所示。取 N 4,试用简 单一次滑动平均法预测第12月份的销售收入，并计算预测的标准误差。月份t销售收入yt1533.82574.63606.94649.85705.16772.0月份t789101112销售收入yt816.4892.7963.91015.11102.7解：首先计算出4,5,11M *993.6 ,预测的标准误差为1 ,、，Mt-(yt yt iyt 3), t4由于？ i Mt，t 4,5,11 ,则夕i211(q

35、yt)2150.5t 511 4计算的Matlab程序如下:y=533.8574.6606.9649.8705.1 772.0.816.4892.7963.91015.11102.7;temp=cumsum(y);mt=(temp(4:11)-0 temp(1:7)/4y12=mt(end)ythat=mt(1:end-1);fangcha=mean(y(5:11)-ythat).A2);sigma=sqrt(fangcha)1.2指数平滑法一次移动平均实际上认为最近N期数据对未来值影响相同，都加权 ；而N期N以前的数据对未来值没有影响，加权为0。但是，二次及更高次移动平均数的权数却不-1是且

36、次数越高，权数的结构越复杂，但永远保持对称的权数，即两端项权数小，N中间项权数大，不符合一般系统的动态性。 一般说来历史数据对未来值的影响是随时间 间隔的增长而递减的。所以，更切合实际的方法应是对各期观测值依时间顺序进行加权 平均作为预测值。指数平滑法可满足这一要求，而且具有简单的递推形式。设观测序列为y1，,yT ,为加权系数，01, 一次指数平滑公式为：Styt (1)St(11St1(ytSt1)(1)假定历史序列无限长，则有St(1)yt(1) yt1 (1居2(1)jyt j j 0(2)式表明St 是全部历史数据 的加权平均，加权系数分别为 2一一，,(1), (1),;显然有(1

37、j 0)j11 (1)由于加权系数序列呈指数函数衰减，加权平均又能消除或减弱随机干扰的影响，所以(2)称为一次指数平滑，类似地，二次指数平滑公式为：StSt(1) (1)St(21)(3)同理，三次指数平滑公式为：St(3)St(1)St(31)(4)一般P次指数平滑公式为：St(P)St(P1) (1)St(P1)(5)不管序列的基本趋势多但计算量很大。因此利用指数平滑公式可以建立指数平滑预测模型。原则上说,么复杂，总可以利用高次指数平滑公式建立一个逼近很好的模型， 用的较多的是几个低阶指数平滑预测模型。(1) 一次指数平滑预测?t 1 St，(2)线性趋势预测模型- Brown单系数线性平

38、滑预测(二次指数平滑预测)?t m at btm, m 1,2,其中 at2St(1) St，bt(St(1)St(2)1(3)二次曲线趋势预测模型Brown单系数二次式平滑预测12?t m at btm -Ctm , m 1,2, 2其中 at3st3st St(3),bt 2(6 5 )St 2(5 4 )S，(4 3 闾3),2(1)指数平滑预测模型是以时刻t为起点，综合历史序列的信息，对未来进行预测的。选择合适的加权系数是提高预测精度的关键环节。根据实践经验，的取值范围一般以0.10.3为宜。值愈大，加权系数序列衰减速度愈快，所以实际上取值大小起着控制参加平均的历史数据的个数的作用。值

39、愈大意味着采用的数据愈少。因此，可以得到选择值的一些基本准则。(1)如果序列的基本趋势比较稳，预测偏差由随机因素造成，则 值应取小一些，以减少修正幅度，使预测模型能包含更多历史数据的信息。(2)如果预测目标的基本趋势已发生系统地变化，则值应取得大一些。这样，可以偏重新数据的信息对原模型进行大幅度修正，以使预测模型适应预测目标的新变 化。另外，由于指数平滑公式是递推计算公式，所以必须确定初始值S01), s02) , s03)。可以取前35个数据的算术平均值作为初始值。例2下表数据是某股票在8个连续交易日的收盘价，试用一次指数平滑法预测第9个交易日的收盘价(初始值S01)Yi,0.4)。时间 t

40、 12345678价格 Yt16.41 17.62 16.15 15.54 17.24 16.83 18.14 17.05解： 由于S01)Yi, StYt (1)St(11, t 1,2,8求得?9S81)17.18预测标准误差为：8(?tYt)2S270.96计算的Matlab程序如下：alpha=0.4;Y=16.4117.62 16.15 15.54 17.24 16.83 18.14 17.05;s1(1)=Y(1)； for i=2:8 s1(i)=alpha*Y(i)+(1-alpha)*s1(i-1);end Yhat9=s1(end)sigma=sqrt(mean(s1(1:

41、end-1)-y(2:end).A2)例3 仍以例 2为例。试用两次指数平滑法预测第9个交易日的收盘价(S01)S02)yi,0.4)。解：a82s8"S8(2)17.38,b8(S：)S,)0.13 ,于是1?9a8b8 17.51。下面进行预测误差分析，取?t1atbt(2StSt(2)-(StSt(2)St丁(StSt)11预测标准误差18(? yt)2 1|2J-T2- 1.21计算的Matlab程序如下:clc,clearalpha=0.4;y=16.4117.62 16.15 15.54 17.24 16.83 18.14 17.05;s1(1)=y(1);for i=2

42、:8s1(i)=alpha*y(i)+(1-alpha)*s1(i-1);ends2=y(1);for i=2:8s2(i)=alpha*s1(i)+(1-alpha)*s2(i-1);enda8=2*s1(8)-s2(8)b8=alpha/(1-alpha)*(s1(8)-s2(8)yhat9=a8+b8yhat(1)=y(1)for i=2:8yhat(i)=s1(i-1)+1/(1-alpha)*(s1(i-1)-s2(i-1);endtemp=sum(yhat-y).A2);sigma=sqrt(temp/6)§ 2平稳时间序列模型这里的平稳是指宽平稳，其特性是序列的统计特性

43、不随时间的平移而变化，即均 值和协方差不随时间的平移而变化。下面自回归模型(Auto Regressive Model )简称AR模型，移动平均模型(Moving Average Model)简称 MA 模型，自回归移动平均模型( Auto Regressive Moving Average Model)简称ARMA模型。下面的 Xt为零均值(即中心化处理的)平稳序列。(1) 一般自回归模型 AR(n)假设时间序列 Xt仅与Xt1,Xt2, ,Xtn有线性关系，而在 Xt 1,Xt 2, ,Xt n 已知条件下，Xt与Xt j (j n 1,n 2,)无关，at是一个独立于 Xt 1,Xt 2

44、, ,Xtn 的白噪声序列，at - N(0, 2)。X t 1Xt 12Xt 2n X t n at上式还可以表不为at Xt1Xt 12Xt 2nXt n可见，AR(n)系统的响应 Xt具有n阶动态性。 AR(n)模型通过把 Xt中的依赖于 Xt 1, Xt 2, , Xt n的部分消除掉之后，使得具有 n阶动态性的序列 Xt转化为独立的 序列at。因此，拟合 AR(n)模型的过程也就是使相关序列独立化的过程。(2)移动平均模型 MA (m)AR (n)系统的特征是系统在t时亥的响应Xt仅与其以前时刻的响应Xt 1,Xt 2, ,Xt n有关，而与其以前时刻进入系统的扰动无关。如果一个系统

45、在t时刻的响应Xt,与其以前时刻t 1,t t 1,t 2, ,t m进入系统的扰动 这一类系统为 MA (m)系统。X t at 1 at 12 at 22, 的响应X，Xt2,无关，而与其以前时刻at 1,at 2, ,at m存在着一定的相关关系，那么,mat m(3)自回归移动平均模型一个系统，如果它在时刻 t的响应Xt ,不仅与其以前时刻的自身值有关，而且还那么，这个系统就是自回归移动平与其以前时刻进入系统的扰动存在一定的依存关系, 均系统。ARMA (n, m)模型为X t Xt 1nXt n at 向 1对于平稳系统来说，由于AR、MA、ARMAmat m(n,m)模型都是ARM

46、A(n,n 1)模型的特例，我们以 ARMA (n,n 1)模型为一般形式来建立时序模型。§ 3 ARMA模型的特性在时间序列的时域分析中，线性差分方程是极为有效的工具。事实上，任何一个 ARMA模型都是一个线性差分方程。3.1 AR(1)系统的格林函数格林函数就是描述系统记忆扰动程度的函数。AR(1)模型为显然Xtt1y(t日1Xt 1 aty(t),则有1) 1y(t)at个一阶非齐次差分方程。(6)由于Xt 1Xt2X1 X t 23X1八t 3atat 1at 21 ( 1 Xt 2at 1 )atatXt依次递推下去，可得到1 at j显然(7)式是差分方程(6)的解。方程

47、解的系数函数/客观地描述了该系统的动态性，故这个系统函数就叫做记忆函数，也叫格林函数。若用 Gj表示，则AR(1)模型的格林函数可以表示为Gjij这样(7)式可等价地写成：XtGjat jj 0上式也可写成tXt Gt kak k这里G00 1。定义后移算子B , BXt(1 iB)Xt at 它的解为1Xt-at (11B11Bat1at 11 at 2Xt 1，B2Xt Xt 2，这样，AR(1)可写成12B2)at(8)G jat j j 0由于格林函数就是差分方程解的系数函数，格林函数的意义可概括如下：(1) Gj是前j个时间单位以前进入系统的扰动at j对系统现在行为(响应)影响的权

48、数。(2) Gj客观地刻画了系统动态响应衰减的快慢程度。(3)对于一个平稳系统来说，在某一时刻由于受到进入系统的扰动at的作用，离开其平衡位置(即平均数-零)，Gj描述系统回到平衡位置的速度，1的值较小，速度较快；1的值较大，回复的速度就较慢。3.2 ARMA (2,1)系统的格林函数(1) ARMA (2,1)系统的格林函数的隐式ARMA (2,1)模型是一个二阶非齐次差分方程Xt 1 Xt 12 Xt 2 at1at 1设该二阶非齐次差分方程的解为(9)(10)XtGjat jj 0 为方便起见，可用 B算子：(11B 2B2)Xt (11B)atXt GjBj at(11)j 0把(11

49、)式代入(10),比较两边B的同次哥的系数得Go 1, Gi 11, Gj £j 12Gj 2,j 3,4,(2) ARMA (2,1)系统的格林函数的显式ARMA (2,1)模型实质上是一个二阶非齐次差分方程：XtXt 12Xt 2 at 向 1欲求其解，必须先求出其相应的齐次差分方程的通解。 齐次方程对应的特征方程为20,124 221特征根为124 22齐次差分方程的通解为G jC1 1C2 2其中C1和C2是任意常数，其值由初始条件唯一地确定。这里的初始条件为:Go 1于是有GoGC21G1G 1C2 211而121 ,即G c2 1G 1c2 2121解之，得1121C1,

50、 C21221则ARMA (2,1)系统的格林函数为:11 j 21 jG j1212213.3逆函数和可逆性前面的格林函数，把 Xt表示为过去at对Xt的影响，或者说系统对过去at的记忆性，也就是用一个 MA模型来逼近 Xt的行为。平稳序列 Xt的这种表达形式称为 Xt的 “传递形式”。同样我们也可以用过去的Xt的一个线性组合来逼近系统现在时刻的行为。即XtIiXt 1I 2X t 2atIjXt j j 1at我们把这种表达形式称为Xt的“逆转形式”。其中的系数函数Ij(l01)称为逆函数。可见它是一个无穷阶的自回归模型。一个过程是否具有逆转形式，也就是说逆函数是否存在的性质，通常称为过程是否具有可逆性，如果一个过程可以用一个无限阶的自回归模型逼近，即逆函数存在，我们就称该过程具有可逆性，否则，就是不可逆的。对于AR (2)模型Xt1 Xt 12Xt 2 at有l11, I