灰色模型介绍及应用

灰色理论基本知识
概言
有关名词概念
建模机理
灰色理论模型应用
（1，1）模型的应用——污染物浓度问题
GM（1，1）残差模型的应用——油菜发病率问题 GM模型在复杂问题中的应用——SARS 疫情问题 GM（1，n）模型的应用——因素相关问题
本章小结
思考题
推荐阅读书目
第十章灰色模型介绍及应用
灰色理论基本知识
概言
客观世界的很多实际问题，其内部的结构、参数以及特征并未全部被人们了解，人们不可能象研究白箱问题那样将其内部机理研究清楚，只能依据某种思维逻辑与推断来构造模型。

对这类部分信息已知而部分信息未知的系统，我们称之为灰色系统。

本章介绍的方法是从灰色系统的本征灰色出发，研究在信息大量缺乏或紊乱的情况下，如何对实际问题进行分析和解决。

灰色系统的研究对象是“部分信息已知、部分信息未知”的“小样本”、“贫信息”不确定性系统，它通过对“部分”已知信息的生成、开发实现对现实世界的确切描述和认识。

信息不完全是“灰”的基本含义。

灰色系统理论建模的主要任务是根据具体灰色系统的行为特征数据，充分开发并利用不多的数据中的显信息和隐信息，寻找因素间或因素本身的数学关系。

通常的办法是采用离散模型，建立一个按时间作逐段分析的模型。

但是，离散模型只能对客观系统的发展做短期分析，适应不了从现在起做较长远的分析、规划、决策的要求。

尽管连续系统的离散近似模型对许多工程应用来讲是有用的，但在某些研究领域中，人们却常常希望使用微分方程模型。

事实上，微分方程的系统描述了我们所希望辨识的系统内部的物理或化学过程的本质。

目前，灰色系统理论已成功地应用于工程控制、经济管理、未来学研究、生态系统及复杂多变的农业系统中，并取得了可喜的成就。

灰色系统理论有可能对社会、经济等抽象系统进行分析、建模、预测、决策和控制，它有可能成为人们认识客观系统改造客观系统的一个新型的理论工具。

有关名词概念
灰数：一个信息不完全的数，称为灰数。

灰元：信息不完全或内容难以穷尽的元素，称为灰元。

灰关系：信息不完全或机制不明确的关系，称为灰关系。

具有灰关系的因素是灰因素，灰因素之间的量化作用，称为灰关联。

灰色系统:含灰数、灰元或灰关系的系统称为信息不完全系统。

如果按照灰色理论去研究它。

则称此系统为灰色系统。

累加生成：由于灰系统对一切随机量都可看作是在一定范围内变化的灰色量，因此，为适应灰系统建模需要，提出“生成”的概念，“生成”即指对原始数据做累加（或累减）处理。

累加生成一般可写成AGO 。

若计(0)
x
为原始数列，()
r x
为r 次累加生成后数列，即
(0)(0)(0)(0){(1),(2),()}x x x x n = ()()()(){(1),(2),
()}r r r r x x x x n =
则r 次累加生成算式为
()(1)
(1)
(1)
(1)1
(1)(1)(1)(1)()(1)()(1)(2)()()
[(1)(2)(1)]()(1)()
k
r r r r r i r r r r r r x k x
x
x
k x i x x x k x k x k x k ----=-----=++==++
-+=-+∑ 一般常用的是一次累加生成，即
(1)
(0)1
(1)(0)()()
(1)()
k
i x k x i x k x k ===-+∑
建模机理
建立GM 模型，实际就是将原始数列经过累加生成后，建立具有微分、差分近似指数规律兼容的方程，成为灰色建模，所建模型称为灰色模型，简记为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 ）模型的特例，其简单的微分方程形式（白化形式的微分方程）是
+=dx
ax u dt
利用常数变易法解得，通解为
()-=+
at u x t ce a
若初始条件为00,()==t x t x ，则可得到微分方程的特解为
0()()-=-+at u u
x t x e a a
或时间响应函数
(1)(1)((1))-+=-+at u u
x t x e a a
其中白化微分方程中的ax 项中的x 为dx
dt
的背景值，也称为初始值; ,a u 为常数（有
时也将u 写成b ）。

