灰色模型介绍及应用

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

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

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

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

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

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

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

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

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

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

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

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

灰数：一个信息不完全的数，称为灰数。

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

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

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

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

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

则称此系统为灰色系统。

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

累加生成一般可写成AGO 。

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

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

当∆t 足够小时，()x t 到()+∆x t t 不会发生突变，因此可取()x t 与()+∆x t t 的平均值作为0∆→t 时的背景值，因此，背景值便可记为或于是白化的微分方程(1)(1)+=dx ax u dt可改写为 或 即因此，上述方程可以改写为矩阵方程形式，即引入下列符号，设 于是便有 令 则 解得将求解得到的代入微分方程的解式（也称时间响应函数），则 由于(0)(1)(1)(1)=xx ，因此求导还原得上述两式便为GM （1，1）的时间响应式，及灰色系统预测模型的基本算式，当然上述两式计算结果只是近似计算值。

为简记，一般可以将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 ）解：第一步，设原始数据为(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因此累加生成数据为第三步，构造矩阵,N B Y(0)(0)(0)[(2),(3),,(6)][4.575 4.968 5.063 5.968 5.507]==T TN Y x x x第四步，计算1ˆ()-=TT NaB B B Y 。

先求1()-TB B ，即根据逆矩阵的求解方法，得再求TN B Y 的值，即进而求得ˆa的值为 计算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) 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==xx ，得第六步，求导还原得第七步，对上述模型进行精度检验。

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

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

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

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

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

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

令(0)()g k 为残差，则即 或利用残差序列(0)g建立新的GM （1，1）模型，求解得时间响应式第三步，结合上两步的GM （1，1）模型，建立残差GM （1，1）模型 结合上两步的GM （1，1）模型，则相应的残差修正时间响应式为 称为导数还原式的残差修正模型。

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

表10.4 某县油菜发病率数据 （%）解：第一步，建立原始数据的GM （1，1）模型 设原始数据为建立GM （1，1）模型，利用GM （1，1）的求解程序得时间响应式为第二步，误差检验 利用时间响应函数模型()(0)0.0648610.368-+=k x k e 计算各时刻值（1,2,,12=k ），并计算相对误差，程序如下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:12y(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; end y' X0' z' w' z*z'sum(abs(w))/12计算结果如表10.5所示表10.5 计算值与实验原始数据值对照表 1,2,,12=k由表可以看出，最大误差高达72.44%，最低的也达到6.06%，模拟误差较大，进一步计算平均相对误差平均相对误差很较大，相对精度约70%。

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

第三步，以部分残差数据为原始数据建立新的GM （1，1）模型取09k =得残差尾端，即取最后5个数据的残差：-0.0790,-0.0253,-0.0374,-0.0103,-0.0190, 用此尾段可建立残差尾段模型，取绝对值，得残差数列以上述的残差数列为原始数据建立新的GM （1，1）模型，得残差的时间响应式第四步，将原始数据和部分残差数据的两个GM （1，1）模型即 和结合，得到修正后的残差GM （1，1）模型第五步，用修正后的模型对8,9,,12=k 的模拟值进行修正，结果为：第六步，精度检验 建立如下程序：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:12y(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; end y' X0' z' w' z*z'sum(abs(w))/5计算结果如表10.6所示表10.6 修正后计算值与实验原始数据值检验结果 8,9,,12=k按此模型，可对9,10,11,12,13=k 五个模拟值进行修正，修正后的平均相对误差13(0)91()19.4%5=∆==∑k q k ，精度有明显的提高。