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）模型的应用因素相关问题本章小结思考题推荐阅读书目2第十章 灰色模型介绍及应用10.1 灰色理论基本知识10.1.1 概言客观世界的很多实际问题，其内部的结构、参数以及特征并未全部被人们了解，人们不可能象研究白箱问题那样将其内部机理研究清楚，只能依据某种思维逻辑与推断来构造模型。对这类部分信息已知而部分信息未知的系统，我们称之为灰色系统。本章介绍的方法是从灰色系统的本征灰色出发，研究在信息大量缺乏或紊乱的情况下，如何对实际问题进行分析和解决。灰色系统的研究对象是“部分信息已知、部分信息未知”的“小样本”、“贫信息”不确定性系统，它通过对“部分”已知信息的生成、开发实现对现实世界的确切描述和认识。信息不完全是“灰”的基本含义。灰色系统理论建模的主要任务是根据具体灰色系统的行为特征数据，充分开发并利用不多的数据中的显信息和隐信息，寻找因素间或因素本身的数学关系。通常的办法是采用离散模型，建立一个按时间作逐段分析的模型。但是，离散模型只能对客观系统的发展做短期分析，适应不了从现在起做较长远的分析、规划、决策的要求。尽管连续系统的离散近似模型对许多工程应用来讲是有用的，但在某些研究领域中，人们却常常希望使用微分方程模型。事实上，微分方程的系统描述了我们所希望辨识的系统内部的物理或化学过程的本质。目前，灰色系统理论已成功地应用于工程控制、经济管理、未来学研究、生态系统及复杂多变的农业系统中，并取得了可喜的成就。灰色系统理论有可能对社会、经济等抽象系统进行分析、建模、预测、决策和控制，它有可能成为人们认识客观系统改造客观系统的一个新型的理论工具。10.1.2 有关名词概念灰数：一个信息不完全的数，称为灰数。灰元：信息不完全或内容难以穷尽的元素，称为灰元。灰关系：信息不完全或机制不明确的关系，称为灰关系。具有灰关系的因素是灰因素，3灰因素之间的量化作用，称为灰关联。灰色系统:含灰数、灰元或灰关系的系统称为信息不完全系统。如果按照灰色理论去研究它。则称此系统为灰色系统。累加生成：由于灰系统对一切随机量都可看作是在一定范围内变化的灰色量，因此，为适应灰系统建模需要，提出“生成”的概念， “生成”即指对原始数据做累加（或累减）处理。累加生成一般可写成 AGO。若计为原始数列，为次累加生成后数列，即(0)x( )rxr(0)(0)(0)(0)(1),(2),( )xxxxnL( )( )( )( )(1),(2),( )rrrrxxxxnL则次累加生成算式为r( )(1)(1)(1)(1)1(1)(1)(1)(1)( )(1)( )(1)(2)( )( )(1)(2)(1)( )(1)( )k rrrrrirrrrrrxkxxxkxixxxkxkxkxkLL一般常用的是一次累加生成，即(1)(0)1(1)(0)( )( )(1)( )kixkxixkxk10.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）模型的特例，其简单的微分方程形式（白化形式的微分方程）是dxaxudt利用常数变易法解得，通解为( )atux tcea4若初始条件为，则可得到微分方程的特解为00, ( )tx tx0( )()atuux txeaa或时间响应函数(1)(1)(1)atuux txeaa其中白化微分方程中的项中的为的背景值，也称为初始值; 为常数axxdx dt,a u（有时也将写成） 。ub按白化导数定义有差分形式的微分方程，即0()( )lim tdxx ttx t dtt显然，当时间密化值定义为 1，即当时，上式可记为1 t1 (1)( )lim tdxx tx tdt记为离散形式(1)( )dxx tx tdt这显然表明是一次累计生成，因此上述方程可改写为dx dt(1)(1)(0)(1)( )(1)dxxtxtxtdt这实际也表明，模型是以生成数（是以的一次累加）为基础的。(1)x(1)x(0)x当足够小时，到不会发生突变，因此可取与的平均t( )x t()x tt( )x t()x tt值作为时的背景值，因此，背景值便可记为0 t(1)(1)(1)1(1)( )2xxtxt或(1)(1)(1)1(1)( )2xxkxk于是白化的微分方程可改写为(1) (1)dxaxudt(0)(1)(1)1(1)(1)( )2xka xkxku5或(0)(1)(1)1(1)(1)( )2 xka xkxku即(0)(1)(1)(0)(1)(1)(0)(1)(1)1(2)(2)(1)2 1(3)(2)(1)21( )( )(1)2 Mxa xxuxa xxuxna xnxnu因此，上述方程可以改写为矩阵方程形式，即(1)(1)(0)(1)(1)(0)(0)(1)(1)1(2)(1)21(2)1(2)(1)1(3)21( )1( )(1)2 MMMa xx x a xxxauxn a xnxn引入下列符号，设(0)(0)(0)(2) (3)( ) MNx xYxn1 11 ME(1)(1)(1)(1)(1)(1)1(2)(1)2 1(2)(1)21( )(1)2 Ma xxa xxXa xnxn于是便有 MNaYaXuEX Eu令6 )aau(1)(1)(1)(1)(1)(1)1(2)(1)12 1(2)(1)121( )(1) 12 M MMa xxa xxBX Ea xnxn则 )MNaYaXuEX EBau解得1() )TT NaaB BB Yu将求解得到的代入微分方程的解式（也称时间响应函数） ，则(1)(1)(1)(1)akuuxkxeaa由于，因此求导还原得(0)(1)(1)(1)xx(0)(0)(1)(1) )akuxka xea上述两式便为 GM（1，1）的时间响应式，及灰色系统预测模型的基本算式，当然上述两式计算结果只是近似计算值。为简记，一般可以将 GM（1，1）的建模过程记为(0)0(1)(0)(1); , )(1)(1)IAGO GM AGOxGM xa uxkxk10.2 灰色理论模型应用10.2.1GM（1，1）模型的应用污染物浓度问题GM（1，1）模型是灰色系统最基本的模型，下面以污染物浓度问题说明 GM（1，1）模型的建立及求解过程。例例 10.110.1 某污染源中某种污染物质量浓度测量值如表 10.1，试建立 GM（1，1）模型表表 10.1 某污染物质量浓度测量值某污染物质量浓度测量值 （mg/L） 年 份2001200220032004200520067( )0x3.9364.5754.9685.0635.9685.507解：第一步，设原始数据为(0)(0)(0)(0)(1),(2),(6)(3.936,4.575,4.968,5.063,5.968,5.507)Lxxxx第二步，对原始数据进行累加生成，即(1)(0)xAGOx(1)(0)(1)(1)3.936xx(1)(1)(0)(2)(1)(2)3.9364.5758.511xxx(1)(1)(0)(3)(2)(3)13.479xxx(1)(1)(0)(4)(3)(4)18.542xxx(1)(1)(0)(5)(4)(5)24.510xxx(1)(1)(0)(6)(5)(6)30.017xxx因此累加生成数据为(1)(0)(3.936,8.511,13.479,18.542,24.510,30.017)xAGOx第三步，构造矩阵,NB Y(1)(1)(1)(1)(1)(1)(1)(1)(1)(1)1(1)(2)12 1-6.2235 1.0000(2)(3)12-10.9950 1.00001-16.0105 1.0000(3)(4)12-21(4)(5)12 1(5)(6)12 xxxxBxxxxxx1.5260 1.0000-27.2635 1.0000 (0)(0)(0)(2),(3),(6)4.575 4.968 5.063 5.968 5.507LTT NYxxx第四步，计算。1()TT NaB BB Y先求，即1()TB B81622.6 -82.0-82.0 5TB B根据逆矩阵的求解方法，得10.0036 0.0592()0.0592 1.1706TB B再求的值，即T NB Y-442.764126.0810T NB Y进而求得的值为 a10.0036 0.0592-442.7641-0.0539a()0.0592 1.1706 26.0810 4.3322u TT NaB BB 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); end B=-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.936xx9(1)(1)0.0539(1)(1)84.326480.3904akkuuxkxeeaa第六步，求导还原得(0)(1)0.0539(1)