按白化导数定义有差分形式的微分方程，即
()()
lim ∆→+∆-=∆t dx x t t x t dt t 显然，当时间密化值定义为1，即当1∆→t 时，上式可记为
1
[(1)()]lim ∆→=+-t dx
x t x t dt 记为离散形式
(1)()=+-dx
x t x t dt
这显然表明
dx
dt
是一次累计生成，因此上述方程可改写为 (1)(1)(0)(1)()(1)=+-=+dx
x t x t x t dt
这实际也表明，模型是以生成数(1)
x
（(1)x 是以(0)x 的一次累加）为基础的。

当∆t 足够小时，()x t 到()+∆x t t 不会发生突变，因此可取()x t 与()+∆x t t 的平均值作为0∆→t 时的背景值，因此，背景值便可记为
(1)(1)
(1)1[(1)()]2
=
++x x t x t 或
(1)(1)
(1)1[(1)()]2
=
++x x k x k 于是白化的微分方程(1)
(1)+=dx ax u dt
可改写为 (0)(1)(1)1
(1)[(1)()]2
++
++=x k a x k x k u 或
(0)(1)(1)1
(1)[(1)()]2
+=-
+++x k a x k x k u 即
(0)(1)(1)(0)(1)(1)(0)(1)(1)1
(2)[(2)(1)]21
(3)[(2)(1)]21
()[()(1)]2
=-++=-++=-
+-
+x a x x u x a x x u
x n a x n x n u
因此，上述方程可以改写为矩阵方程形式，即
(1)(1)(0)
(1)(1)(0)(0)(1)(1)
1[(2)(1)]21(2)1[(2)(1)]1(3)21()1[()(1)]2⎡⎤-+⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥-+-⎢⎥⎣
⎦
a x x x a x x x a u x n a x n x n 引入下列符号，设
(0)
(0)(0)(2)(3)()⎡⎤⎢⎥⎢
⎥=⎢⎥⎢⎥⎣⎦N x x Y x n 111⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦E (1)(1)
(1)(1)(1)(1)
1[(2)(1)]21[(2)(1)]21[()(1)]2⎡⎤-+⎢⎥⎢⎥
⎢⎥-+⎢⎥=⎢⎥⎢⎥
⎢⎥-+-⎢⎥⎣
⎦
a x x a x x X a x n x n 于是便有
[]⎡⎤
=+=⎢⎥⎣⎦
N a Y aX uE X E u
令
⎡⎤
=⎢⎥⎣⎦
a a u (1)(1)
(1)(1)(1)(1)1[(2)(1)]121[(2)(1)]1[]2
1[()(1)]12⎡⎤-+⎢⎥⎢⎥
⎢⎥-+⎢⎥==⎢⎥⎢⎥⎢⎥-+-⎢⎥⎣⎦
a x x a x x B X E a x n x n
则
[]⎡⎤
=+==⎢⎥⎣⎦
N a Y aX uE X E Ba u
解得
1()-⎡⎤
==⎢⎥⎣⎦
T T N a a B B B Y u
将求解得到的代入微分方程的解式（也称时间响应函数），则
(1)(1)(1)((1))-+=-+ak u u
x k x e a a
由于(0)
(1)(1)(1)=x
x ，因此求导还原得
(0)(0)(1)((1))-+=--ak u
x k a x e a
上述两式便为GM （1，1）的时间响应式，及灰色系统预测模型的基本算式，当然上述两式计算结果只是近似计算值。

为简记，一般可以将GM （1，1）的建模过程记为
(0)0(1)(0)((1);,)(1)(1)⋅⋅⇒⇒+⇒+IAGO GM AGOx GM x a u x k x k
灰色理论模型应用
（1，1）模型的应用——污染物浓度问题
GM （1，1）模型是灰色系统最基本的模型，下面以污染物浓度问题说明GM （1，1）模型的建立及求解过程。

