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

jj)E(i)，接受该状态被转换。j)E(i)，则状态转换以如下概率被接受：

E(i)Ejj)E(i)，接受该状态被转换。j)E(i)，则状态转换以如下概率被接受：

E(i)E(j)KT

eE(i)KT

exE(j)KT

§1现代优化算法简介现代优化算法是80年代初兴起的启发式算法。这些算法包括禁忌搜索（tabusearch），模拟退火（simulatedannealing），遗传算法（geneticalgorithms），人工神经网络（neuralnetworks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标－求NP-hard组合优化问题的全局最优解。虽然有这些目标，但NP-hard理论限制它们只能以启发式的算法去求解实际问题。启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（AntColonyAlgorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。现代优化算法解决组合优化问题，如TSP（TravelingSalesmanProblem）问题，QAP（QuadraticAssignmentProblem）问题，JSP（Job-shopSchedulingProblem）问题等效果很好。本章我们只介绍模拟退火算法，初步介绍一下蚁群算法，其它优化算法可以参看相关的参考资料。

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

e其中K是物理学中的波尔兹曼常数，T是材料温度。在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态i的概率满足波尔兹曼分布：

P(xi)

jS其中表示材料当前状态的随机变量，S表示状态空间集合。显然

-271-

e

S|表示集合S

eeE(i)EeeejSeE(i)E

minE(j)且S{i|E(i)E}。minmin

