灰色系統(tǒng)預測GM1,1模型及其Matlab實現(xiàn)

1、灰色系統(tǒng)預測GM(1,1)模型及其Matlab實現(xiàn)預備知識(i)灰色系統(tǒng)白色系統(tǒng)是指系統(tǒng)內(nèi)部特征是完全已知的；黑色系統(tǒng)是指系統(tǒng)內(nèi)部信息完全未知的；而灰色系統(tǒng)是介于白色系統(tǒng)和黑色系統(tǒng)之間的一種系統(tǒng)，灰色系統(tǒng)其內(nèi)部一部分信息已知，另一部分信息未知或不確定。(2)灰色預測灰色預測，是指對系統(tǒng)行為特征值的發(fā)展變化進行的預測，對既含有已知信息又含有不確定信息的系統(tǒng)進行的預測，也就是對在一定范圍內(nèi)變化的、與時間序列有關的灰過程進行預測。盡管灰過程中所顯示的現(xiàn)象是隨機的、雜亂無章的，但畢竟是有序的、有界的，因此得到的數(shù)據(jù)集合具備潛在的規(guī)律?；疑A測是利用這種規(guī)律建立灰色模型對灰色系統(tǒng)進行預測。目前使用最廣泛

2、的灰色預測模型就是關于數(shù)列預測的一個變量、一階微分的GM(1,1)模型。它是基于隨機的原始時間序列，經(jīng)按時間累加后所形成的新的時間序列呈現(xiàn)的規(guī)律可用一階線性微分方程的解來逼近。經(jīng)證明，經(jīng)一階線性微分方程的解逼近所揭示的原始時間序列呈指數(shù)變化規(guī)律。因此，當原始時間序列隱含著指數(shù)變化規(guī)律時，灰色模型GM(1,1)的預測是非常成功的。1灰色系統(tǒng)的模型GM(1,1)1.1GM(1,1)的一般形式設有變量X(0)=X(0)(i),i=1,2,.,n為某一預測對象的非負單調(diào)原始數(shù)據(jù)列，為建立灰色預測模型：首先對X(0)進行一次累加(1AGO,AcumulatedGeneratingOperator)生成一