例 某污染源中某种污染物质量浓度测量值如表，试建立GM （1，1）模型
表 某污染物质量浓度测量值 （mg/L ）
解：第一步，设原始数据为
(0)(0)(0)(0)((1),(2),,(6))(3.936,4.575,4.968,5.063,5.968,5.507)
==x x x x 第二步，对原始数据进行累加生成，即(1)(0)=x AGOx
(1)(0)(1)(1) 3.936==x x
(1)(1)(0)(2)[(1)(2)] 3.936 4.5758.511=+=+=x x x
(1)(1)(0)(3)[(2)(3)]13.479=+=x x x (1)(1)(0)(4)[(3)(4)]18.542=+=x x x (1)(1)(0)(5)[(4)(5)]24.510=+=x x x (1)(1)(0)(6)[(5)(6)]30.017=+=x x x
因此累加生成数据为
(1)(0)(3.936,8.511,13.479,18.542,24.510,30.017)==x AGOx
第三步，构造矩阵,N B Y
(1)(1)(1)(1)(1)(1)(1)
(1)(1)(1)
1[(1)(2)]121 -6.2235 1.0000[(2)(3)]12 -10.9950 1.00001
-16.0105 1.0000[(3)(4)]
12 -21[(4)(5)]121[(5)(6)]12⎡⎤-+⎢⎥⎢
⎥
⎢⎥-+⎢⎥⎢⎥⎢⎥==-+⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥-+⎢⎥⎣⎦
x x x x B x x x x x x 1.5260 1.0000 -27.2635 1.0000
⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎣
⎦ (0)(0)(0)[(2),(3),
,(6)][4.575 4.968 5.063 5.968 5.507]==T T
N Y x x x 第四步，计算1ˆ()-=T T N a
B B B Y 。

先求1
()
-T
B B ，即
1622.6 -82.0 -82.0 5⎡⎤
=⎢⎥
⎣⎦
T B B 根据逆矩阵的求解方法，得
1 0.0036 0.0592() 0.059
2 1.1706-⎡⎤
=⎢⎥
⎣⎦
T B B 再求T
N B Y 的值，即
-442.7641 26.0810⎡⎤
=⎢⎥
⎣⎦
T N B Y 进而求得ˆa
的值为 1
0.0036 0.0592-442.7641-0.0539a ˆ() 0.0592 1.1706 26.0810 4.3322u -⎡⎤⎡⎤⎡⎤⎡⎤====⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦
T
T
N a B B B Y 计算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); end B=.*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) B
b1=B.'*B
b2=inv(B.'*B) b3=B.'*YN.' b4=u/a
b5=X1(1)-b4 b6=-a*b5
第五步，将,a u 的值代入微分方程的时间响应函数，
令(1)
(1)ˆ(1)(0) 3.936==x
x ，得
(1)(1)0.0539ˆ(1)((1))84.326480.3904-+=-+=-ak k u u
x
k x e e a a
第六步，求导还原得
(0)(1)0.0539ˆ(1)((1)) 4.5443-+=--=ak k u
x
k a x e e a
第七步，对上述模型进行精度检验。

常用的方法是回代检验，即分别用(1)
(1)ˆ(1),(0)x
x 模型求出各时刻值，然后求相对误差。

先利用时间响应函数模型(1)
0.0539ˆ(1)84.326480.3904+=-k x
k e 求各时刻值
（1,2,
,5=k ）
，并计算相对误差，结果如表所示. 表 精度检验实测值、残差值表 1,2,
,5=k
再利用时间响应函数模型(0)
0.0539ˆ(1) 4.5443+=k x
k e 求各时刻值（1,2,
,5=k ）
，并计算相对误差，结果如表所示.
表 计算值与实验原始数据值对照表 1,2,
,5=k
从残差检验结果看，累计生成数列曲线拟合较好，相对误差在即1%左右；而还原数列的相对误差较大，其原因是累加生成数据将原始数据的随机性弱化，正负误差有抵消的，当数据再被还原回来时便表现出来。

GM （1，1）残差模型的应用——油菜发病率问题
当ＧＭ（１，１）模型的精度不符合要求时，可用残差序列建立ＧＭ（１，１）模型，对原来的模型进行修正，以提高精度，即建立残差ＧＭ（１，１）模型，步骤如下
第一步，利用原始数据建立GM （1，1）模型，得时间响应式
(1)(1)(1)((1))-+=-+ak u u
x k x e a a
(0)(0)(1)((1))-+=--ak u
x k a x e a
其中第二个式子也成为导数还原值。

鉴于导数还原值与原始数据（累减还原值）不一致，为减少往返运算造成的误差，往往用原始数据与导数还原值的残差修正的(0)
x
模拟值(0)
x。

第二步，利用残差数列建立新的GM （1，1）模型。

建立残差模型的过程和计算方法同于GM （1，1）建模过程，只不过建立残差模型所用的原始数列采用的是残差数据。

令(0)
()g
k 为残差，则
(0)(0)(0)()()()=-g k x k x k
即
(0)(0)(0)(0)((),(1),
())(,1,
,)=+=+g g k g k g n k i i n
或
(0)(0)(0)(0)((1),(2),
(1))=-+g g g g n i
利用残差序列(0)
g
建立新的GM （1，1）模型，求解得时间响应式
(1)(1)'0''(1)(())''
-+=-
+a k u u g k g k e a a (0)(1)'0'
(1)'(())'
-+=--a k u g k a g k e a
第三步，结合上两步的GM （1，1）模型，建立残差GM （1，1）模型 结合上两步的GM （1，1）模型，则相应的残差修正时间响应式为
()()(0)
(0)
(0)(1)'0
(1)(1)'(1)(')(())'---⎧--<⎪⎪+=⎨
⎪--+--≥⎪⎩
ak ak a k
u a x e k k a x k u u a x e a g k e k k a a
称为导数还原式的残差修正模型。

