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

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

2、的灰色預(yù)測(cè)模型就是關(guān)于數(shù)列預(yù)測(cè)的一個(gè)變量、一階微分的GM(1,1)模型。它是基于隨機(jī)的原始時(shí)間序列，經(jīng)按時(shí)間累加后所形成的新的時(shí)間序列呈現(xiàn)的規(guī)律可用一階線性微分方程的解來(lái)逼近。經(jīng)證明，經(jīng)一階線性微分方程的解逼近所揭示的原始時(shí)間序列呈指數(shù)變化規(guī)律。因此，當(dāng)原始時(shí)間序列隱含著指數(shù)變化規(guī)律時(shí)，灰色模型GM(1,1)的預(yù)測(cè)是非常成功的。1灰色系統(tǒng)的模型GM(1,1)1.1GM(1,1)的一般形式設(shè)有變量X(0)=X(0)(i),i=1,2,.,n為某一預(yù)測(cè)對(duì)象的非負(fù)單調(diào)原始數(shù)據(jù)列，為建立灰色預(yù)測(cè)模型：首先對(duì)X(0)進(jìn)行一次累加(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)模型。上述白化微分方程的解為(離散響應(yīng))：X、(k+1)=(X(0)(1)-)ekaX(k)=(X(0)(1)ujeWkf+u(4)aa式中：k為時(shí)間序列，可取年、季或月。1.2辯識(shí)算法記參數(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預(yù)測(cè)值的還原由于GM模型得到的是一次累加量，kWn+1,n+2,時(shí)刻的預(yù)測(cè)值，必須將GM模型所得數(shù)據(jù)X(1)(k+1)(或X(k)經(jīng)過(guò)逆生成，即累減生成(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因?yàn)閄(k-1)=£X(i),所以X(k)=X(k)X、(k-1)。i=12 應(yīng)用舉例取某高校1998年2003年的某專業(yè)招生數(shù)據(jù)建模，見(jiàn)表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對(duì)X(o)進(jìn)行一次累加(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)得灰色預(yù)測(cè)模型X(k)=(x(0)(1)=(132+277.0137483)e0.205(kJ)277.0137483=409.0137483e0.205(k-1)_277.0137483預(yù)測(cè)值及預(yù)測(cè)精度見(jiàn)表2。表2某高校專業(yè)招生預(yù)測(cè)值及預(yù)測(cè)精度表年GM(1,1)模型計(jì)算值1AGO還原值實(shí)際值誤差擬合相對(duì)誤差()2000132132132132002001225.06087962249392112002339.295441834211411843.382

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

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

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

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

11、，X?(0)為相應(yīng)的模擬誤差序列，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)為殘差萬(wàn)差，nk4.一s2._.1.稱c=&為均方差比值；對(duì)于給定的C0>0,當(dāng)CCC0時(shí)，稱模型為均方差比合格模Si型。2.稱p=pq*(k)司C0.6745S1)為小誤差概率，對(duì)于給定的P0>0,當(dāng)pAP0時(shí)，稱模型為小誤差概率合格模型。精度檢驗(yàn)等級(jí)參照表一-指際臨界性精度等級(jí)相對(duì)誤差關(guān)聯(lián)度均方差比

12、值小誤差概率一級(jí)10.010.900.350.95二級(jí)0.050.800.500.80三級(jí)0.100.700.650.70四級(jí)：0.200.600.800.60般情況下，最常用的是相對(duì)誤差檢驗(yàn)指標(biāo)。應(yīng)用舉例2、設(shè)原始序列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)模型，并進(jìn)行檢驗(yàn)。解：1)對(duì)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)對(duì)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及時(shí)間響應(yīng)式姆(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)誤差檢驗(yàn)序號(hào)實(shí)際數(shù)據(jù)模擬數(shù)據(jù)殘差相對(duì)誤差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平均相對(duì)誤差15二'4k4=1.0625%-0.09110.06621(1.41%0.51%2.69%1.80%)4計(jì)算X與必的灰色關(guān)聯(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精度為一級(jí)，可以用針(k1)=85.3728卻37156k-82.4986欠(k+1)=*(k+1)-?(k)預(yù)測(cè)。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