1、第十章 灰色模型介紹及應用 （徐利艷 天津農學院 2.4萬字） 10.1灰色理論基本知識 10.1.1概言 10.1.2有關名詞概念 10.1.3GM建模機理 10.2灰色理論模型應用 10.2.1GM（1，1）模型的應用污染物濃度問題 10.2.2 GM（1，1）殘差模型的應用油菜發病率問題 10.2.3 GM模型在復雜問題中的應用SARS 疫情問題 10.2.4 GM（1，n）模型的應用因素相關問題 本章小結思考題 推薦閱讀書目第十章 灰色模型介紹及應用10.1灰色理論基本知識10.1.1概言客觀世界的很多實際問題，其內部的結構、參數以及特征并未全部被人們了解，人們不可能象研究白箱問題那樣

2、將其內部機理研究清楚，只能依據某種思維邏輯與推斷來構造模型。對這類部分信息已知而部分信息未知的系統，我們稱之為灰色系統。本章介紹的方法是從灰色系統的本征灰色出發，研究在信息大量缺乏或紊亂的情況下，如何對實際問題進行分析和解決。灰色系統的研究對象是“部分信息已知、部分信息未知”的“小樣本”、“貧信息”不確定性系統，它通過對“部分”已知信息的生成、開發實現對現實世界的確切描述和認識。信息不完全是“灰”的基本含義。灰色系統理論建模的主要任務是根據具體灰色系統的行為特征數據，充分開發并利用不多的數據中的顯信息和隱信息，尋找因素間或因素本身的數學關系。通常的辦法是采用離散模型，建立一個按時間作逐段分析的

3、模型。但是，離散模型只能對客觀系統的發展做短期分析，適應不了從現在起做較長遠的分析、規劃、決策的要求。盡管連續系統的離散近似模型對許多工程應用來講是有用的，但在某些研究領域中，人們卻常常希望使用微分方程模型。事實上，微分方程的系統描述了我們所希望辨識的系統內部的物理或化學過程的本質。目前，灰色系統理論已成功地應用于工程控制、經濟管理、未來學研究、生態系統及復雜多變的農業系統中，并取得了可喜的成就。灰色系統理論有可能對社會、經濟等抽象系統進行分析、建模、預測、決策和控制，它有可能成為人們認識客觀系統改造客觀系統的一個新型的理論工具。10.1.2有關名詞概念灰數：一個信息不完全的數，稱為灰數。灰元

4、：信息不完全或內容難以窮盡的元素，稱為灰元。灰關系：信息不完全或機制不明確的關系，稱為灰關系。具有灰關系的因素是灰因素，灰因素之間的量化作用，稱為灰關聯。灰色系統:含灰數、灰元或灰關系的系統稱為信息不完全系統。如果按照灰色理論去研究它。則稱此系統為灰色系統。累加生成：由于灰系統對一切隨機量都可看作是在一定范圍內變化的灰色量，因此，為適應灰系統建模需要，提出“生成”的概念，“生成”即指對原始數據做累加（或累減）處理。累加生成一般可寫成AGO。若計為原始數列，為次累加生成后數列，即則次累加生成算式為 一般常用的是一次累加生成，即10.1.3GM建模機理建立GM模型，實際就是將原始數列經過累加生成后

5、，建立具有微分、差分近似指數規律兼容的方程，成為灰色建模，所建模型稱為灰色模型，簡記為GM（Grey Model）。如GM（m,n）稱為m階n個變量的灰色模型，其中GM（1，1）模型是GM（1，n）模型的特例，是灰色系統最基本的模型，也是常用的預測模型，因此本章重點介紹幾種GM（1，1）模型的建模過程和計算方法，并簡單介紹GM（1，n）建模過程。GM（1，1）的建模機理GM（1，1）模型是GM（1，N）模型的特例，其簡單的微分方程形式（白化形式的微分方程）是利用常數變易法解得，通解為若初始條件為，則可得到微分方程的特解為或時間響應函數其中白化微分方程中的項中的為的背景值，也稱為初始值; 為常數