例 某县油菜发病率数据如表所示，试建立残差GM 模型并进行求解。

表 某县油菜发病率数据 （%）
解：第一步，建立原始数据的GM （1，1）模型
设原始数据为
()
(0)(0)(0)(0)((1),(2),,(13))
0.01*6,20,40,25,40,45,35,21,14,18,15.5,17,15==x x x x
建立GM （1，1）模型，利用GM （1，1）的求解程序得时间响应式为
()(1)
0.064861 5.680 5.740-+=-+k x k e ()(0)
0.0648610.368-+=k x
k e
第二步，误差检验 利用时间响应函数模型()(0)
0.0648610.368-+=k x k e 计算各时刻值（1,2,
,12=k ）
，并计算相对误差，程序如下
function 10toliti02(X0) %format long ;
%X0=*[6 20 40 25 40 45 35 21 14 18 17 15]; [m,n]=size(X0); s(1)=1; for i=1:12 y(i+1)=*exp*i); z(i+1)=X0(i+1)-y(i+1); w(i+1)=z(i+1)/X0(i+1); s(i+1)=i+1; end y' X0' z' w' z*z'
sum(abs(w))/12
计算结果如表所示
表 计算值与实验原始数据值对照表 1,2,
,12=k
由表可以看出，最大误差高达%，最低的也达到%，模拟误差较大，进一步
计算平均相对误差
13(0)
2
1()28.01%12=∆==∑k q k 平均相对误差很较大，相对精度约70%。

因此为了提高远原点（即现时）精度，即将最后一个误差减小，需采用残差模型进行修正。

第三步，以部分残差数据为原始数据建立新的GM （1，1）模型
取09k =得残差尾端，即取最后5个数据的残差：,,,,,
用此尾段可建立残差尾段模型，取绝对值，得残差数列
()(0)(0)
0.0790,0.0253,0.0374,0.0103,0.0190==g q
以上述的残差数列为原始数据建立新的GM （1，1）模型，得残差的时间响应式
()(1)
0.189410.17320.2522-+=-+k g
k e
()(0)
0.189410.0328-+=k g
k e
第四步，将原始数据和部分残差数据的两个GM （1，1）模型即
()(0)
0.0648610.368-+=k x
k e 和
()(0)
0.145710.1876-+=k g
k e
结合，得到修正后的残差GM （1，1）模型
00.064860.064860.18940.368,9
(1)0.3680.0328,9-∧
--⎧<⎪
+=⎨⎪-≥⎩
k k k e k x k e e k
第五步，用修正后的模型对8,9,
,12=k 的模拟值进行修正，结果为：
(0)(0)(0)ˆˆˆ((9),(10),,(13))0.2118,0.1993,0.1874,0.1762,0.1656=x
x x
第六步，精度检验 建立如下程序：
function 10toliti021(X0) %format long ;
%X0=*[6 20 40 25 40 45 35 21 14 18 17 15]; [m,n]=size(X0); s(1)=1; for i=8:12
y(i+1)=*exp*i)*exp*i); z(i+1)=X0(i+1)-y(i+1); w(i+1)=z(i+1)/X0(i+1); s(i+1)=i+1; end y' X0' z' w'
z*z'
sum(abs(w))/5
计算结果如表所示
表 修正后计算值与实验原始数据值检验结果 8,9,
,12=k
按此模型，可对9,10,11,12,13=k 五个模拟值进行修正，修正后的平均相对误差
13(0)
9
1()19.4%5=∆==∑k q k ，精度有明显的提高。

尤其对于原点附近的两个数据：，
相对误差分别降低为%和%，低于允许误差要求。

