灰色模型介绍及应用

第十章灰色模型介绍及应用（徐利艳天津农学院 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.3GM模型在复杂问题中的应用——SARS 疫情问题 10.2.4 GM（1，n）模型的应用——因素相关问题本章小结思考题推荐阅读书目第十章灰色模型介绍及应用10.1灰色理论基本知识10.1.1概言客观世界的很多实际问题，其内部的结构、参数以及特征并未全部被人们了解，人们不可能象研究白箱问题那样将其内部机理研究清楚，只能依据某种思维逻辑与推断来构造模型。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

则称此系统为灰色系统。

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

累加生成一般可写成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)()kr r r r r i r r r r r r x k xxxk x i x x x k x k x k x k ----=-----=++==++-+=-+∑ 一般常用的是一次累加生成，即(1)(0)1(1)(0)()()(1)()ki x k x i x k x k ===-+∑10.1.3GM 建模机理建立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 ）模型的特例，其简单的微分方程形式（白化形式的微分方程）是+=dxax u dt利用常数变易法解得，通解为()-=+at u x t ce a若初始条件为00,()==t x t x ，则可得到微分方程的特解为0()()-=-+at u ux t x e a a或时间响应函数(1)(1)((1))-+=-+at u ux t x e a a其中白化微分方程中的ax 项中的x 为dxdt的背景值，也称为初始值; ,a u 为常数（有时也将u 写成b ）。

按白化导数定义有差分形式的微分方程，即()()lim ∆→+∆-=∆t dx x t t x t dt t 显然，当时间密化值定义为1，即当1∆→t 时，上式可记为1[(1)()]lim ∆→=+-t dxx t x t dt 记为离散形式(1)()=+-dxx t x t dt这显然表明dxdt是一次累计生成，因此上述方程可改写为 (1)(1)(0)(1)()(1)=+-=+dxx 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 ux 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[]21[()(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 ux k x e a a由于(0)(1)(1)(1)=xx ，因此求导还原得(0)(0)(1)((1))-+=--ak ux 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 k10.2灰色理论模型应用10.2.1GM （1，1）模型的应用——污染物浓度问题GM （1，1）模型是灰色系统最基本的模型，下面以污染物浓度问题说明GM （1，1）模型的建立及求解过程。

例10.1 某污染源中某种污染物质量浓度测量值如表10.1，试建立GM （1，1）模型表10.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 TN Y x x x 第四步，计算1ˆ()-=T T N aB B B Y 。

先求1()-TB B ，即1622.6 -82.0 -82.0 5⎡⎤=⎢⎥⎣⎦T B B 根据逆矩阵的求解方法，得1 0.0036 0.0592() 0.0592 1.1706-⎡⎤=⎢⎥⎣⎦T B B 再求TN B Y 的值，即-442.7641 26.0810⎡⎤=⎢⎥⎣⎦T N B Y 进而求得ˆa的值为 10.0036 0.0592-442.7641-0.0539a ˆ() 0.0592 1.1706 26.0810 4.3322u -⎡⎤⎡⎤⎡⎤⎡⎤====⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦TTN a B B B Y 计算GM1_1的程序如下function 10toliti01(X0) [m,n]=size(X0); X1=cumsum(X0); X2=[];for i=1:n-1X2(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)-b4 b6=-a*b5第五步，将,a u 的值代入微分方程的时间响应函数，令(1)(1)ˆ(1)(0) 3.936==xx ，得(1)(1)0.0539ˆ(1)((1))84.326480.3904-+=-+=-ak k u uxk x e e a a第六步，求导还原得(0)(1)0.0539ˆ(1)((1)) 4.5443-+=--=ak k uxk a x e e a第七步，对上述模型进行精度检验。