6、（有時也將寫成）。按白化導數定義有差分形式的微分方程，即顯然，當時間密化值定義為1，即當時，上式可記為記為離散形式這顯然表明是一次累計生成，因此上述方程可改寫為這實際也表明，模型是以生成數（是以的一次累加）為基礎的。當足夠小時，到不會發生突變，因此可取與的平均值作為時的背景值，因此，背景值便可記為或于是白化的微分方程可改寫為或即因此，上述方程可以改寫為矩陣方程形式，即引入下列符號，設 于是便有令 則解得將求解得到的代入微分方程的解式（也稱時間響應函數），則由于，因此求導還原得上述兩式便為GM（1，1）的時間響應式，及灰色系統預測模型的基本算式，當然上述兩式計算結果只是近似計算值。為簡記，一般可

7、以將GM（1，1）的建模過程記為10.2灰色理論模型應用10.2.1GM（1，1）模型的應用污染物濃度問題GM（1，1）模型是灰色系統最基本的模型，下面以污染物濃度問題說明GM（1，1）模型的建立及求解過程。例10.1 某污染源中某種污染物質量濃度測量值如表10.1，試建立GM（1，1）模型 表10.1 某污染物質量濃度測量值 （mg/L） 年 份2001200220032004200520063.9364.5754.9685.0635.9685.507解：第一步，設原始數據為 第二步，對原始數據進行累加生成，即因此累加生成數據為第三步，構造矩陣 第四步，計算。先求，即根據逆矩陣的求解方法，得

