灰色系统预测
重点内容：灰色系统理论的产生和发展动态，灰色系统的基本概念，灰色系统与模糊数学、黑箱方法的区别，灰色系统预测GM （1，1）模型，GM （1，N ）模型，灰色系统模型的检验，应用举例。
1灰色系统理论的产生和发展动态
1982邓聚龙发表第一篇中文论文《灰色控制系统》标志着灰色系统这一学科诞生。
1985灰色系统研究会成立，灰色系统相关研究发展迅速。 1989海洋出版社出版英文版《灰色系统论文集》，同年，英文版国际刊物《灰色系统》杂志正式创刊。目前，国际、国内200多种期刊发表灰色系统论文，许多国际会议把灰色系统列为讨论专题。国际著名检索已检索我国学者的灰色系统论著500多次。灰色系统理论已应用范围已拓展到工业、农业、社会、经济、能源、地质、石油等众多科学领域，成功地解决了生产、生活和科学研究中的大量实际问题，取得了显著成果。
2灰色系统的基本原理
2.1灰色系统的基本概念
我们将信息完全明确的系统称为白色系统，信息未知的系统称为黑色系统，部分信息明确、部分信息不明确的系统称为灰色系统。 系统信息不完全的情况有以下四种： 1. 元素信息不完全 2. 结构信息不完全 3. 边界信息不完全 4. 运行行为信息不完全
2.2灰色系统与模糊数学、黑箱方法的区别 主要在于对系统内涵与外延处理态度不同； 研究对象内涵与外延的性质不同。
灰色系统着重外延明确、内涵不明确的对象，模糊数学着重外延不明确、内涵明确的对象。
“黑箱”方法着重系统外部行为数据的处理方法，是因果关系的两户方法，使扬外延而弃内涵的处理方法，而灰色系统方法是外延内涵均注重的方法。
2.3灰色系统的基本原理 公理1：差异信息原理。“差异”是信息，凡信息必有差异。
公理2：解的非唯一性原理。信息不完全，不明确地解是非唯一的。 公理3：最少信息原理。灰色系统理论的特点是充分开发利用已有的“最少信息”。
公理4：认知根据原理。信息是认知的根据。
公理5：新信息优先原理。新信息对认知的作用大于老信息。 公理6：灰性不灭原理。“信息不完全”是绝对的。
2.4灰色系统理论的主要内容
灰色系统理论经过10多年的发展，已基本建立起了一门新兴学科的结构体系，其主要内容包括以“灰色朦胧集”为基础的理论体系、以晦涩关联空间为依托的分析体系、以晦涩序列生成为基础的方法体系，以灰色模型（G ，M ）为核心的模型体系。以系统分析、评估、建模、预测、决策、控制、优化为主体的技术体系。
灰色关联分析 灰色统计 灰色聚类
3灰色系统预测模型
灰色预测方法的特点表现在：首先是它把离散数据视为连续变量在其变化过程中所取的离散值，从而可利用微分方程式处理数据；而不直接使用原始数据而是由它产生累加生成数，对生成数列使用微分方程模型。这样，可以抵消大部分随机误差，显示出规律性。
3.1灰色系统理论的建模思想
下面举一个例子，说明灰色理论的建模思想。考虑4个数据，记为(
0) (1), (0) (2), (0) (3), (0) (4)
上图表明原始数据X 没有明显的规律性，其发展态势是摆动的。如果将原始数据作累加生成，记第K 个累加生成为X (1) (K ) , 并且
X (1) (1) =X (0) (1) =1
X (1) (2) =X (0) (1) +X (0) (2) =1+2=3
X (1) (3) =X (0) (1) +X (0) (2) +X (0) (3) =1+2+1. 5=4. 5
X (1) (4) =X (0) (1) +X (0) (2) +X (0) (3) +X (0) (4) =1+2+1. 5+3=7. 5
(0)
上图表明生成数列X 是单调递增数列。
3.2灰色系统预测模型建立 1. 数列预测GM （1，1）模型
灰色系统理论的微分方程成为Gm 模型，G 表示gray （灰色），m 表示model （模型），Gm （1，1）表示1阶的、1个变量的微分方程模型。
Gm （1，1）建模过程和机理如下： 记原始数据序列X (0) 为非负序列
X (0)=x (0)(1), x (0)(2), x (0)(3),..., x (0)(n )
{}
其中，x (0) (k ) ≥0, k =1, 2, , n 其相应的生成数据序列为X (1)
X (1)=x (1)(1), x (1)(2), x (1)(3),..,. x (1)(n )
(1)
k i =1
{}
其中，x (k ) =∑x (0) (i ) , k =1, 2, , n
Z (1) 为X (1) 的紧邻均值生成序列 Z (1) =z (1) (1), z (1) (2), , z (1) (n )
{}
其中，Z (1) (k ) =0. 5x (1) (k ) +0. 5x (1) (k -1), k =1, 2, n
称x (0) (k ) +az (1) (k ) =b 为Gm(1,1)模型，其中a ，b 是需要通过建模求解的参数，若a =(a , b ) T为参数列，且
⎡-z (1) (2) 1⎤⎡x (0) (2) ⎤
⎢-z (1) (3) 1⎥⎢x (0) (3) ⎥
Y =⎢⎥，B =⎢-z (1) (4) 1⎥ ⎢⎥⎢(0) ⎥(1)
⎢⎣x (n ) ⎦⎣-z (5) ⎥⎦
则求微分方程x (0) (k ) +az (1) (k ) =b 的最小二乘估计系数列，满足
ˆ=(B T B ) -1B T Y a dx (1)
+ax (1) =b 为灰微分方程，x (0) (k ) +az (1) (k ) =b 的白化方程，称dt
也叫影子方程。
如上所述，则有
dx (1)
+ax (1) =b 的解或称时间响应函数为 1. 白化方程dt
b b
ˆ(1) (t ) =(x (1) (0) -) e -at + x
a a
2.Gm(1,1)灰微分方程x (0) (k ) +az (1) (k ) =b 的时间响应序列为
b b
ˆ(1) (k +1) =(x (1) (0) -) e -ak +, k =1, 2, , n x
a a
3. 取x (1) (0) =x (0) (1) ，则
b b
ˆ(1) (k +1) =(x (0) (1) -) e -ak +, k =1, 2, , n x
a a
4. 还原值
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ), k =1, 2, , n x
2. 系统综合预测GM （1，N ）模型P134
4灰色系统模型的检验
定义1.
设原始序列
X (0) =x (0) (1), x (0) (2), , x (0) (n )
{}
相应的模型模拟序列为
ˆ(0) =x ˆ(0) (1), x ˆ(0) (2), , x ˆ(0) (n ) X
{}
残差序列
ˆ(0) (1), x (0) (2) -x ˆ(0) (2), , x (0) (n ) -x ˆ(0) (n ) } ={x (0) (1) -x 相对误差序列
⎧ε(1) ε(2) ε(n ) ⎫∆=⎨(0) , (0) , , (0) ⎬
x (1) x (2) x (n ) ⎩⎭
n
={∆k }1
ε(0) ={ε(1), ε(2), ε(n ) }
1. 对于k ＜n, 称∆k =
ε(k )
x (0) (k )
为k 点模拟相对误差，称∆n =
ε(n )
x (0) (n )
1n
为滤波相对误差，称=∑∆k 为平均模拟相对误差；
n k =1
2. 称1-为平均相对精度，1-∆n 为滤波精度；
3. 给定α，当
定义2
ˆ(0) 为相应的模拟误差序列，ˆ(0) X 设X (0) 为原始序列，ε为X (0) 与X
的绝对关联度，若对于给定的ε0>0, ε>ε0，则称模型为关联合格模型。
定义3
ˆ(0) 为相应的模拟误差序列，ε(0) 为残差序设X (0) 为原始序列，X
列。
1n (0)
=∑x (k ) 为X (0) 的均值， n k =1
1n (0) 2
s 1=∑(x (k ) -) 2为x (0) 的方差，
n k =11n
=∑ε(k ) 为残差均值，
n k =11n 2
s 2=∑(ε(k ) -) 2为残差方差，
n k =1s
1. 称c =2为均方差比值；对于给定的c 0>0，当c
s 1
为均方差比合格模型。
2. 称p =p ε(k ) -0, 当p >p 0时，称模型为小误差概率合格模型。
5应用举例
例 1 设原始序列
X (0) =x (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-AGO ，得
[D为X (0) 的一次累加生成算子，记为1-AGO ，A cumulated Generating Operator]
X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4), x (1) (5) }
=(2. 874, 6. 152, 9. 489, 12. 579, 16. 558)
2) 对X (1) 作紧邻均值生成，令
Z (1) (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 (1) (5)
=(2. 874, 4. 513, 7. 820, 11. 84, 14. 718)
{}
于是，
(0)
⎡x (2) ⎤1⎤⎡3. 278⎤(0) ⎢⎥1⎥，Y =⎢x (3) ⎥=⎢3. 337⎥ (0) 1⎥x (4) ⎥⎢3. 390⎥⎢1⎦⎥⎢3. 679⎦⎥(0) ⎣⎢⎥x (5) ⎣⎦
⎡-4. 5131⎤
. 513-7. 820-11. 184-14. 718⎤∙⎢-7. 8201⎥ B TB =⎡-41111⎥⎢⎣⎦⎢-11. 841⎥
⎢⎣-14. 7181⎥⎦
423. 221-38. 235⎤ =⎡4⎢⎥⎣-38. 235⎦
⎡-z (1) (2) ⎢-z (1) (3) B =⎢(1)
⎢-z (4)
(1) ⎢⎣-z (5) 1⎤
⎡-4. 513⎥1⎢-7. 820⎥=
1⎥⎢-11. 84
⎢-14. 718⎣1⎥⎦
423. 221-38. 235⎤==⎡0. 0173180. 165542⎤ (B B ) =⎡-4⎥⎢⎢⎥⎣38. 235⎦⎣0. 16655421. 832371⎦1⎧38. 235⎤⎫⎡4=⎪423. 221⨯4-38. 2352⎢⎥⎣38. 235423. 221⎦⎪
⎪⎪⎨⎬
138. 235⎤⎡4⎪=⎪
38. 235423. 221⎥⎪230. 969⎢⎪⎣⎦
⎩⎭
ˆ=(B TB ) -1B TY a
T
-1
-1
⎡3. 278⎤0. 0173180. 165542⎤∙⎡-4. 513-7. 820-11. 184-14. 718⎤∙⎢3. 337⎥=⎡0111⎥⎢⎥⎣. 16655421. 832371⎦⎢⎣1⎦⎢3. 390⎥
⎢⎣3. 679⎥⎦
⎡3. 278⎤
0. 0873860. 030115-0. 028143-0. 089344⎤∙⎢3. 337⎥ =⎡1⎢⎣. 0852800. 537833-0. 01905110. 604076⎥⎦⎢3. 390⎥
⎢⎣3. 679⎥⎦
. 037156⎤ =⎡-30
⎣. 065318⎥⎦ ⎢
3）确定模型
dx (1)
-0. 037156x (1) =3. 065318 dt
及时间响应式
b b
ˆ(1) (k +1) =(x (0) (1) -) e -ak + x
a a
e 0. 037156k -82. 4986 =85. 3728
4) 求X (1) 的模拟值
ˆ(1) =x ˆ(1) (1), (1) (2), x ˆ(1) (3), x ˆ(1) (4), x ˆ(1) (5) X x
{}
=(2.8740,6.1058,9.4599,12.9410,16.5538) 5) 还原出X (0) 的模拟值，由
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) x
ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4), x ˆ(0) (5) } 得 X
=（2.8740，3.2318，3.3541，3.4811，3.6128）
⎡ε(2) ⎤⎢ε(3) ⎥
s =εTε=[ε(2) ε(3) ε(4) ε(5) ]∙⎢
ε(4) ⎥⎢ε(5) ⎥⎣⎦
⎡0. 0462⎤
0. 0171⎥ =[0. 0462-0. 0171-0. 09110. 0662]∙⎢-
⎢-0. 0911⎥⎢⎣0. 0662⎥⎦
=0.0151085
平均相对误差
151
∆=∑∆k =(1. 41%+0. 51%+2. 69%+1. 80%)
4k =14
=1.0625%
ˆ的灰色关联度 计算X 与X
S =
1
(x (k ) -x (1) +(x (5) -x (1)) ∑2k =2
1
2
4
=(3. 278-2. 874) +(3. 337-2. 874) +(3. 390-2. 874) +(3. 679-2. 874)
=0.404+0.463+0.516+
=1.7855
ˆ=S
1
ˆˆˆ(5) -x ˆ(1) (x (k ) -x (1) +(x ∑2k =2
4
1
=(3. 2318-2. 874) +(3. 3541-2. 874) +(3. 4811-2. 874) +(3. 6128-2. 874)
2
=0. 3578+0. 4801+0. 6071+0.
=1.8144
ˆ-S =S
1
ˆˆˆ(5) -x ˆ(1)) [][(x (5) -x (1)) -(x (x (k ) -x (1)) -(x (k ) -x (1)) +∑
k =24
2
1
=(0. 3578-0. 404) +(0. 4801-0. 463) +(0. 6071-0. 516) +(0. 3694-0. 4025)
2
=-0. 0462+0. 0171+0. 091-0.
=0.04535
ε=
ˆ1+S +S ˆ+S ˆ-S 1+S +S
=
1+1. 7855+1. 81444. 5999
=
1+1. 7855+1. 8144+0. 045354. 64525
=0.9902＞0.90
精度为一级，可以用
ˆ(1) (k +1) =85. 3728x e 0. 037156k -82. 4986 ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) 预测。 x
例 2
试建立Gm(1,1)模型的白化方程及时间响应式，并对Gm(1,1)模型进行检验，预测该企业2001-2005年产值。 解：设时间序列为
=（27260，29547，62411，35388） X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4) }
) =(27260, 56807, 89218, 124606
X (0) =x (0) (1), x (0) (2), x (0) (3), x (0) (4)
{}
2) 对X (1) 作紧邻均值生成，令
Z (1) (k ) =0. 5x (1) (k ) +0. 5x (1) (k -1) Z (1) =z (1) (1), z (1) (2), z (1) (3), z (1) (4)
) =(27260, 42033. 5, 73012. 5, 106912
{}
于是，
⎡-z (1) (2) 1⎤⎡-42033⎡x (0) (2) ⎤⎡29547. 51⎤⎤(1) (0) ⎢⎥⎢⎥B =-z (3) 1=⎢-73012. 51，Y =x (3) =⎢32411⎥ ⎢(1) ⎥⎢1069121⎥⎢⎥(0) ⎥⎢35388⎥⎣⎦⎦-z (4) 1x ⎢⎥⎢⎣⎦⎣(4) ⎥⎦⎣
ˆ=[a , b ]T作最小二乘估计，得 对参数列α
ˆ=(B TB ) -1B TY =⎡-0. 089995⎤ a
⎢⎣25790. 28⎥⎦
dx (1)
-ax (1) =b 设dt
由于
可得Gm(1,1)模型的白化方程
dx (1)
--0. 089995x (1) =25790. 28dt
a =-0. 089995, b =25790. 28
其时间响应式为
b b ⎧(1)
ˆ(k +1) =(x (0) (1) -) e -ak +=313834e 0. 089995k -286574⎪x
⎨a a
(0) (1) (1)
⎪ˆ(k +1) =x ˆ(k +1) -x ˆ(k ) x ⎩
由此得模拟序列 ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4) } X
=（27260，29553，32336，35381） 检验：
残差序列为
ε(0) =(ε(0) (1), ε(0) (2), ε(0) (3), ε(0) (4))
=（0，-6，75，7）
⎧ε(0) (1) ε(0) (2) ε(0) (3) ε(0) (4) ⎫∆=⎨(0) , (0) , (0) , (0) ⎬ ⎩x (1) x (2) x (3) x (4) ⎭
=(0, 0. 0002, 0. 00231, 0. 0002) ∆=(∆1, ∆2, ∆3, ∆4)
平均相对误差
14
∆=∑∆k =0. 00068=0. 068%
模拟误差∆4=0. 0002=0. 02%
ˆ(0) 的灰色关联度ε 计算X (0) 与X
S =
ˆ=S
ˆ-S =S ˆ∑[(x
k =23(0) ∑(x k =233(0) 1(k ) -x (0) (1) +(x (0) (4) -x (0) (1)) =11502 21ˆˆˆ(4) -x ˆ(1) =11429.5 (x (k ) -x (1) +(x ∑2k =2ˆ(0) (1)) -(x (0) (k ) -x (0) (1)) +(k ) -x ˆ]1[(x 2(0) ˆ(0) (1)) -(x (0) (4) -x (0) (1)) (4) -x
=72.5
ε=ˆ1+S +S
ˆ+S ˆ-S 1+S +S =1+11502+11429. 5=0. 997>0. 90 1+11502+11429. 5+72. 5
精度为一级
计算均方差比
14(0) =∑x (k ) =31151. 5 4k =1
14(0) 2S 1=∑(x (k ) -) 2=9313116. 25, S 1=3051. 74 4k =1
14=∑ε(K ) =19 4K =1
142S 2=∑(ε(k ) -) 2=1066. 5, s 2=32. 66 4k =1
S 32. 66所以，c =2==0. 011
计算小误差概率
0. 6745S 1=2058. 40
ε(1) -=19, (2) -=25 ε(3) -=56ε(4) -=12
所以，P =P (ε(k ) -0. 95
小误差概率为一级，故可用
⎧x ˆ(1) (k +1) =313834e 0. 089995k -286574 ⎨ˆ(0) (1) (1) ˆˆx (k +1) =x (k +1) -x (k ) ⎩
进行预测，2001-2005年预测值为
ˆ(0) =x ˆ(0) (5), x ˆ(0) (6), x ˆ(0) (7), x ˆ(0) (8), x ˆ(0) (9) X {}
=（38713，42359，46318，50712，55488）
例3
解： X (0) ={x (0) (1), x (0) (2), x (0) (3), x (0) (4), x (0) (5) }
=(1. 67, 1. 51, 1. 03, 2. 14, 1. 99)
X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4), x (1) (5) }
=(1. 67, 3. 18, 4. 21, 6. 35, 8. 43)
对X (1) 作紧邻均值生成，令
Z (1) (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 (1) (5) }
=(1. 67, 2. 425, 3. 695, 5. 28, 7. 345)
于是，
(0) ⎡x (2) ⎤1⎤⎡1. 51⎤(0) ⎢⎥1⎥，Y =⎢x (3) ⎥=⎢1. 03⎥ (0) 1⎥x (4) ⎥⎢2. 14⎥⎢1⎥⎢(0) ⎦⎣1. 99⎥⎦⎢⎥x (5) ⎣⎦
⎡-2. 4251⎤. 425-3. 695-5. 28-7. 345⎤∙⎢-3. 6951⎥ B TB =⎡-21111⎥⎢⎣⎦⎢-5. 281⎥⎢⎣-7. 3451⎥⎦
101. 361-18. 745⎤ =⎡4⎥⎢⎣-18. 745⎦
0. 073980. 34669⎤ (B TB ) -1=⎡0⎢⎣. 346691. 8747⎥⎦
⎡1. 51⎤. 425-3. 695-5. 28-7. 345⎤∙⎢1. 03⎥ B TY =⎡-21111⎥⎢⎣⎦⎢2. 14⎥⎢⎣1. 99⎥⎦
. 38335⎤ =⎡-33⎢⎣6. 67⎥⎦
ˆ=⎡a ⎤=(B TB ) -1B TY =⎡0. 0379800. 34669⎤∙⎡-33. 38335⎤ a ⎢⎢⎣b ⎥⎦⎣0. 346691. 8747⎥⎦⎢⎣6. 67⎥⎦⎡-z (1) (2) ⎢-z (1) (3) B =⎢(1) ⎢-z (4) (1) ⎢⎣-z (5) 1⎤⎡-2. 425⎥1⎢-3. 695⎥=1⎥⎢-5. 28⎢⎣-7. 3451⎥⎦
. 157⎤ =⎡-00⎢⎣. 931⎥⎦
dx (1)
-ax (1) =b 方程为dt (1) dx -0. 157x (1) =0. 931 dt
及时间响应式
b b ˆ(1) (k +1) =(x (0) (1) -) e -ak + x a a
=(1. 67+5. 93) e 0. 157k -5. 93
=7. 6e 0. 157k -5. 93
求X (1) 的模拟值
ˆ(1) =x ˆ(1) (1), (1) (2), x ˆ(1) (3), x ˆ(1) (4), x ˆ(1) (5) X x {}
=(1.67,2.962,4.474,6.202,8.311)
还原出X (0) 的模拟值，由
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) x
ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4), x ˆ(0) (5) } 得 X
=（1.67，1.292，1.512，1.728，2.109）
ˆ的灰色关联度 计算X 与X
S =1(x (k ) -x (1) +(x (5) -x (1)) ∑2k =2
1=(-0. 61-0. 64+0. 47+⨯0. 32=0. 17 2
11ˆˆˆˆ=-0. 378-0. 1580. 058+⨯0. (x (k ) -x (1) +(x (5) -x (1) ∑22k =244ˆ=S
=0.2585
ˆ-S =S 1ˆˆˆ(5) -x ˆ(1)) [][(x (5) -x (1)) -(x (x (k ) -x (1)) -(x (k ) -x (1)) +∑k =242
=(-0. 378+0. 16+(-0. 158+0. 64) +(0. 058-0. 47+(0. 439-0. 32)
=-0. 218+0. 482-0. 417+0. 12
=0.0885
ε=ˆ1+S +S
ˆ+S ˆ-S 1+S +S =1+0. 17+0. 25851. 4285= 1+0. 17+0. 2585+0. 08851. 517
=0.94＞0.90
精度为一级，关联度为一级，可以用
⎧x ˆ(1) (k +1) =7. 6e 0. 157k -5. 93 ⎨ˆ(0) (1) (1) ˆˆ⎩x (k +1) =x (k +1) -x (k )
进行预测。
ˆ(1) =x ˆ(1) (6), x ˆ(1) (7), x ˆ(1) (8), x ˆ(1) (9), x ˆ(1) (10), X
ˆ(1) (11), x ˆ10) (12), x ˆ(1) (13), x ˆ(1) (14), x ˆ(1) (15) } x
=(10. 732, 13. 565, 16. 879, 20. 756, 25. 293, 30. 601, 36. 811, 44. 076, 52. 577, 62. 523) ˆ(0) =(x ˆ(0) (6), x ˆ(0) (7), x ˆ(0) (8), x ˆ(0) (9), x ˆ(0) (10), , x ˆ(0) (15) )X {
=(2. 42, 21. 83, 3. 31, 43. 87, 74. 53, 7 5. 308, 6. 21, 7. 265, 8. 501, 9. 946)
灰色系统预测
重点内容：灰色系统理论的产生和发展动态，灰色系统的基本概念，灰色系统与模糊数学、黑箱方法的区别，灰色系统预测GM （1，1）模型，GM （1，N ）模型，灰色系统模型的检验，应用举例。
1灰色系统理论的产生和发展动态
1982邓聚龙发表第一篇中文论文《灰色控制系统》标志着灰色系统这一学科诞生。
1985灰色系统研究会成立，灰色系统相关研究发展迅速。 1989海洋出版社出版英文版《灰色系统论文集》，同年，英文版国际刊物《灰色系统》杂志正式创刊。目前，国际、国内200多种期刊发表灰色系统论文，许多国际会议把灰色系统列为讨论专题。国际著名检索已检索我国学者的灰色系统论著500多次。灰色系统理论已应用范围已拓展到工业、农业、社会、经济、能源、地质、石油等众多科学领域，成功地解决了生产、生活和科学研究中的大量实际问题，取得了显著成果。
2灰色系统的基本原理
2.1灰色系统的基本概念
我们将信息完全明确的系统称为白色系统，信息未知的系统称为黑色系统，部分信息明确、部分信息不明确的系统称为灰色系统。 系统信息不完全的情况有以下四种： 1. 元素信息不完全 2. 结构信息不完全 3. 边界信息不完全 4. 运行行为信息不完全
2.2灰色系统与模糊数学、黑箱方法的区别 主要在于对系统内涵与外延处理态度不同； 研究对象内涵与外延的性质不同。
灰色系统着重外延明确、内涵不明确的对象，模糊数学着重外延不明确、内涵明确的对象。
“黑箱”方法着重系统外部行为数据的处理方法，是因果关系的两户方法，使扬外延而弃内涵的处理方法，而灰色系统方法是外延内涵均注重的方法。
2.3灰色系统的基本原理 公理1：差异信息原理。“差异”是信息，凡信息必有差异。
公理2：解的非唯一性原理。信息不完全，不明确地解是非唯一的。 公理3：最少信息原理。灰色系统理论的特点是充分开发利用已有的“最少信息”。
公理4：认知根据原理。信息是认知的根据。
公理5：新信息优先原理。新信息对认知的作用大于老信息。 公理6：灰性不灭原理。“信息不完全”是绝对的。
2.4灰色系统理论的主要内容
灰色系统理论经过10多年的发展，已基本建立起了一门新兴学科的结构体系，其主要内容包括以“灰色朦胧集”为基础的理论体系、以晦涩关联空间为依托的分析体系、以晦涩序列生成为基础的方法体系，以灰色模型（G ，M ）为核心的模型体系。以系统分析、评估、建模、预测、决策、控制、优化为主体的技术体系。
灰色关联分析 灰色统计 灰色聚类
3灰色系统预测模型
灰色预测方法的特点表现在：首先是它把离散数据视为连续变量在其变化过程中所取的离散值，从而可利用微分方程式处理数据；而不直接使用原始数据而是由它产生累加生成数，对生成数列使用微分方程模型。这样，可以抵消大部分随机误差，显示出规律性。
3.1灰色系统理论的建模思想
下面举一个例子，说明灰色理论的建模思想。考虑4个数据，记为(
0) (1), (0) (2), (0) (3), (0) (4)
上图表明原始数据X 没有明显的规律性，其发展态势是摆动的。如果将原始数据作累加生成，记第K 个累加生成为X (1) (K ) , 并且
X (1) (1) =X (0) (1) =1
X (1) (2) =X (0) (1) +X (0) (2) =1+2=3
X (1) (3) =X (0) (1) +X (0) (2) +X (0) (3) =1+2+1. 5=4. 5
X (1) (4) =X (0) (1) +X (0) (2) +X (0) (3) +X (0) (4) =1+2+1. 5+3=7. 5
(0)
上图表明生成数列X 是单调递增数列。
3.2灰色系统预测模型建立 1. 数列预测GM （1，1）模型
灰色系统理论的微分方程成为Gm 模型，G 表示gray （灰色），m 表示model （模型），Gm （1，1）表示1阶的、1个变量的微分方程模型。
Gm （1，1）建模过程和机理如下： 记原始数据序列X (0) 为非负序列
X (0)=x (0)(1), x (0)(2), x (0)(3),..., x (0)(n )
{}
其中，x (0) (k ) ≥0, k =1, 2, , n 其相应的生成数据序列为X (1)
X (1)=x (1)(1), x (1)(2), x (1)(3),..,. x (1)(n )
(1)
k i =1
{}
其中，x (k ) =∑x (0) (i ) , k =1, 2, , n
Z (1) 为X (1) 的紧邻均值生成序列 Z (1) =z (1) (1), z (1) (2), , z (1) (n )
{}
其中，Z (1) (k ) =0. 5x (1) (k ) +0. 5x (1) (k -1), k =1, 2, n
称x (0) (k ) +az (1) (k ) =b 为Gm(1,1)模型，其中a ，b 是需要通过建模求解的参数，若a =(a , b ) T为参数列，且
⎡-z (1) (2) 1⎤⎡x (0) (2) ⎤
⎢-z (1) (3) 1⎥⎢x (0) (3) ⎥
Y =⎢⎥，B =⎢-z (1) (4) 1⎥ ⎢⎥⎢(0) ⎥(1)
⎢⎣x (n ) ⎦⎣-z (5) ⎥⎦
则求微分方程x (0) (k ) +az (1) (k ) =b 的最小二乘估计系数列，满足
ˆ=(B T B ) -1B T Y a dx (1)
+ax (1) =b 为灰微分方程，x (0) (k ) +az (1) (k ) =b 的白化方程，称dt
也叫影子方程。
如上所述，则有
dx (1)
+ax (1) =b 的解或称时间响应函数为 1. 白化方程dt
b b
ˆ(1) (t ) =(x (1) (0) -) e -at + x
a a
2.Gm(1,1)灰微分方程x (0) (k ) +az (1) (k ) =b 的时间响应序列为
b b
ˆ(1) (k +1) =(x (1) (0) -) e -ak +, k =1, 2, , n x
a a
3. 取x (1) (0) =x (0) (1) ，则
b b
ˆ(1) (k +1) =(x (0) (1) -) e -ak +, k =1, 2, , n x
a a
4. 还原值
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ), k =1, 2, , n x
2. 系统综合预测GM （1，N ）模型P134
4灰色系统模型的检验
定义1.
设原始序列
X (0) =x (0) (1), x (0) (2), , x (0) (n )
{}
相应的模型模拟序列为
ˆ(0) =x ˆ(0) (1), x ˆ(0) (2), , x ˆ(0) (n ) X
{}
残差序列
ˆ(0) (1), x (0) (2) -x ˆ(0) (2), , x (0) (n ) -x ˆ(0) (n ) } ={x (0) (1) -x 相对误差序列
⎧ε(1) ε(2) ε(n ) ⎫∆=⎨(0) , (0) , , (0) ⎬
x (1) x (2) x (n ) ⎩⎭
n
={∆k }1
ε(0) ={ε(1), ε(2), ε(n ) }
1. 对于k ＜n, 称∆k =
ε(k )
x (0) (k )
为k 点模拟相对误差，称∆n =
ε(n )
x (0) (n )
1n
为滤波相对误差，称=∑∆k 为平均模拟相对误差；
n k =1
2. 称1-为平均相对精度，1-∆n 为滤波精度；
3. 给定α，当
定义2
ˆ(0) 为相应的模拟误差序列，ˆ(0) X 设X (0) 为原始序列，ε为X (0) 与X
的绝对关联度，若对于给定的ε0>0, ε>ε0，则称模型为关联合格模型。
定义3
ˆ(0) 为相应的模拟误差序列，ε(0) 为残差序设X (0) 为原始序列，X
列。
1n (0)
=∑x (k ) 为X (0) 的均值， n k =1
1n (0) 2
s 1=∑(x (k ) -) 2为x (0) 的方差，
n k =11n
=∑ε(k ) 为残差均值，
n k =11n 2
s 2=∑(ε(k ) -) 2为残差方差，
n k =1s
1. 称c =2为均方差比值；对于给定的c 0>0，当c
s 1
为均方差比合格模型。
2. 称p =p ε(k ) -0, 当p >p 0时，称模型为小误差概率合格模型。
5应用举例
例 1 设原始序列
X (0) =x (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-AGO ，得
[D为X (0) 的一次累加生成算子，记为1-AGO ，A cumulated Generating Operator]
X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4), x (1) (5) }
=(2. 874, 6. 152, 9. 489, 12. 579, 16. 558)
2) 对X (1) 作紧邻均值生成，令
Z (1) (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 (1) (5)
=(2. 874, 4. 513, 7. 820, 11. 84, 14. 718)
{}
于是，
(0)
⎡x (2) ⎤1⎤⎡3. 278⎤(0) ⎢⎥1⎥，Y =⎢x (3) ⎥=⎢3. 337⎥ (0) 1⎥x (4) ⎥⎢3. 390⎥⎢1⎦⎥⎢3. 679⎦⎥(0) ⎣⎢⎥x (5) ⎣⎦
⎡-4. 5131⎤
. 513-7. 820-11. 184-14. 718⎤∙⎢-7. 8201⎥ B TB =⎡-41111⎥⎢⎣⎦⎢-11. 841⎥
⎢⎣-14. 7181⎥⎦
423. 221-38. 235⎤ =⎡4⎢⎥⎣-38. 235⎦
⎡-z (1) (2) ⎢-z (1) (3) B =⎢(1)
⎢-z (4)
(1) ⎢⎣-z (5) 1⎤
⎡-4. 513⎥1⎢-7. 820⎥=
1⎥⎢-11. 84
⎢-14. 718⎣1⎥⎦
423. 221-38. 235⎤==⎡0. 0173180. 165542⎤ (B B ) =⎡-4⎥⎢⎢⎥⎣38. 235⎦⎣0. 16655421. 832371⎦1⎧38. 235⎤⎫⎡4=⎪423. 221⨯4-38. 2352⎢⎥⎣38. 235423. 221⎦⎪
⎪⎪⎨⎬
138. 235⎤⎡4⎪=⎪
38. 235423. 221⎥⎪230. 969⎢⎪⎣⎦
⎩⎭
ˆ=(B TB ) -1B TY a
T
-1
-1
⎡3. 278⎤0. 0173180. 165542⎤∙⎡-4. 513-7. 820-11. 184-14. 718⎤∙⎢3. 337⎥=⎡0111⎥⎢⎥⎣. 16655421. 832371⎦⎢⎣1⎦⎢3. 390⎥
⎢⎣3. 679⎥⎦
⎡3. 278⎤
0. 0873860. 030115-0. 028143-0. 089344⎤∙⎢3. 337⎥ =⎡1⎢⎣. 0852800. 537833-0. 01905110. 604076⎥⎦⎢3. 390⎥
⎢⎣3. 679⎥⎦
. 037156⎤ =⎡-30
⎣. 065318⎥⎦ ⎢
3）确定模型
dx (1)
-0. 037156x (1) =3. 065318 dt
及时间响应式
b b
ˆ(1) (k +1) =(x (0) (1) -) e -ak + x
a a
e 0. 037156k -82. 4986 =85. 3728
4) 求X (1) 的模拟值
ˆ(1) =x ˆ(1) (1), (1) (2), x ˆ(1) (3), x ˆ(1) (4), x ˆ(1) (5) X x
{}
=(2.8740,6.1058,9.4599,12.9410,16.5538) 5) 还原出X (0) 的模拟值，由
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) x
ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4), x ˆ(0) (5) } 得 X
=（2.8740，3.2318，3.3541，3.4811，3.6128）
⎡ε(2) ⎤⎢ε(3) ⎥
s =εTε=[ε(2) ε(3) ε(4) ε(5) ]∙⎢
ε(4) ⎥⎢ε(5) ⎥⎣⎦
⎡0. 0462⎤
0. 0171⎥ =[0. 0462-0. 0171-0. 09110. 0662]∙⎢-
⎢-0. 0911⎥⎢⎣0. 0662⎥⎦
=0.0151085
平均相对误差
151
∆=∑∆k =(1. 41%+0. 51%+2. 69%+1. 80%)
4k =14
=1.0625%
ˆ的灰色关联度 计算X 与X
S =
1
(x (k ) -x (1) +(x (5) -x (1)) ∑2k =2
1
2
4
=(3. 278-2. 874) +(3. 337-2. 874) +(3. 390-2. 874) +(3. 679-2. 874)
=0.404+0.463+0.516+
=1.7855
ˆ=S
1
ˆˆˆ(5) -x ˆ(1) (x (k ) -x (1) +(x ∑2k =2
4
1
=(3. 2318-2. 874) +(3. 3541-2. 874) +(3. 4811-2. 874) +(3. 6128-2. 874)
2
=0. 3578+0. 4801+0. 6071+0.
=1.8144
ˆ-S =S
1
ˆˆˆ(5) -x ˆ(1)) [][(x (5) -x (1)) -(x (x (k ) -x (1)) -(x (k ) -x (1)) +∑
k =24
2
1
=(0. 3578-0. 404) +(0. 4801-0. 463) +(0. 6071-0. 516) +(0. 3694-0. 4025)
2
=-0. 0462+0. 0171+0. 091-0.
=0.04535
ε=
ˆ1+S +S ˆ+S ˆ-S 1+S +S
=
1+1. 7855+1. 81444. 5999
=
1+1. 7855+1. 8144+0. 045354. 64525
=0.9902＞0.90
精度为一级，可以用
ˆ(1) (k +1) =85. 3728x e 0. 037156k -82. 4986 ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) 预测。 x
例 2
试建立Gm(1,1)模型的白化方程及时间响应式，并对Gm(1,1)模型进行检验，预测该企业2001-2005年产值。 解：设时间序列为
=（27260，29547，62411，35388） X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4) }
) =(27260, 56807, 89218, 124606
X (0) =x (0) (1), x (0) (2), x (0) (3), x (0) (4)
{}
2) 对X (1) 作紧邻均值生成，令
Z (1) (k ) =0. 5x (1) (k ) +0. 5x (1) (k -1) Z (1) =z (1) (1), z (1) (2), z (1) (3), z (1) (4)
) =(27260, 42033. 5, 73012. 5, 106912
{}
于是，
⎡-z (1) (2) 1⎤⎡-42033⎡x (0) (2) ⎤⎡29547. 51⎤⎤(1) (0) ⎢⎥⎢⎥B =-z (3) 1=⎢-73012. 51，Y =x (3) =⎢32411⎥ ⎢(1) ⎥⎢1069121⎥⎢⎥(0) ⎥⎢35388⎥⎣⎦⎦-z (4) 1x ⎢⎥⎢⎣⎦⎣(4) ⎥⎦⎣
ˆ=[a , b ]T作最小二乘估计，得 对参数列α
ˆ=(B TB ) -1B TY =⎡-0. 089995⎤ a
⎢⎣25790. 28⎥⎦
dx (1)
-ax (1) =b 设dt
由于
可得Gm(1,1)模型的白化方程
dx (1)
--0. 089995x (1) =25790. 28dt
a =-0. 089995, b =25790. 28
其时间响应式为
b b ⎧(1)
ˆ(k +1) =(x (0) (1) -) e -ak +=313834e 0. 089995k -286574⎪x
⎨a a
(0) (1) (1)
⎪ˆ(k +1) =x ˆ(k +1) -x ˆ(k ) x ⎩
由此得模拟序列 ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4) } X
=（27260，29553，32336，35381） 检验：
残差序列为
ε(0) =(ε(0) (1), ε(0) (2), ε(0) (3), ε(0) (4))
=（0，-6，75，7）
⎧ε(0) (1) ε(0) (2) ε(0) (3) ε(0) (4) ⎫∆=⎨(0) , (0) , (0) , (0) ⎬ ⎩x (1) x (2) x (3) x (4) ⎭
=(0, 0. 0002, 0. 00231, 0. 0002) ∆=(∆1, ∆2, ∆3, ∆4)
平均相对误差
14
∆=∑∆k =0. 00068=0. 068%
模拟误差∆4=0. 0002=0. 02%
ˆ(0) 的灰色关联度ε 计算X (0) 与X
S =
ˆ=S
ˆ-S =S ˆ∑[(x
k =23(0) ∑(x k =233(0) 1(k ) -x (0) (1) +(x (0) (4) -x (0) (1)) =11502 21ˆˆˆ(4) -x ˆ(1) =11429.5 (x (k ) -x (1) +(x ∑2k =2ˆ(0) (1)) -(x (0) (k ) -x (0) (1)) +(k ) -x ˆ]1[(x 2(0) ˆ(0) (1)) -(x (0) (4) -x (0) (1)) (4) -x
=72.5
ε=ˆ1+S +S
ˆ+S ˆ-S 1+S +S =1+11502+11429. 5=0. 997>0. 90 1+11502+11429. 5+72. 5
精度为一级
计算均方差比
14(0) =∑x (k ) =31151. 5 4k =1
14(0) 2S 1=∑(x (k ) -) 2=9313116. 25, S 1=3051. 74 4k =1
14=∑ε(K ) =19 4K =1
142S 2=∑(ε(k ) -) 2=1066. 5, s 2=32. 66 4k =1
S 32. 66所以，c =2==0. 011
计算小误差概率
0. 6745S 1=2058. 40
ε(1) -=19, (2) -=25 ε(3) -=56ε(4) -=12
所以，P =P (ε(k ) -0. 95
小误差概率为一级，故可用
⎧x ˆ(1) (k +1) =313834e 0. 089995k -286574 ⎨ˆ(0) (1) (1) ˆˆx (k +1) =x (k +1) -x (k ) ⎩
进行预测，2001-2005年预测值为
ˆ(0) =x ˆ(0) (5), x ˆ(0) (6), x ˆ(0) (7), x ˆ(0) (8), x ˆ(0) (9) X {}
=（38713，42359，46318，50712，55488）
例3
解： X (0) ={x (0) (1), x (0) (2), x (0) (3), x (0) (4), x (0) (5) }
=(1. 67, 1. 51, 1. 03, 2. 14, 1. 99)
X (1) ={x (1) (1), x (1) (2), x (1) (3), x (1) (4), x (1) (5) }
=(1. 67, 3. 18, 4. 21, 6. 35, 8. 43)
对X (1) 作紧邻均值生成，令
Z (1) (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 (1) (5) }
=(1. 67, 2. 425, 3. 695, 5. 28, 7. 345)
于是，
(0) ⎡x (2) ⎤1⎤⎡1. 51⎤(0) ⎢⎥1⎥，Y =⎢x (3) ⎥=⎢1. 03⎥ (0) 1⎥x (4) ⎥⎢2. 14⎥⎢1⎥⎢(0) ⎦⎣1. 99⎥⎦⎢⎥x (5) ⎣⎦
⎡-2. 4251⎤. 425-3. 695-5. 28-7. 345⎤∙⎢-3. 6951⎥ B TB =⎡-21111⎥⎢⎣⎦⎢-5. 281⎥⎢⎣-7. 3451⎥⎦
101. 361-18. 745⎤ =⎡4⎥⎢⎣-18. 745⎦
0. 073980. 34669⎤ (B TB ) -1=⎡0⎢⎣. 346691. 8747⎥⎦
⎡1. 51⎤. 425-3. 695-5. 28-7. 345⎤∙⎢1. 03⎥ B TY =⎡-21111⎥⎢⎣⎦⎢2. 14⎥⎢⎣1. 99⎥⎦
. 38335⎤ =⎡-33⎢⎣6. 67⎥⎦
ˆ=⎡a ⎤=(B TB ) -1B TY =⎡0. 0379800. 34669⎤∙⎡-33. 38335⎤ a ⎢⎢⎣b ⎥⎦⎣0. 346691. 8747⎥⎦⎢⎣6. 67⎥⎦⎡-z (1) (2) ⎢-z (1) (3) B =⎢(1) ⎢-z (4) (1) ⎢⎣-z (5) 1⎤⎡-2. 425⎥1⎢-3. 695⎥=1⎥⎢-5. 28⎢⎣-7. 3451⎥⎦
. 157⎤ =⎡-00⎢⎣. 931⎥⎦
dx (1)
-ax (1) =b 方程为dt (1) dx -0. 157x (1) =0. 931 dt
及时间响应式
b b ˆ(1) (k +1) =(x (0) (1) -) e -ak + x a a
=(1. 67+5. 93) e 0. 157k -5. 93
=7. 6e 0. 157k -5. 93
求X (1) 的模拟值
ˆ(1) =x ˆ(1) (1), (1) (2), x ˆ(1) (3), x ˆ(1) (4), x ˆ(1) (5) X x {}
=(1.67,2.962,4.474,6.202,8.311)
还原出X (0) 的模拟值，由
ˆ(0) (k +1) =x ˆ(1) (k +1) -x ˆ(1) (k ) x
ˆ(0) ={x ˆ(0) (1), x ˆ(0) (2), x ˆ(0) (3), x ˆ(0) (4), x ˆ(0) (5) } 得 X
=（1.67，1.292，1.512，1.728，2.109）
ˆ的灰色关联度 计算X 与X
S =1(x (k ) -x (1) +(x (5) -x (1)) ∑2k =2
1=(-0. 61-0. 64+0. 47+⨯0. 32=0. 17 2
11ˆˆˆˆ=-0. 378-0. 1580. 058+⨯0. (x (k ) -x (1) +(x (5) -x (1) ∑22k =244ˆ=S
=0.2585
ˆ-S =S 1ˆˆˆ(5) -x ˆ(1)) [][(x (5) -x (1)) -(x (x (k ) -x (1)) -(x (k ) -x (1)) +∑k =242
=(-0. 378+0. 16+(-0. 158+0. 64) +(0. 058-0. 47+(0. 439-0. 32)
=-0. 218+0. 482-0. 417+0. 12
=0.0885
ε=ˆ1+S +S
ˆ+S ˆ-S 1+S +S =1+0. 17+0. 25851. 4285= 1+0. 17+0. 2585+0. 08851. 517
=0.94＞0.90
精度为一级，关联度为一级，可以用
⎧x ˆ(1) (k +1) =7. 6e 0. 157k -5. 93 ⎨ˆ(0) (1) (1) ˆˆ⎩x (k +1) =x (k +1) -x (k )
进行预测。
ˆ(1) =x ˆ(1) (6), x ˆ(1) (7), x ˆ(1) (8), x ˆ(1) (9), x ˆ(1) (10), X
ˆ(1) (11), x ˆ10) (12), x ˆ(1) (13), x ˆ(1) (14), x ˆ(1) (15) } x
=(10. 732, 13. 565, 16. 879, 20. 756, 25. 293, 30. 601, 36. 811, 44. 076, 52. 577, 62. 523) ˆ(0) =(x ˆ(0) (6), x ˆ(0) (7), x ˆ(0) (8), x ˆ(0) (9), x ˆ(0) (10), , x ˆ(0) (15) )X {
=(2. 42, 21. 83, 3. 31, 43. 87, 74. 53, 7 5. 308, 6. 21, 7. 265, 8. 501, 9. 946)