3、次累加序列：X=X(1)(k),k=1,2,，n其中X(k)=£X(0)(i)i1=X(1)(k-1)+X(k)XX可建立下述白化形式的微分方程：1 1)(2)d+aX(1)=udt即GM(1,1)模型。上述白化微分方程的解為(離散響應)：X、(k+1)=(X(0)(1)-)ekaX(k)=(X(0)(1)ujeWkf+u(4)aa式中：k為時間序列，可取年、季或月。1.2辯識算法記參數(shù)序列為a,a=a,uT,a'可用下式求解：a=(BTB)-1BTYn式中：b數(shù)據(jù)陣；Yn一數(shù)據(jù)列_L(X(1)(1)十X(2)1(6)2 B=1- 2(X(2)+X(3)1.- -1(X(1)

4、(n-1)+X(n)1.2(7)- n=(X(0)(2),X(0)(3),X(0)(1.3預測值的還原由于GM模型得到的是一次累加量，kWn+1,n+2,時刻的預測值，必須將GM模型所得數(shù)據(jù)X(1)(k+1)(或X(k)經(jīng)過逆生成，即累減生成(IAGO還原為X(0)(k+1)(或X(0)(k)才能用。X(1)(k)=ZX(0)(i)1 =1k1=EX(0)(i)+X(0)(k)i1X(k)=X(1)(k)9X(0)(i)i1因為X(k-1)=£X(i),所以X(k)=X(k)X、(k-1)。i=12 應用舉例取某高校1998年2003年的某專業(yè)招生數(shù)據(jù)建模，見表1。表1某高校專業(yè)招生

5、數(shù)據(jù)表年招生人數(shù)20001322001922002118200313020041872005207以表1中的數(shù)據(jù)構(gòu)造原始數(shù)據(jù)列x(0),即X=X(0)(1),X(0)(2),X(0)(3),X(0)(4),X(0)(5),X(6)=132,92,118,130,187,207對X(o)進行一次累加(1AGO),生成數(shù)列：kx(1)(k)=XX(0)(i)即x(1)=X(1)(1),X(1)(2),X(1)(3),X(1)(4),X(1)(5),X(6)=132,224,342,472,659,866和數(shù)據(jù)陣日數(shù)據(jù)列Y|一1781,Yn=(92,118,130,187,207)D_2831B=1

6、-4071-565.51'-762.51由式得Ta=a,u0.205一IGM(1,1)為56.7878184ua(k二)ea由式(4)得灰色預測模型X(k)=(x(0)(1)=(132+277.0137483)e0.205(kJ)277.0137483=409.0137483e0.205(k-1)_277.0137483預測值及預測精度見表2。表2某高校專業(yè)招生預測值及預測精度表年GM(1,1)模型計算值1AGO還原值實際值誤差擬合相對誤差()2000132132132132002001225.06087962249392112002339.295441834211411843.382

7、003479.52123472140130107.692004651.65196591721871582005862.9466129866211207一41.9320061122.316167259252一72.78由表2知預測精度較高。2006年某專業(yè)招生人數(shù)預測值為259人。由于人數(shù)為整數(shù),所以結(jié)果取整數(shù)部分。GM(1,1)是一種長期預測模型，在沒有大的市場波動及政策性變化的前提下，該預測值應是可信的。如前所述，影響招生人數(shù)的因素很多且難以預測。因此，在采用灰色系統(tǒng)理論進行定量預測時，如果存在對預測對象影響較大的因素，就要在定性分析的基礎上，尋找原始數(shù)據(jù)信息的突變點的量化值，然后再對預測值

8、進行必要的修正，使預測值更接近實際情況，提高預測值的可信度，為科學決策提供可靠的數(shù)據(jù)。另外，若作長期預測，要考慮對上限值的約束條件。應用灰色預測模型GM(1,1),對某專業(yè)招生人數(shù)進行了預測，具有較高的預測精度。應用灰色模型進行預測較之其它常規(guī)的時間序列預測法有以下顯著的特點。(1) 灰色模型是一種長期預測模型，將預測系統(tǒng)中的隨機元素作為灰色數(shù)據(jù)進行處理，而找出數(shù)據(jù)的內(nèi)在規(guī)律。進行預測所需原始數(shù)據(jù)量小，預測精度較高，無須像其它預測法要么需要數(shù)據(jù)量大且規(guī)律性強，要么需要憑經(jīng)驗給出系數(shù)。(2) 理論性強，計算方便，籍助計算機及其程序設計語言，使得數(shù)據(jù)處理簡便、快速、準確性好。(3) 用有限的表征系

9、統(tǒng)行為特征的外部元素，分析系統(tǒng)的內(nèi)在規(guī)律。灰色系統(tǒng)理論采用對系統(tǒng)的行為特征數(shù)據(jù)進行生成的方法，對雜亂無章的系統(tǒng)的行為特征數(shù)據(jù)進行處理，從雜亂無章的現(xiàn)像中發(fā)現(xiàn)系統(tǒng)的內(nèi)在規(guī)律，這是該方法的獨特之處。(4) 適用性強。用灰色模型既可對周期性變化的系統(tǒng)行為進行預測，亦可對非周期性變化的系統(tǒng)行為進行預測；既可進行宏觀長期的預測，亦可用于微觀短期的預測。2灰色系統(tǒng)模型的檢驗定義1.設原始序列X(0)(0)(1),x(0)(2),x(0)(n)相應的模型模擬序列為,(0)=3(1),?(2),外0)(n)：'殘差序列；(0)=；(1),；(2),；(n”x(0)(1)->?(0)(1),x(0

10、)(2)一/)(2),x(0)(n)-姆)相對誤差序列、x(0)(1)x(0)(2)'；(n)x(n)1.對于kvn,稱Ak；(k)x(0)(k)為k點模擬相對誤差，稱An=，一1n，一，*，*稱=1£Ak為平均模擬相對誤差;k2.3.；(n)x(0)(n)為濾波相對誤差,nkd稱1-為平均相對精度，1-An為濾波精度；給定口，當2<“,且An<ot成立時，稱模型為殘差合格模型。定義2設X(0)為原始序列，義(0)為相應的模擬誤差序列，注為X(0)與父(0)的絕對關聯(lián)度,若對于給定的zo>Q,z>0o,則稱模型為關聯(lián)合格模型。定義3設X(0)為原始序列

11、，X?(0)為相應的模擬誤差序列，0)為殘差序列。1n(H)一(H)x=_£x(0)(k)為X(0"q均值，nk4ns2=l工(x(k)_x)2為x的方差，nk4_1八、(=-1lw(k)為殘差均值，n21J-2,一、*S2=一£(s(k)-s)為殘差萬差，nk4.一s2._.1.稱c=&為均方差比值；對于給定的C0>0,當CCC0時，稱模型為均方差比合格模Si型。2.稱p=pq*(k)司C0.6745S1)為小誤差概率，對于給定的P0>0,當pAP0時，稱模型為小誤差概率合格模型。精度檢驗等級參照表一-指際臨界性精度等級相對誤差關聯(lián)度均方差比

12、值小誤差概率一級10.010.900.350.95二級0.050.800.500.80三級0.100.700.650.70四級：0.200.600.800.60般情況下，最常用的是相對誤差檢驗指標。應用舉例2、設原始序列X(0)1x(0)(1),x(0)(2),x(0)(3),x(0)(4),x(0)(5)=2.874,3.278,3.337,3.390,3.679建立GM(1,1)模型，并進行檢驗。解：1)對X(0)作1-AGQ得D為X(0)的一次累加生成算子，記為1-AGOX(1)_1(1)(1)x(2)x(3)x(1)(4)x(1)(5).JX-x(%x(2),x(3),x(4),x(5

13、)=2.874,6.152,9.489,12.579,16.5582)對X(1)作緊鄰均值生成，令Z(k)=0.5x(1)(k)0.5x(1)(k-1)Z(1)='z(1)(1),z(1)(2),z(1)(3),z(1)(4),z(5)=：2874,4.513,7.820,11.84,14.718于是，一-z(2)_z-z-z(4)1111I-4.5131I-7.8201-11.841'-14.7181"x(0)(2)"x(0)x(0)(4).x(0)(5)J3.27813.3373.3903.679rTr-4.513BB=|1-7.820-11.184-1

14、4.718111-4.5131-7.8201-11.841-14.7181-423.221-38.2351=1-38.2354(BB)J-423.2211-38.235-38.235=0.0173184-|0.16655420.1655421.832371=1一438.2351423.221x4-38.2352|I38-235423.221_1_438.2351=230.969138.235423.221J3.278-14.71813.3371*3.390-3.679i?=(BB),BY0.0173180.16554214.513-7.820-11.184|0.16655421.8323711

15、110.0873860.030115Jl.0852800.5378333.278-0.028143-0.089344|3.337-0.01905110.604076*13.390-3.679_一0.0371561一I3.065318_3)確定模型dx-0.037156x()=3.065318dt及時間響應式姆(k1)=(x(0)(1)-)ek-baa=85.3728e0.037156k-82.49864)求X的模擬值潭二京必共曲必)=(2.8740,6.1058,9.4599,12.9410,16.5538)5)還原出X(0)的模擬值，由婷(k1)=，(k1)-，)(k)得*(0)=歹0)(1

16、),00)(2),。0)(3),?(0)(4),。0)(5)=(2.8740,3.2318,3.3541,3.4811,3.6128)6)誤差檢驗序號實際數(shù)據(jù)模擬數(shù)據(jù)殘差相對誤差x(0)(k)x(0)(k)«k)=x(0)(k)-?(0)(k)k=鞏k)(0)/ix(k)23.2783.23180.04621.41%33.3373.3541-0.01710.51%43.3903.4811-0.09112.69%53.6793.61280.06621.80%殘差平方和H2)1武3）名（4）45）】斐）；K)?(5)-0.0462=0.0462-0.0171-0.09110.0662】1