:xR，其中xS{y|e

S|表示集合S

eeE(i)EeeejSeE(i)E

minE(j)且S{i|E(i)E}。minmin

:xR，其中xS{y|yR,y0}，SS

x(0)生成下一个x'N(x(0))x'作为一个新解x(1)依赖于下面概率：

x'的函数值比前一个解的函数值更小，则接受x(1)x'作为T0

x(k1)

1

eTii1ii1

x(k1)

ejS

E(E(E(KT

emin

x

0

x'作为一个新解。

i

若f(x')f(x(k))

E(i)KT

E(E(i)EminKTKT

j)Eminj)Eminj)EminKTKTKTminS|E(ef(x')f(x(k))T0j)KTmin

|j)EminKT（1）其它1|S|1minf(x')f(x(0))T0min若iSmin0其它jSmin

其它T

其中|中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，

limlimT0T0

jS

limT0

jS其中EminjS上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。考虑这样一个组合优化问题：优化函数为F，它表示优化问题的一个可行解，R表示函数的定义域。N(x)表示的一个邻域集合。首先给定一个初始温度T和该优化问题的一个初始解x(0)，并由解，是否接受1若f(x')f(x(0))P(x(0)x')

换句话说，如果生成的解

f(x')f(x(0))一个新解。否则以概率e接受泛泛地说，对于某一个温度T和该优化问题的一个解x(k)，可以生成x'。接受x'作为下一个新解的概率为：

P(x(k)x')

在温度T下，经过很多次的转移之后，降低温度T，得到T。在T过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。我们注意到，在每个T下，所得到的一个新状态

-272-

,x(k1)无关，因此这是一个马尔可夫过程。使用马x'的N(x(k))中是均匀分布的，且新状态x'被接受的概率满足式（1），那么经过xieiex

P*1min已知敌方100个目标的经度、纬度如下：经度51.175810.81981.945144.0356,x(k1)无关，因此这是一个马尔可夫过程。使用马x'的N(x(k))中是均匀分布的，且新状态x'被接受的概率满足式（1），那么经过xieiex

P*1min已知敌方100个目标的经度、纬度如下：经度51.175810.81981.945144.035632.191016.561822.10750.773253.352451.961215.254834.1688f(xi)Ti若x|纬度0.032216.25290.205713.54015.869923.614318.55690.465626.725622.851127.211122.7571f(xi)Ti经度46.325322.789126.495128.983636.486310.55970.121547.413430.816512.79386.20709.4402（2）minmin0纬度28.275323.104522.122125.987929.728415.117818.872623.778313.459515.73075.14423.9200其它经度30.331310.158431.484738.47220.971850.211148.207741.867127.71334.956849.243011.5812纬度6.934812.48198.964020.173128.147710.294416.88893.56675.07068.366916.704414.5677尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态x(k)生成概率，在有限次的转换，在温度T下的平衡态

P(T)i

jS当温度T降为0时，1PS|i

并且ixiS这说明如果温度下降十分缓慢而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找到全局最优解。在模拟退火算法中应注意以下问题：（1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。（2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续m次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值T（3）选择初始温度和确定某个可行解的邻域的方法也要恰当。2.2应用举例例经度纬度53.712115.304656.543221.418820.105015.456226.241818.176028.269429.00118.958624.66358.15199.532531.949917.630943.54743.906123.92227.630621.505124.090917.116820.0354

-273-

0.408816.847413.636827.14858.46481.951014.297829.511427.670516.835919.569723.18043.1853(d)djj1,2,,102，这里Dij120.408816.847413.636827.14858.46481.951014.297829.511427.670516.835919.569723.18043.1853(d)djj1,2,,102，这里Dij122XOZ6370

21212101,102}的所有固定起点和终点的循环排列集合，即j(1,2,,102)9.555935.661919.866039.516851.818134.057458.828947.50999.155649.981611.511838.339240.1400102102表示表示i,两点的距离，i,为实对称矩阵。则问题是

，(xXOY2222

OAOBOAOB)为{2,3,,101}的循环排列,102ji，本文中我们使用MonteCarlo方法求得一个较好的初始解。11.42199.933315.122416.937123.015923.396014.522924.066414.13046.082817.388419.995020.3030)，过A,B两点的大圆的劣弧长即为两点2，102}24.450924.46543.161656.50898.998323.062418.663510.112153.798919.363544.039824.654323.98766.56343.16444.242813.709023.64408.43196.743627.26620.219917.662216.263519.60579.403026.72130.777518.524552.521150.115619.985752.842328.781233.649036.954539.713936.998041.108428.56676.957614.359815.795723.78165.790227.288027.66590.398023.026528.420324.399227.714937.584814.470358.684938.430013.790940.880139.94948.08311.349615.73206.99094.1591

我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵D其中求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A,B两点的地理坐标分别为(x,y),的实际距离。以地心为坐标原点O，以赤道平面为平面，以0度经线圈所在的平面为平面建立三维直角坐标系。则A,B两点的直角坐标分别为：A(Rcosxcosy,Rsinxcosy,Rsiny)B(Rcosx1cosy1,Rsinx1cosy1,Rsiny1)其中R为地球半径。A,B两点的实际距离dRarccos

化简得dRarccos[(xcosx)cosycosysinysiny求解的模拟退火算法描述如下：（1）解空间解空间S可表为{1,2,,

S{(,,)|1,(,,110212101其中每一个循环排列表示侦察100个目标的一个回路，表示在第i次侦察点，初始解可选为（2）目标函数-274-

102

102ii1

u,v1u1vv1u1uv1102u,v

u1v1wuvw1102

u1vv11f0/T)f00，则接受新的路径。否则，以概率exp(f/T)0到1之间的随机数则接受。

0.9999

1030e，算法结束，输出101

i1

v）交换u与v之间的顺序，此时的新路径为：

w，将uvw

uv1u1uv

TT。minf(,,,)2而一次迭代由下列三步构成：（3）新解的产生①2变换法任选序号（u②3变换法任选序号和和之间的路径插到之后，对应的新路径为（设uvw）（4）代价函数差对于2变换法，路径差可表示为f(dd)(dd（5）接受准则Pexp(f如果f接受新的路径，即若exp(f/T)大于（6）降温利用选定的降温系数进行降温即：，得到新的温度，这里我们取（7）结束条件用选定的终止温度e，判断退火过程是否结束。若T当前状态。我们编写如下的matlab程序如下：clc,clearloadsj.txt%加载敌方100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中x=sj(:,1:2:8);x=x(:);y=sj(:,2:2:8);y=y(:);sj=[xy];d1=[70,40];sj=[d1;sj;d1];sj=sj*pi/180;%距离矩阵dd=zeros(102);fori=1:101forj=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);end

-275-

endd=d+d';S0=[];Sum=inf;rand('state',sum(clock));forj=1:1000S=[11+randperm(100),102];temp=0;fori=1:101temp=temp+d(S(i),S(i+1));endiftemp<SumS0=S;Sum=temp;endende=0.1^30;L=20000;at=0.999;T=1;%退火过程fork=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));%接受准则ifdf<0S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;elseifexp(-df/T)>rand(1)S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;endT=T*at;ifT<ebreak;endend%输出巡航路径及路径长度S0,Sum

计算结果为44小时左右。其中的一个巡航路径如下图所示：

-276-

0102010203040506070

j

bj)的轨迹强度（即ij

(0i)1c（c为常数）i,j1,2,,n，ij；(t)为t时刻边弧(i,j)ijij

1,2,,m)在运动过程中根据各条路径上的信息量决定

n个时刻，蚂蚁完成一次循环，各路径上信息量要作调整。i(t)为tij(t)

jn

35

30

25

20

15

10

5

0

§3蚁群算法3.1蚁群算法简介蚁群是自然界中常见的一种生物，人们对蚂蚁的关注大都是因为“蚁群搬家，天要下雨”之类的民谚。然而随着近代仿生学的发展，这种似乎微不足道的小东西越来越多地受到学者们地关注。1991年意大利学者M.Dorigo等人首先提出了蚁群算法，人们开始了对蚁群的研究：相对弱小，功能并不强大的个体是如何完成复杂的工作的（如寻找到食物的最佳路径并返回等）。在此基础上一种很好的优化算法逐步发展起来。蚁群算法的特点是模拟自然界中蚂蚁的群体行为。科学家发现，蚁群总是能够发现从蚁巢到食物源的最短路径。经研究发现，蚂蚁在行走过的路上留下一种挥发性的激素，蚂蚁就是通过这种激素进行信息交流。蚂蚁趋向于走激素积累较多的路径。找到最短路径的蚂蚁总是最早返回巢穴，从而在路上留下了较多的激素。由于最短路径上积累了较多的激素，选择这条路径的蚂蚁就会越来越多，到最后所有的蚂蚁都会趋向于选择这条最短路径。基于蚂蚁这种行为而提出的蚁群算法具有群体合作，正反馈选择，并行计算等三大特点，并且可以根据需要为人工蚁加入前瞻、回溯等自然蚁所没有的特点。在使用蚁群算法求解现实问题时，先生成具有一定数量蚂蚁的蚁群，让每一只蚂蚁建立一个解或解的一部分，每只人工蚁从问题的初始状态出发，根据“激素”浓度来选择下一个要转移到的状态，直到建立起一个解，每只蚂蚁根据所找到的解的好坏程度在所经过的状态上释放与解的质量成正比例的“激素”。之后，每只蚂蚁又开始新的求解过程，直到寻找到满意解。为避免停滞现象，引入了激素更新机制。3.2解决TSP问题的蚁群算法描述现以TSP问题的求解为例说明蚁群系统模型。首先引进如下记号：n为城市的个数；m为蚁群中蚂蚁的数量；d为两城市i和时刻位于城市i的

蚂蚁的个数，m为t时刻边弧(i,连线上残留的

信息量），且设的能见度，反映由城市i转移到城市的期望程度。根据上述原理，蚂蚁k(k转移方向。与真实蚁群系统不同，人工蚁群系统具有一定的记忆功能。随着时间的推移，以前留下的信息逐渐消逝，经

-277-

j

00j

00

m个人工蚁按（3）式找到了可行解，则将各边的信息量用下式修改。即调

(t1)(t)(0,1)

kj)上的信息量；j)上的信息量的增量；参数1

0（NCij

1,2,,m）pk转移至下一顶点jjij1,,m)，记录当前的最好解。

NCNC1，若NC

(t1)(t)(1

k

BEijij（3）

，（4）

m

k1kij

，

d(t)(t)

,(i,j)BE,jS

iviv0,jS

vS

ij0,其它1）设人工蚁群在并行地搜索TSP的解，并通过一种信息素做媒介相互通信，在每个结点上且和该结点相连的边上以信息素量做搜索下一结点的试探依据，直到找到一个TSP问题的可行解。2）在时刻t人工蚁k由位置i转移至位置的转移概率为(t)(t)pk(t)ij

其中参数为轨迹的相对重要性（）；为能见度的相对重要性（）；S为可行点集，即蚂蚁k下一步允许选择的城市。,分别反映了蚂蚁在运动过程中所积累的信息及启发式因子在蚂蚁选择路径中所起的不同作用。3）当整信息量的轨迹强度更新方程为ijijijkij其中为第只蚂蚁在本次循环中留在路径(i,径(i,为轨迹的持久性；为轨迹衰减度，表示信息消逝程度。对上述系统模型，采用人工蚁群方法求解的算法步骤可归结为：step1：NC为迭代步数或搜索次数）；各和的初始化；将m蚂蚁置于n个顶点上。step2：将各蚂蚁的初始出发点置于当前解集中；对每个蚂蚁k（k按概率；将顶点置于当前解集。step3：计算各蚂蚁的目标函数值zk(kstep4：按更新方程修改轨迹强度。step5：预定的迭代次数且无退化行为（即找到的都是相同解），则转step2。若为了简化计算，增加处理较大规模的TSP问题的能力，则可将（4）式修改为：)ijij其中1ij

此处为本次最优路线上的边集。3.3人工蚁群算法性能的讨论人工蚁群算法是一种基于种群的进化算法。作为一个新兴的研究领域，虽它还远未像GA、SA等算法那样形成系统的分析方法和坚实的数学基础，但目前已有一些基-278-

j)上信息量的增量ij不同。

k

Q,若第只蚂蚁在时刻t和j)上信息量的增量ij不同。

k

Q,若第只蚂蚁在时刻t和t1之间经过ijk

k

t0时，(0)c，ij,,Q等参数对算法性能也有很大的影响。值的大小表明留

Q过大会使算法收敛于局部极小值，过小又会影响算法的Q的值也需要随之变化；蚂蚁的数目越多，算法的全局

t无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的t满足：d

0,

kk,若第k只蚂蚁在时刻t和t1之间经过ij

QL

其它

,若第k只蚂蚁在本次循环中过ijij0,

k0,其它

其它在M.Dorigo三种不同的模型中，循环路径(i,1）Ant-quantitysystem模型中，Qij

2）在Ant-densitysystem模型中，ij3）在Ant-cyclesystem模型中，

ij

其中Q是反映蚂蚁所留轨迹数量的常数，Lk表示第k只蚂蚁在本次循环中所走路径的长度；且型3）利用的是整体信息。人工蚁群算法中，在每个结点上的信息量受重视的程度，值越大，蚂蚁选择以前选过的点的可能性越大，但过大会使搜索过早陷于局部极小点；的大小表明启发式信息受重视的程度；Q值会影响算法的收敛速度，收敛速度，随问题规模的增大搜索能力越强，但数目加大将使算法的收敛速度减慢。第二十四章时间序列模型

时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序列的方法构成数据分析的一个重要领域，即时间序列分析。时间序列根据所研究的依据不同，可有不同的分类。1．按所研究的对象的多少分，有一元时间序列和多元时间序列。2．按时间的连续性可将时间序列分为离散时间序列和连续时间序列两种。3．按序列的统计特性分，有平稳时间序列和非平稳时间序列。如果一个时间序列的概率分布与时间一、二阶矩存在，而且对任意时刻（1）均值为常数（2）协方差为时间间隔的函数。则称该序列为宽平稳时间序列，也叫广义平稳时间序列。我们以后所研究的时间序列主要是宽平稳时间序列。4．按时间序列的分布规律来分，有高斯型时间序列和非高斯型时间序列。

§1确定性时间序列分析方法概述时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势-279-

ttt

tttt

tttt

tttt

tttt

yE(R)0，E(R2)t2

y,,T

N111NNN

(2)(1)1)(2)(1)1)

(1)(ˆy)2tttN1

NN的Nttttt2

T

tt1tN1

(1)t1tNttNt1ttN

11NN

1NTN，，（1）长期趋势变动。它是指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，它反映了客观事物的主要变化趋势。（2）季节变动。（3）循环变动。通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相似的波动。（4）不规则变动。通常它分为突然变动和随机变动。通常用T表示长期趋势项，S表示季节变动趋势项，C表示随机干扰项。常见的确定性时间序列模型有以下几种类型：（1）加法模型yTSC（2）乘法模型yTSC（3）混合模型yTSRytStTtCtR其中如果在预测时间范围以内，无突然变动且随机变动的方差较小，并且有理由认为过去和现在的演变趋势将继续发展到未来时，可用一些经验方法进行预测，具体方法如下：1.1移动平均法设观测序列为，取移动平均的项数N。一次移动平均值计算公式为：M(yyy

(yy)(yy)M(yy)

二次移动平均值计算公式为：M(MM)M(MM)tttN1t1ttN当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型：yˆM(yˆyˆ)N,N1,，t1tt其预测标准误差为：

TS

最近期序列值的平均值作为未来各期的预测结果。一般取值范围：5N200。当历史序列的基本趋势变化不大且序列中随机变动成分较多时，N取值应较大一些。否则的取值应小一些。在有确定的季节变动周期的资料中，移动

-280-

N

MM(MMN

MM(MM)。TTTT4，试用简

yt

yt

(1)

yˆMt4,5,,11yˆM993.6，预测的标准误差为t1t121111

ttt5

NN

1N

y,,1，一次指数平滑公式为：T(1)(2)(1)(2)

14

(1)，(1)

(yˆy)2

114

1N2N1

，t差。均方预测误差最小者为好。当预测目标的基本趋势与某一线性模型相吻合时，常用二次移动平均法，但序列同时存在线性趋势与周期波动时，可用趋势移动平均法建立预测模型：yˆabm，m1,2,TmTT其中a，b

例1某企业1月～11月份的销售收入时间序列如下表所示。取N单一次滑动平均法预测第12月份的销售收入，并计算预测的标准误差。

月份t123456销售收入533.8574.6606.9649.8705.1772.0月份t789101112销售收入816.4892.7963.91015.11102.7

解：首先计算出M(yyy)4,5,,11ttt1t3由于，则

S150.5

计算的Matlab程序如下：y=[533.8574.6606.9649.8705.1772.0...816.4892.7963.91015.11102.7];temp=cumsum(y);mt=(temp(4:11)-[0temp(1:7)])/4y12=mt(end)ythat=mt(1:end-1);fangcha=mean((y(5:11)-ythat).^2);sigma=sqrt(fangcha)

1.2指数平滑法

一次移动平均实际上认为最近期数据对未来值影响相同，都加权；而期

以前的数据对未来值没有影响，加权为0。但是，二次及更高次移动平均数的权数却不

是，且次数越高，权数的结构越复杂，但永远保持对称的权数，即两端项权数小，

中间项权数大，不符合一般系统的动态性。一般说来历史数据对未来值的影响是随时间间隔的增长而递减的。所以，更切合实际的方法应是对各期观测值依时间顺序进行加权平均作为预测值。指数平滑法可满足这一要求，而且具有简单的递推形式。设观测序列为，为加权系数，0

-281-

(1)(1)tt1

(1)(1)St,(1),(1)2,；显然有

S(2)tt1S(3)PS(P)tt1）线性趋势预测模型－Brown单系数线性平滑预测（二次指数平滑预测）abm，m1,2,

1abmtt3S(1)(1)(1)tt1

(1)(1)St,(1),(1)2,；显然有

S(2)tt1S(3)PS(P)tt1）线性趋势预测模型－Brown单系数线性平滑预测（二次指数平滑预测）abm，m1,2,

1abmtt3S(1)3S(2)ttt

2(1)的取值范围一值愈大，加权系数序列衰减速度愈快，所以实际上取值大小值愈大意味着采用的数据愈少。因此，值的一些基本准则。值应取小一些，

值应取得大一些。这样，(yS)t1tt1（2）(1)

(1)（3）tt1次指数平滑公式为：（5）t1C1,2,t(2)2（1）j（4）t1(43)St，m(3)t1假定历史序列无限长，则有Sy(1)[y(1)S](1)jtt1t2tjj0（2）式表明是全部历史数据的加权平均，加权系数分别为

由于加权系数序列呈指数函数衰减，加权平均又能消除或减弱随机干扰的影响，所以（2）称为一次指数平滑，类似地，二次指数平滑公式为：S(2)S(1)同理，三次指数平滑公式为：S(3)S(2)一般S(P)S(P1)利用指数平滑公式可以建立指数平滑预测模型。原则上说，不管序列的基本趋势多么复杂，总可以利用高次指数平滑公式建立一个逼近很好的模型，但计算量很大。因此用的较多的是几个低阶指数平滑预测模型。（1）一次指数平滑预测yˆS(1)，（2yˆtmtt2S(1)S(2)，b(S(1)ttttt（3）二次曲线趋势预测模型－Brown单系数二次式平滑预测yˆtm其中ab[(65)S(1)2(54)S

指数平滑预测模型是以时刻t为起点，综合历史序列的信息，对未来进行预测的。选择合适的加权系数是提高预测精度的关键环节。根据实践经验，般以0.1～0.3为宜。起着控制参加平均的历史数据的个数的作用。可以得到选择（1）如果序列的基本趋势比较稳，预测偏差由随机因素造成，则以减少修正幅度，使预测模型能包含更多历史数据的信息。（2）如果预测目标的基本趋势已发生系统地变化，则可以偏重新数据的信息对原模型进行大幅度修正，以使预测模型适应预测目标的新变-282-

S(1),S,S00

S(1)yS(1),S,S00

S(1)y0.401

yt

01tt1

988

ttt2

S(1)Sy0.4001a2S(1)S17.38，b(S(1)S(2))888

9881(ˆy)2ttt1(2)(3)

，S(1)(yˆy)2

7

(2)

182另外，由于指数平滑公式是递推计算公式，所以必须确定初始值可以取前3～5个数据的算术平均值作为初始值。例2下表数据是某股票在8个连续交易日的收盘价，试用一次指数平滑法预测第9个交易日的收盘价（初始值，。

时间t12345678价格16.4117.6216.1515.5417.2416.8318.1417.05

解：由于S(1)yyS(1)，t1,2,,8求得yˆS(1)17.18预测标准误差为：

S0.96

计算的Matlab程序如下：alpha=0.4;y=[16.4117.6216.1515.5417.2416.8318.1417.05];s1(1)=y(1);fori=2:8s1(i)=alpha*y(i)+(1-alpha)*s1(i-1);endyhat9=s1(end)sigma=sqrt(mean((s1(1:end-1)-y(2:end)).^2))

例3仍以例2为例。试用两次指数平滑法预测第9个交易日的收盘价（，。

解：，于是yˆab17.51。下面进行预测误差分析，取yˆab(2S(1)S(2))(S(1))1(S(1)S(2))t1ttttttttt预测标准误差

8S1.21

计算的Matlab程序如下：clc,clear

-283-

X

XXX

XX,X,,XX,X,,XXtX(jn1,n2,)att的白噪声序列，aXXXXtat1t12t2ntnt

AR(n)XnAR(n)Xtt的部分消除掉之后，使得具有转化为独立的at

系统的特征是系统在时刻的响应仅与其以前时刻的响应有关，而与其以前时刻进入系统的扰动无关。如果一个系统在t时X1,t2,的响应X,Xtt2

MA(m)t1t2tnt1t2tn

t

t2tm2

存在着一定的相关关系，那么，y=[16.4117.6216.1515.5417.2416.8318.1417.05];s1(1)=y(1);fori=2:8s1(i)=alpha*y(i)+(1-alpha)*s1(i-1);ends2=y(1);fori=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)fori=2:8yhat(i)=s1(i-1)+1/(1-alpha)*(s1(i-1)-s2(i-1));endtemp=sum((yhat-y).^2);sigma=sqrt(temp/6)

§2平稳时间序列模型这里的平稳是指宽平稳，其特性是序列的统计特性不随时间的平移而变化，即均值和协方差不随时间的平移而变化。下面自回归模型（AutoRegressiveModel）简称AR模型，移动平均模型（MovingAverageModel）简称MA模型，自回归移动平均模型（AutoRegressiveMovingAverageModel）简称ARMA模型。下面的（1）一般自回归模型AR(n)假设时间序列仅与有线性关系，而在已知条件下，与无关，是一个独立X,X,,~N(0,)。t1tna

上式还可以表示为aXXXXtt1t12t2ntn可见，系统的响应具有阶动态性。模型通过把中的依赖于X,X,,Xtnn阶动态性的序列Xt1t2t序列。因此，拟合AR(n)模型的过程也就是使相关序列独立化的过程。（2）移动平均模型MA(m)AR(n)tXX,X,,Xtnt1t2刻的响应，与其以前时刻tt1,t2,,tm进入系统的扰动a,a这一类系统为系统。

-284-

tt1t12t2mtm

X

ARMAt1(n,m)模型都是ARMA(n,tt1t12t2mtm

X

ARMAt1(n,m)模型都是ARMA(n,n1)模1)模型为一般形式来建立时序模型。

Xty(t)，则有

1t

Xaa11t31t21t1t

jGjj

t

1。01t1t1t

（6）

2t21t1t3

（8）

k0ja（7）（3）自回归移动平均模型一个系统，如果它在时刻t的响应与其以前时刻进入系统的扰动存在一定的依存关系，那么，这个系统就是自回归移动平均系统。ARMA(n,m)模型为XXXaaat1t1ntnt1mtm对于平稳系统来说，由于AR、MA、型的特例，我们以ARMA(n,n

§3ARMA模型的特性在时间序列的时域分析中，线性差分方程是极为有效的工具。事实上，任何一个ARMA模型都是一个线性差分方程。3.1AR(1)系统的格林函数格林函数就是描述系统记忆扰动程度的函数。AR(1)模型为XXa设y1(t1)y(t)a显然是一个一阶非齐次差分方程。由于XXa(Xa)at1t1t11t2t1t

X2aaa依次递推下去，可得到Xt1tjj0显然（7）式是差分方程（6）的解。方程解的系数函数函数，也叫格林函数。若用表示，则AR(1)模型的格林函数可以表示为Gj1这样（7）式可等价地写成：XGatjtjj0上式也可写成XGattkk这里G

-285-

BXX,B2XBXX,B2XX,，这样，AR(1)可写成tt1tt2

1tt11B2a

jtj

Gj

ARMA(2,1)系统的格林函数

t1t12t2t1t1

（10）

ARMA(2,1)系统的格林函数的显式

t1t12t2t1t1

0，1211t

a

t

（9）

111j1j12j2j

，G（11）(1B)Xa它的解为Xa(1B2B2

ata1t11t2Gtjj0由于格林函数就是差分方程解的系数函数，格林函数的意义可概括如下：（1）G是前个时间单位以前进入系统的扰动响的权数。（2）G（3）对于一个平稳系统来说，在某一时刻由于受到进入系统的扰动a的作用，离开其平衡位置（即平均数－零），描述系统回到平衡位置的速度，度较快；1的值较大，回复的速度就较慢。3.2（1）ARMA(2,1)系统的格林函数的隐式ARMA(2,1)模型是一个二阶非齐次差分方程XXXaa设该二阶非齐次差分方程的解为XGatjtjj0为方便起见，可用B算子：(1BB2)X(1B)a12t1tXGttj0把（11）式代入（10），比较两边B的同次幂的系数得G1，GGG,j（2）ARMA(2,1)模型实质上是一个二阶非齐次差分方程：XXXaa欲求其解，必须先求出其相应的齐次差分方程的通解。齐次方程对应的特征方程为2特征根为

-286-

jj

2

jj

2

121

ARMA(12,1)系统的格林函数为：

1221XXXXttt

t

X(Ij0

t1t12t2t

，，I22

0G

012Gcc1

121

221

ttt

X

j1221212

G11111

cc1112211

，c

1j21j12，

12c

2

c122121

12齐次差分方程的通解为Gccj112其中c和cG1

于是有Gcc

而，即

解之，得c

则G

3.3逆函数和可逆性前面的格林函数，把表示为过去a对的影响，或者说系统对过去a性，也就是用一个MA模型来逼近的行为。平稳序列X的这种表达形式称为的“传递形式”。同样我们也可以用过去的的一个线性组合来逼近系统现在时刻的行为。即XIXIXaIXat1t12t2tjtjt我们把这种表达形式称为的“逆转形式”。其中的系数函数I数。可见它是一个无穷阶的自回归模型。一个过程是否具有逆转形式，也就是说逆函数是否存在的性质，通常称为过程是否具有可逆性，如果一个过程可以用一个无限阶的自回归模型逼近，即逆函数存在，我们就称该过程具有可逆性，否则，就是不可逆的。对于AR（2）模型XXXa有II0,j122j可见，所谓可逆性，是指移动平均模型可以用AR模型表示。MA(1)模型：

-287-

t1t

X1tB1

t1t

X1tB1

I|1时，才有j0，故j1

1，拟合ARMA(2n,2n1)p3，拟和AR(p)模型。XXXXapIj

XXXaa1t12t2t1t1

jj1

I3II1,I2,I3替代I1,I2,I3是一种近似代替。通1的绝对值若大于1，则取其倒数作为初始值，以满足可逆性条件。I1,I2,I3及1，再用下式来确定ARMA(2,1)模型中的,III12121n1，拟合