这说明，对原点数据GM （1，1）模型修正是有必要的。

模型在复杂问题中的应用——SARS 疫情问题
例 SARS 疫情问题
2003年的SARS 疫情对我国的发展产生了一定影响，尤其在经济发展方面产生了很大的影响，下面就SARS 疫情对我国经济的影响问题建立GM 模型并求解。

问题的提出
2003年的SARS 疫情对中国部分行业的经济发展产生了一定影响，特别是对部分疫情较严重的省市的相关行业所造成的影响是显著的，经济影响主要分为直接经济影响和间接影响。

直接经济影响涉及商品零售业、旅游业、综合服务等行业。

很多方面难以进行定量地评估，现仅就SARS 疫情较重的某市商品零售业、旅游业和综合服务业的影响进行定量的评估分析。

究竟 SARS 疫情对商品零售业、旅游业和综合服务业的影响有多大，已知某市从1997 年
1 月到2003 年1
2 月的商品零售额、接待旅游人数和综合服务收入的统计数据如表—所示
表商品的零售额（亿元）
表接待海外旅游人数（万人）
表综合服务业累计数额
（亿元）
业、旅游业和综合服务业所造成的影响。

模型的假设与分析
模型假设：
（1）假设该市的统计数据都是可靠准确的；
（2）假设该市在SARS 疫情流行期间和结束之后，数据的变化只与SARS 疫情的影响有关，不考虑其它随机因素的影响。

模型分析：
根据所掌握的历史统计数据可以看出，在正常情况下，全年的平均值较好地反映了相关指标的变化规律，这样可以把预测评估分成两部分：
（1）利用灰色理论建立GM(1,1)模型，由1997－2002 年的平均值预测2003 年平均值； （2）通过历史数据计算每个月的指标值与全年总值的关系，从而可预测出正常情况下2003 年每个月的指标值，再与实际值比较可以估算出SARS 疫情实际造成的影响。

建立灰色预测模型GM(1,1) 第一步,数据处理
（1）原始数据：
根据表中的已知数据，计算1997－2002年某项指标的年平均值，作为原始数据，记为
(0)(0)(0)(0)((1),(2),
,(6))=x x x x
并要求级比(0)
(0)()(1)/()(0.7515,1.3307)(1,2,
,6)λ=-∈=i x
i x i i
（2）数据的累加生成： 对原始数据(0)
x
进行一次累加生成，
(1)(0)(1)(1)=x x ，
1
(1)(1)(0)
(0)1
(1)[()(1)](),1,2,
,1+=+=++==-∑k i x k x k x
k x i k n
因此累加生成数据(1)
x 记为
(1)(0)(1)(1)(1)((1),(2),
,(6))==x AGOx x x x
（2）背景值的选择：
取累加生成数据(1)
x 的加权平均值为背景值(1)z ，即
(1)(1)(1)(1)(1)(1)(),(1,2,
,5)αα+=++-=z k x k x k k
其中α为确定参数。

背景值(1)z 记为
(1)(1)(1)(1)((2),(3),
,(6))=z x x x
第二步，GM （1，1）模型的建立
（1）建立GM （1，1）的白化微分方程模型
(1)
(1)+=dx ax u dt
其中是a 发展灰度，u 是内生控制灰度。

（2）转化为灰微分方程
(0)(1)()(),(2,3,
,6)+==x k az k u k
或
(0)(1)()(),(2,3,
,6)=-+=x k az k u k
即矩阵形式为
⎡⎤
==⎢⎥⎣⎦
a Y B BA u
其中
(1)(1)(0)(0)(0)(1)(2)1(3)1((2),(3),
,(6)),,(6)1⎡⎤-⎢⎥-⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦⎢⎥-⎣⎦
z a z Y x x x B A u z (3)转化为时间响应函数
利用最小二乘法得到参数的估计值,a u ，进而得到灰微分方程的解(1)
x ，对(1)
x
求导
还原得(0)
x。

即参数的估计值,a u 为
1()-⎡⎤
==⎢⎥⎣⎦
T T a A B B B Y u
微分方程的解式（也称时间响应函数）为
(1)(1)(1)((1))-+=-+ak u u
x k x e a a
(0)(0)(1)((1))-+=--ak u
x k a x e a
其中(0)
(1)(1)(1)=x
x ，(0)(1)+x k 称为还原值。