8、再求的值，即進而求得的值為計算GM1_1的程序如下function 10toliti01(X0)m,n=size(X0);X1=cumsum(X0); X2=;for i=1:n-1 X2(i,:)=X1(i)+X1(i+1);endB=-0.5.*X2;t=ones(n-1,1);B=B,t; YN=X0(2:end);P_t=YN./X1(1:(length(X0)-1)A=inv(B.'*B)*B.'*YN.'a=A(1) u=A(2) Bb1=B.'*Bb2=inv(B.'*B)b3=B.'*YN.'b4=u/ab5=X1(1)-

9、b4b6=-a*b5第五步，將的值代入微分方程的時間響應函數，令，得第六步，求導還原得第七步，對上述模型進行精度檢驗。常用的方法是回代檢驗，即分別用模型求出各時刻值，然后求相對誤差。先利用時間響應函數模型求各時刻值（），并計算相對誤差，結果如表10.2所示. 表10.2 精度檢驗實測值、殘差值表 GM計算值實測值殘差相對殘差8.605913.534418.735924.225430.01908.5110 13.4790 18.5420 24.5100 30.0170-0.0949 -0.0554 -0.1939 0.2846 -0.0020-0.0112-0.0041-0.01050.0116

10、-0.0001再利用時間響應函數模型求各時刻值（），并計算相對誤差，結果如表10.3所示. 表10.3 計算值與實驗原始數據值對照表 GM計算值實測值殘差相對殘差4.79605.06165.34195.63775.94994.57504.96805.06305.96805.5070-0.2210-0.0936-0.27890.3303-0.4429-0.0483 -0.0188 -0.0551 0.0553-0.0804從殘差檢驗結果看，累計生成數列曲線擬合較好，相對誤差在0.01即1%左右；而還原數列的相對誤差較大，其原因是累加生成數據將原始數據的隨機性弱化，正負誤差有抵消的，當數據再被還原

11、回來時便表現出來。10.2.2 GM（1，1）殘差模型的應用油菜發病率問題當（，）模型的精度不符合要求時，可用殘差序列建立（，）模型，對原來的模型進行修正，以提高精度，即建立殘差（，）模型，步驟如下第一步，利用原始數據建立GM（1，1）模型，得時間響應式其中第二個式子也成為導數還原值。鑒于導數還原值與原始數據（累減還原值）不一致，為減少往返運算造成的誤差，往往用原始數據與導數還原值的殘差修正的模擬值。第二步，利用殘差數列建立新的GM（1，1）模型。建立殘差模型的過程和計算方法同于GM（1，1）建模過程，只不過建立殘差模型所用的原始數列采用的是殘差數據。令為殘差，則即或利用殘差序列建立新的GM（

12、1，1）模型，求解得時間響應式第三步，結合上兩步的GM（1，1）模型，建立殘差GM（1，1）模型結合上兩步的GM（1，1）模型，則相應的殘差修正時間響應式為稱為導數還原式的殘差修正模型。例10.2 某縣油菜發病率數據如表10.4所示，試建立殘差GM模型并進行求解。 表10.4 某縣油菜發病率數據 （%） 序 號12345678910111213620402540453521141815.51715解：第一步，建立原始數據的GM（1，1）模型設原始數據為建立GM（1，1）模型，利用GM（1，1）的求解程序得時間響應式為第二步，誤差檢驗利用時間響應函數模型計算各時刻值（），并計算相對誤差，程序如下

13、function 10toliti02(X0)%format long ;%X0=0.01*6 20 40 25 40 45 35 21 14 18 15.5 17 15;m,n=size(X0);s(1)=1;for i=1:12 y(i+1)=0.368*exp(-0.06486*i); z(i+1)=X0(i+1)-y(i+1); w(i+1)=z(i+1)/X0(i+1); s(i+1)=i+1;endy'X0'z'w'z*z'sum(abs(w)/12計算結果如表10.5所示 表10.5 計算值與實驗原始數據值對照表 GM計算值實測值殘差相對殘

14、差0.34490.32320.30290.28390.26610.24940.23370.21900.20530.19240.18030.16900.20.40.250.40.450.350.210.140.180.1550.170.15-0.14490.0768-0.05290.11610.18390.1006-0.0237-0.0790-0.0253-0.0374-0.0103-0.0190-0.72440.1919-0.21170.29020.40870.2875-0.1129-0.5645-0.1404-0.2412-0.0606-0.1265 由表可以看出，最大誤差高達72.44%，

15、最低的也達到6.06%，模擬誤差較大，進一步計算平均相對誤差 平均相對誤差很較大，相對精度約70%。因此為了提高遠原點（即現時）精度，即將最后一個誤差減小，需采用殘差模型進行修正。第三步，以部分殘差數據為原始數據建立新的GM（1，1）模型取得殘差尾端，即取最后5個數據的殘差：-0.0790,-0.0253,-0.0374,-0.0103,-0.0190,用此尾段可建立殘差尾段模型，取絕對值，得殘差數列以上述的殘差數列為原始數據建立新的GM（1，1）模型，得殘差的時間響應式第四步，將原始數據和部分殘差數據的兩個GM（1，1）模型即和結合，得到修正后的殘差GM（1，1）模型第五步，用修正后的模型對

16、的模擬值進行修正，結果為：第六步，精度檢驗建立如下程序：function 10toliti021(X0)%format long ;%X0=0.01*6 20 40 25 40 45 35 21 14 18 15.5 17 15;m,n=size(X0);s(1)=1;for i=8:12 y(i+1)=0.368*exp(-0.06486*i)-0.0328*exp(-0.1894*i); z(i+1)=X0(i+1)-y(i+1); w(i+1)=z(i+1)/X0(i+1); s(i+1)=i+1;endy'X0'z'w'z*z'sum(abs(w

17、)/5計算結果如表10.6所示 表10.6 修正后計算值與實驗原始數據值檢驗結果 GM計算值實測值殘差相對殘差0.21180.19930.18740.17620.16560.140.180.1550.170.15-0.0718-0.0193-0.0324-0.0062-0.0156-0.5130-0.1073-0.2093-0.0366-0.1040按此模型，可對五個模擬值進行修正，修正后的平均相對誤差，精度有明顯的提高。尤其對于原點附近的兩個數據：0.17，0.15相對誤差分別降低為3.66%和10.4%，低于允許誤差要求。這說明，對原點數據GM（1，1）模型修正是有必要的。10.2.3GM

18、模型在復雜問題中的應用SARS 疫情問題例10.3 SARS 疫情問題2003年的SARS疫情對我國的發展產生了一定影響，尤其在經濟發展方面產生了很大的影響，下面就SARS疫情對我國經濟的影響問題建立GM模型并求解。10.2.3.1 問題的提出 2003年的SARS 疫情對中國部分行業的經濟發展產生了一定影響，特別是對部分疫情較嚴重的省市的相關行業所造成的影響是顯著的，經濟影響主要分為直接經濟影響和間接影響。直接經濟影響涉及商品零售業、旅游業、綜合服務等行業。很多方面難以進行定量地評估，現僅就SARS 疫情較重的某市商品零售業、旅游業和綜合服務業的影響進行定量的評估分析。究竟 SARS 疫情對

19、商品零售業、旅游業和綜合服務業的影響有多大，已知某市從1997 年1 月到2003 年12 月的商品零售額、接待旅游人數和綜合服務收入的統計數據如表10.710.9所示 表10.7 商品的零售額 （億元） 月份年份1月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月199719981999200020012002200383.0 79.8 78.1 85.1 86.6 88.2 90.3 86.7 93.3 92.5 90.9 96.9101.7 85.1 87.8 91.6 93.4 94.5 97.4 99.5 104.2 102.3 101.

20、0 123.592.2 114.0 93.3 101.0 103.5 105.2 109.5 109.2 109.6 111.2 121.7 131.3105.0 125.7 106.6 116.0 117.6 118.0 121.7 118.7 120.2 127.8 121.8 121.9139.3 129.5 122.5 124.5 135.7 130.8 138.7 133.7 136.8 138.9 129.6 133.7137.5 135.3 133.0 133.4 142.8 141.6 142.9 147.3 159.6 162.1 153.5 155.9163.2 159.7

21、 158.4 145.2 124.0 144.1 157.0 162.6 171.8 180.7 173.5 176.5 表10.8 接待海外旅游人數 （萬人）月份年份1月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月19971998199920002001200220039.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.69.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.910.1 12.9 17.7 21.0 21.

22、0 20.4 21.9 25.8 29.3 29.8 23.6 16.511.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.511.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.713.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.915.4 17.1 23.5 11.6 1.78 2.61 8.8 16.2 20.1 24.9 26.5 21.8 表10.9 綜合服務業累計數額 （億元）月份年份2 月

23、3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月199719981999200020012002200396 144 194 276 383 466 554 652 747 832 972111 169 235 400 459 565 695 805 881 1011 1139151 238 335 425 541 641 739 866 975 1087 1238164 263 376 531 600 711 913 1038 1173 1296 1497182 318 445 576 708 856 1000 1145 1292 1435 1667216

24、 361 504 642 818 979 1142 1305 1479 1644 1920241 404 584 741 923 1114 1298 1492 1684 1885 2218 試根據這些歷史數據建立預測評估模型，評估 2003 年SARS 疫情給該市的商品零售業、旅游業和綜合服務業所造成的影響。10.2.3.2模型的假設與分析模型假設：（1）假設該市的統計數據都是可靠準確的；（2）假設該市在SARS 疫情流行期間和結束之后，數據的變化只與SARS 疫情的影響有關，不考慮其它隨機因素的影響。模型分析：根據所掌握的歷史統計數據可以看出，在正常情況下，全年的平均值較好地反映了相關指標的

25、變化規律，這樣可以把預測評估分成兩部分：（1）利用灰色理論建立GM(1,1)模型，由19972002 年的平均值預測2003 年平均值；（2）通過歷史數據計算每個月的指標值與全年總值的關系，從而可預測出正常情況下2003 年每個月的指標值，再與實際值比較可以估算出SARS 疫情實際造成的影響。10.2.3.3建立灰色預測模型GM(1,1)第一步,數據處理（1）原始數據：根據表中的已知數據，計算19972002年某項指標的年平均值，作為原始數據，記為并要求級比（2）數據的累加生成：對原始數據進行一次累加生成，因此累加生成數據記為（2）背景值的選擇：取累加生成數據的加權平均值為背景值，即其中為確定

26、參數。背景值記為第二步，GM（1，1）模型的建立（1）建立GM（1，1）的白化微分方程模型其中是a發展灰度，u是內生控制灰度。（2）轉化為灰微分方程或即矩陣形式為其中(3)轉化為時間響應函數利用最小二乘法得到參數的估計值，進而得到灰微分方程的解，對求導還原得。即參數的估計值為微分方程的解式（也稱時間響應函數）為其中，稱為還原值。第三步，利用模型預測指標值根據時間響應函數可以預測出正常情況下2003年的平均值，則，則預測2003年的總值為。根據歷史數據，可以計算出2003 年第i個月的指標值占全年總值的比例為，即記為，于是可以可到2003年每一個月的指標值10.2.4.4模型求解及結果分析（1）

27、商品零售額根據商品零售額的數據表，計算得年平均值（即原始數據）和一次累加生成值，分別為顯然的所有級比都在可容區域內，這里取，計算可得背景值計算得參數的估計值為，進而得到時間響應函數再根據時間響應函數預測可得，2003年的月平均值為億元；年總值為億元。又根據比例的表達式計算出每月的比例為 因此2003 年112 月的預測值(單位：億元)為將預測值與實際值進行比較，結果如表10.10所示 表10.10 2003年商品的零售額比較表 （億元） 月份1月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月實際值預測值163.2 159.7 158.4 145.

28、2 124.0 144.1 157.0 162.6 171.8 180.7 173.5 176.5152.9 155.3 144.2 151.2 157.7 157.4 162.6 161.3 168.0 170.5 166.7 177.1 圖形如圖10-1：（藍線為實際值，紅線為預測值）圖10-1 2003年商品的零售額實際值與預測值比較圖通過圖形可以直觀的看出：（1）預測值波動比較小，真實值波動比較劇烈；（2）5月份左右真實值遠遠低于預測值，年初和年末都高于預測值。這是由實際情況造成的，年初當SARS 疫情剛剛開始的時候，人們儲備保健藥品和保健食物等，拉動了零售額的增長；5月份左右，SAR

29、S 疫情比較猖獗，此時好多學校和單位等實行封閉管理，大大限制了人們的消費，因此零售額明顯降低；年末SARS 疫情慢慢遠去，此前被限制的消費得以充分實現，又促進了零售額的增長。當然可以根據模型所得數據，對SARS 疫情給該市的商品零售業造成的影響進行定量分析，這里不再詳述。計算的MATLAB 程序如下：function 10toliti03clcclc,clearload shuju1.txt %把原始數據保存在純文本文件shuju1.txt中han1(end,:)=;x0=mean(han1,2)m=size(han1,2);n=size(x0,1);z1=;x1=cumsum(x0)alph

30、a=0.5;for i=1:n-1z1(i,:)=(1-alpha)*x1(i)+alpha*x1(i+1);endz1Y=x0(2:n);B=-z1,ones(n-1,1);A=inv(B'*B)*B'*Y;a=A(1)u=A(2)b4=u/ab5=x1(1)-b4b6=-a*b5k=6;x7hat=(x0(1)-u/a)*(exp(-a*k)-exp(-a*(k-1)z=m*x7hatu=sum(han1)/sum(sum(han1)v=z*u（2）接待海外旅游人數根據接待海外旅游人數的數據表，計算得年平均值（即原始數據）和一次累加生成值，分別為顯然的所有級比都在可容區域內

31、，這里取，計算可得背景值計算得參數的估計值為，進而得到時間響應函數再根據時間響應函數預測可得，2003年的月平均值為萬人；年總值為萬人。又根據比例的表達式計算出每月的比例為 因此2003 年112 月的預測值(單位：萬人)為將預測值與實際值進行比較，結果如表10.11所示 表10.11 2003年接待海外旅游人數 （萬人） 月份1月 2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月實際值預測值15.4 17.1 23.5 11.6 1.78 2.61 8.8 16.2 20.1 24.9 26.5 21.814.8 26.6 25.5 31.9 33

32、.0 30.8 30.4 37.1 36.7 37.8 33.2 25.5 圖形如圖10-2：（藍線為實際值，紅線為預測值）圖10-2 2003年接待海外旅游人數實際值與預測值比較圖通過圖形可以直觀的看出：（1）預測值波動比較小，真實值波動比較劇烈；（2）真實值低于預測值，尤其5月份左右真實值遠遠低于預測值，年初和年末相差不太大。這是由實際情況造成的， 5月份左右SARS 疫情比較猖獗，此時好多學校和單位等實行封閉管理，大大限制了人們的出行，同時人們也基于自身安全的因素，能不出門就不出門，因此旅游人數大大降低，旅游業處于低谷；年初和年末SARS 疫情對人們出行的影響不大，因此年初和年末年末海外

33、旅游人數的實值略低于預測值。（3）綜合服務業累計數據根據綜合服務業累計數據的數據表，計算得年平均值（即原始數據）和一次累加生成值，分別為顯然的所有級比都在可容區域內，這里取，計算可得背景值計算得參數的估計值為。進而得到時間響應函數再根據時間響應函數預測可得，2003年的月平均值為億元；年總值為億元。又根據比例的表達式計算出每月的比例為 因此2003 年112 月的預測值(單位：億元)為將預測值與實際值進行比較，結果如表10.12所示 表10.12 2003年綜合服務業累計數據的比較表 （億元） 月份2 月 3 月 4 月 5 月 6 月 7 月 8 月 9 月 10 月 11 月 12 月實際

34、值預測值241 404 584 741 923 1114 1298 1492 1684 1885 2218240 390 545 744 916 1101 1316 1517 1709 1907 2201圖形如圖10-3：（藍線為實際值，紅線為預測值）圖10-3 2003年綜合服務業累計數據實際值與預測值比較圖通過圖形可以直觀的看出：（1）預測值與真實值相差不大；（2）5月份左右真實值與預測值最接近。這是由實際情況造成的，SARS 疫情對于綜合服務業中的部分行業影響較大，如航空交通運輸、賓館餐飲等，但有些行業影響不大，如電信、通訊等，因此總平均來看，影響還不算太大。（4）模型的評價從三方面的結

35、果分析，可以看出模型的結論與實際情況相符，這說明了模型的正確性和可靠性。雖然該模型是就某經濟指標的發展規律進行評估預測而建立的，但類似地也適用于其它方面的一些數據規律的評估預測問題，即該模型具有很廣泛的應用性。10.2.4 GM（1，n）模型的應用因素相關問題GM（1，n）模型表示一階的含有n個變量灰色模型，適合于建立系統的狀態模型與各變量的動態分析，與GM（1，1）模型不同，不適合預測用；但建模與計算過程與GM（1，1）模型類似。GM（1，2）模型是GM（1，n）模型的基礎，因此下面以GM（1，2）模型為例說明GM（1，n）模型的建模過程。GM（1，2）模型表示一階的含有兩個變量灰色模型，其

36、相應的白化微分方程模型為時間相應函數(即解)為還原值為例10.4 某系統中兩因素和相互之間存在關系，其中因素為系統特征因素，因素為相關因素，兩因素的關系如表10.13所示表10.13 兩因素的關系 序號1 2 3 4 5因素因素10.7 12.3 15.4 19.2 21.5 40.3 50.5 60.3 65.5 75.2解：第一步，數據處理取原始數據為對原始數據累加生成，得第二步，建立GM（1，2）模型白化的微分方程模型為轉化為灰微分方程或其中。即矩陣形式為其中(3)時間響應函數利用最小二乘法得到參數的估計值，即進而得到微分方程的解即時間響應函數及還原值，即 第三步，模型的求解建立M文件1

37、0toliti04.m如下：clc,clearx10=10.7,12.3,15.4,19.2,21.5;x20=40.3,50.5,60.3,65.5,75.2;n=length(x10);x11=cumsum(x10)x21=cumsum(x20)for i=2:nz11(i)=0.5*(x11(i)+x11(i-1);endB=-z11(2:n)',x21(2:n)'Y=x10(2:n)'A=inv(B'*B)*B'*Yx=dsolve('Dx+a*x=b*x2','x(0)=x0');x=subs(x,'a','b','x0','x2',A(1),A(2),x10(1),'x21');digits(6),x=vpa(x);x=simple(x)x=subs(x,'t','x21',0:n-1,x21(1:n)xha