17、-0.0171=0.0151085平均相對誤差15二'4k4=1.0625%-0.09110.06621(1.41%0.51%2.69%1.80%)4計算X與必的灰色關聯(lián)度二1S=£(x(k)-x(1)+-(x(5)-x(1)kz221小一八一、=(3.278-2.874)+(3.337-2.874)+(3.390-2.874)十一(3.6792.874)2=0.4040.4630.5160.4025=1.7855二八八1八八£(x(k)x(1)+(x(5)欠(1)y21=(3.2318-2.874)+(3.3541-2.874)+(3.4811-2.874)+-(

18、3.6128-2.874)=0.35780.48010.60710.3694=1.8144。一r一一11r一一nS=£(x(k)-x(1)-(?(k)-洌1)】+1(x(5)-x(1)-(?(5)-41)1=(0.3578-0.404)+(0.4801-0.463)+(0.60710.516)+(0.36940.4025)2=-0.0462+0.0171+0.091-0.01655=0.04535右_1+S+§_1+1.7855+1.8144_4.59991+|S|+同+§_S1+1.7855+1.8144+0.045354.64525=0.9902>0.9

19、0精度為一級，可以用針(k1)=85.3728卻37156k-82.4986欠(k+1)=*(k+1)-?(k)預測。functiony,p,e=gm_1_1(X,k)%graymodel:GM(1,1)%Exampley,p=gm_1_1(200250300350,2)ifnargout>3,error('Toomanyoutputargument.');endifnargin=1,k=1;x_orig=X;elseifnargin=0|nargin>2error('Wrongnumberofinputarguments.');endx_orig=

20、X;predict=k;%AGOprocessx=cumsum(x_orig);%computethecoefficient(aandu)n=length(x_orig);%firstgeneratethematrixBfori=1:(n-1);B(i)=-(x(i)+x(i+1)/2;endB=B'ones(n-1,1);%thengeneratethematrixYfori=1:(n-1);y(i)=x_orig(i+1);endY=y'%getthecoefficient.a=au(1)u=au(2)au=(inv(B'*B)*(B'*Y);%changethegreymodeltosymbolicexpressioncoef1=au(2)/au(1);coef2=x_orig(1)-coef1;coef3=0-au(1);costr1=num2str(coef1);costr2=num2str(abs(coef2);costr3=num2str(coef3);eq=strcat