第三步，利用模型预测指标值
根据时间响应函数可以预测出正常情况下2003年的平均值，则x ，则预测2003年的总值为12=Z x 。

根据历史数据，可以计算出2003 年第i 个月的指标值占全年总值的比例
为i u ，即
6
112
6
1
1
,(1,2,
,12)====
=∑∑∑ij
j i ij
i j a
u i a
记为1212(,,,)=u u u u ，于是可以可到2003年每一个月的指标值=V Zu
模型求解及结果分析 （1）商品零售额
根据商品零售额的数据表，计算得年平均值（即原始数据(0)
(1)x
）和一次累加生成值
(1)(1)x ，分别为
(0)(1)(87.6167, 98.5,108.475,118.4167,132.8083,145.4083)=x (1)(1)(87.6167,186.1167, 294.5917, 413.0083, 545.8167, 691.225)=x
显然(0)
x
的所有级比都在可容区域内，这里取0.5α=，计算可得背景值
(1)(136.8667,240.3542,353.8000,479.4125,618.5208)=z
计算得参数的估计值为0.0983,84.7563=-=a u ，进而得到时间响应函数
(1)0.0983(1)949.6443862.0276+=-k x k e
(0)0.0983(1)93.3710+=k x k e
再根据时间响应函数预测可得，2003年的月平均值为160.4135=x 亿元；年总值为
121925.0==Z x 亿元。

又根据比例i u 的表达式计算出每月的比例为
(0.0794, 0.0807, 0.0749, 0.0786, 0.0819, 0.0818, 0.0845, 0.0838, 0.0872, 0.0886, 0.0866, 0.092)
=i u
因此2003 年1~12 月的预测值(单位：亿元)为
(152.8654,155.3486,144.1859,151.2177,157.7157,157.4140,
162.5660,161.3128,167.9501,170.5260,166.7433,177.1169)
==V Zu
将预测值与实际值进行比较，结果如表所示
表 2003年商品的零售额比较表 （亿元）
图形如图10-1：（蓝线为实际值，红线为预测值）
024681012
120
130
140
150
160
170
180
190
图10-1 2003年商品的零售额实际值与预测值比较图
通过图形可以直观的看出：（1）预测值波动比较小，真实值波动比较剧烈；（2）5月份左右真实值远远低于预测值，年初和年末都高于预测值。

这是由实际情况造成的，年初当
SARS 疫情刚刚开始的时候，人们储备保健药品和保健食物等，拉动了零售额的增长；5月份左右，SARS 疫情比较猖獗，此时好多学校和单位等实行封闭管理，大大限制了人们的消费，因此零售额明显降低；年末SARS 疫情慢慢远去，此前被限制的消费得以充分实现，又促进了零售额的增长。

当然可以根据模型所得数据，对SARS 疫情给该市的商品零售业造成的影响进行定量分析，这里不再详述。

计算的MATLAB 程序如下：
function 10toliti03
clc
clc,clear
load %把原始数据保存在纯文本文件中
han1(end,:)=[];
x0=mean(han1,2)
m=size(han1,2);
n=size(x0,1);
z1=[];
x1=cumsum(x0)
alpha=;
for i=1:n-1
z1(i,:)=(1-alpha)*x1(i)+alpha*x1(i+1);
end
z1
Y=x0(2:n);
B=[-z1,ones(n-1,1)];
A=inv(B'*B)*B'*Y;
a=A(1)
u=A(2)
b4=u/a
b5=x1(1)-b4
b6=-a*b5
k=6;
x7hat=(x0(1)-u/a)*(exp(-a*k)-exp(-a*(k-1))) z=m*x7hat
u=sum(han1)/sum(sum(han1)) v=z*u
（2）接待海外旅游人数
根据接待海外旅游人数的数据表，计算得年平均值（即原始数据(0)
(1)x ）和一次累加
生成值(1)
(1)x ，分别为
(0)(1)(19.1000,18.1083,20.8333,24.3917,24.7500,27.1750)=x (1)(1)(19.1000,37.2083,58.0417,82.4333,107.1833,134.3583)=x
显然(0)
x
的所有级比都在可容区域内，这里取0.5α=，计算可得背景值
(1)( 28.1542,47.6250,70.2375,94.8083,120.7708)=z
计算得参数的估计值为0.0938,16.2671=-=a u ，进而得到时间响应函数
(1)0.0938(1)192.4955173.3955+=-k x k e
(0)0.0938(1)192.4955+=k x k e
再根据时间响应函数预测可得，2003年的月平均值为30.2649=x 万人；年总值为
12363.1785==Z x 万人。

又根据比例i u 的表达式计算出每月的比例为
(0.0407, 0.0732, 0.0703, 0.0878, 0.0907, 0.0848, 0.0836, 0.1022, 0.101, 0.1041, 0.0914, 0.0701)
=i u
因此2003 年1~12 月的预测值(单位：万人)为
( 14.7992,26.5801,25.5439,31.8961,32.9548,30.7923,
30.3644,37.1220,36.6715,37.7978,33.1800,25.4763)
==V Zu
将预测值与实际值进行比较，结果如表所示
表 2003年接待海外旅游人数 （万人）
图形如图10-2：（蓝线为实际值，红线为预测值）
024681012
5101520253035
40
图10-2 2003年接待海外旅游人数实际值与预测值比较图
通过图形可以直观的看出：（1）预测值波动比较小，真实值波动比较剧烈；（2）真实值低于预测值，尤其5月份左右真实值远远低于预测值，年初和年末相差不太大。

这是由实际情况造成的， 5月份左右SARS 疫情比较猖獗，此时好多学校和单位等实行封闭管理，大大限制了人们的出行，同时人们也基于自身安全的因素，能不出门就不出门，因此旅游人数大大降低，旅游业处于低谷；年初和年末SARS 疫情对人们出行的影响不大，因此年初和年末年末海外旅游人数的实值略低于预测值。

（3）综合服务业累计数据
根据综合服务业累计数据的数据表，计算得年平均值（即原始数据(0)
(1)x ）和一次累
加生成值(1)
(1)x ，分别为
(0)(1)(483.3,588.2,657.8,778.4,874.9,1000.9)=x (1)(1)(483.3,1071.5,1729.3,2507.6,3382.5,4383.5)=x
显然(0)
x
的所有级比都在可容区域内，这里取0.5α=，计算可得背景值
(1)(777.4,1400.4,2118.5,2945.1,3883.0)=z
计算得参数的估计值为0.1343,481.2013=-=a u 。

进而得到时间响应函数
(1)0.1343(1)406.59358.26+=-k x k e
(0)0.1343(1)546.1129+=k x k e
再根据时间响应函数预测可得，2003年的月平均值为1144.0=x 亿元；年总值为
111258.4==Z x 亿元。

又根据比例i u 的表达式计算出每月的比例为
(0.0191, 0.031, 0.0433, 0.0591, 0.0728, 0.085, 0.1046, 0.1205, 0.1358, 0.1515, 0.1749)
=i u
因此2003 年1~12 月的预测值(单位：亿元)为
(240.1,389.7,545.2,743.8,915.8,1100.9,
1316.2,1516.6,1708.7,1906.5,2200.9)
==V Zu
将预测值与实际值进行比较，结果如表所示
表 2003年综合服务业累计数据的比较表 （亿元）
图形如图10-3：（蓝线为实际值，红线为预测值）
图10-3 2003年综合服务业累计数据实际值与预测值比较图
通过图形可以直观的看出：（1）预测值与真实值相差不大；（2）5月份左右真实值与预测值最接近。

这是由实际情况造成的，SARS 疫情对于综合服务业中的部分行业影响较大，如航空交通运输、宾馆餐饮等，但有些行业影响不大，如电信、通讯等，因此总平均来看，影响还不算太大。

（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）模型表示一阶的含有两个变量灰色模型，其相应的白化微分方程模型为
(1)
(1)(1)112+=dx ax ux dt
时间相应函数(即解)为
(1)(0)(1)(1)1122(1)[(1)(1)](1)-+=-
+++ak u u
x k x x k e x k a a
还原值为
(0)(1)(0)111(1)(1)()+=+-x k x k x k
例 某系统中两因素1x 和2x 相互之间存在关系，其中因素1x 为系统特征因素，因素
2x 为相关因素，两因素的关系如表所示
表 两因素的关系
解：第一步，数据处理 取原始数据为
(0)1(10.7,12.3,15.4,19.2,21.5)=x
(0)2(40.3,50.5,60.3,65.5,75.2)=x
对原始数据累加生成，得
(1)(0)11(10.7,23,38.4,57.6,79.1)==x AGOx (1)(0)22(40.3,90.8,151.1,216.6,291.8)==x AGOx
第二步，建立GM （1，2）模型 白化的微分方程模型为
(1)
(1)(1)112+=dx ax ux dt
转化为灰微分方程
(0)(1)(1)111()()(),(2,3,
,5)+==x k az k uz k k
或
(0)(1)(1)111()()(),(2,3,
,5)=-+=x k az k uz k k
其中
(1)(1)(1)1
()(1)(),(2,3,4,5)2
+-==x k x k z
k k 。

即矩阵形式为
⎡⎤
==⎢⎥⎣⎦
a Y B BA u
其中
(1)(1)12(1)
(1)(0)(0)(0)12(1)
(1)12(1)
(1)12(2)(2)(3)
(3)((2),(3),
,(5)),,(4)(4)(6)(5)⎡⎤
-⎢⎥
-⎡⎤
⎢⎥===⎢⎥⎢⎥-⎣⎦
⎢⎥
-⎢⎥
⎣⎦z x a z x Y x x x B A u z x z x
(3)时间响应函数
利用最小二乘法得到参数的估计值,a u ，即
1()-⎡⎤
==⎢⎥⎣⎦
T T a A B B B Y u
进而得到微分方程的解即时间响应函数(1)
1x 及还原值(1)
0x ，即
(1)(0)(1)(1)1122(1)[(1)(1)](1)-+=-
+++ak u u
x k x x k e x k a a
(0)(1)(0)111(1)(1)()+=+-x k x k x k
第三步，模型的求解
建立M 文件如下：
clc,clear x10=[,,,,]; x20=[,,,,]; n=length(x10); x11=cumsum(x10) x21=cumsum(x20) for i=2:n
z11(i)=*(x11(i)+x11(i-1)); end
B=[-z11(2:n)',x21(2:n)']; Y=x10(2:n)'; A=inv(B'*B)*B'*Y
x=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)}) xhat=[x(1),diff(x)] epsilon1=x11-x
delta1=abs(epsilon1./x11) epsilon0=x10-xhat delta0=abs(epsilon0./x10)
利用程序10toliti04求解,得 1.0975,0,3307==a u ，进而得到时间响应函数及还原值
(1)(1) 1.0975(1)122(1)(10.70.0310(1)0.3013(1)-+=-+++k x k x k e x k
(0)1(10.7,21.7991,41.6482,63.2342,86.9622)=x
第四步，精度检验
先对累加生成数据进行的误差检验，如果进度较高可直接求解还原值，否则要进行模型的修正，检验结果如表所示
表累加生成数据的误差检验表
根据误差检验表，可知：累加生成数据相对误差都小于10%，最大的相对误差为%，模型的相对误差比较小，精度较高；可以直接用来求解还原值。

利用累加生成数据进行累减生成，得到还原值(0)
1()
x k，即
(0) 1()(10.7000,11.0991,19.8491,21.5860,23.7280) =
x k
并对还原值进行的误差检验，检验结果如表所示
表原始数据的误差检验表
根据误差检验表，可知：原始数据的相对误差集中在10%左右，最大的相对误差为%，模型的相对误差比较小，精度较高。

GM（1，n）模型建模过程与GM（1，2）模型类似，不再详细介绍。

本章小结
本章主要介绍了灰色系统的基本概念和三个基本的灰色模型：GM（1，1）模型、残差
GM（1，1）模型和GM（1，2）模型，其中GM（1，1）模型是其他模型的理论基础，会准确熟练的建立GM（1，1）并求解，并在此基础上进行深入学习。

思考题
思考题1 某江段连续四年氨氮质量浓度如表所示，试建立GM（1，1）模型并求解。

表某江段氨氮质量浓度 mg/L
思考题2 某市某月1-15日早晨7点左右的大气能见度值如表所示，试建立GM(1,1)残差模型并求解。

表某市大气能见度
思考题3 设系统特征数据序列为(0)
1(2.874, 3.278, 3.307, 3.39, 3.679)
x=
相关因素数据序列为(0)
2(7.04, 7.645, 8.075, 8.53, 8.774)
x=，试建立GM(1,2)模型并求解。

推荐阅读书目
1．刘思峰等．2004．灰色系统理论及应用［Ｍ］．3版．北京：科学出版社．2．王治祯等．2001．灰色系统及模糊数学在环境保护中的应用［Ｍ］．1版．哈尔滨：哈尔滨工业大学出版社．
3．邓聚龙．1990．多维灰色规划［Ｍ］．1版．北京：华中理工大学出版社．。