灰色模型介绍及应用

灰色模型介绍及应用PAGEPAGE13第十章灰色模型介绍及应用（徐利艳天津农学院2.4万字）10.1灰色理论基本知识10.1.1概言10.1.2有关名词概念10.1.3GM建模机理10.2灰色理论模型应用10.2.1GM（1，1）模型的应用——污染物浓度问题10.2.2GM（1，1）残差模型的应用——油菜发病率问题10.2.3GM模型在复杂问题中的应用——SARS疫情问题10.2.4GM（1，n）模型的应用——因素相关问题本章小结思考题推荐阅读书目第十章灰色模型介绍及应用10.1灰色理论基本知识10.1.1概言客观世界的很多实际问题，其内部的结构、参数以及特征并未全部被人们了解，人们不可能象研究白箱问题那样将其内部机理研究清楚，只能依据某种思维逻辑与推断来构造模型。对这类部分信息已知而部分信息未知的系统，我们称之利用常数变易法解得，通解为若初始条件为，则可得到微分方程的特解为或时间响应函数其中白化微分方程中的项中的为的背景值，也称为初始值;为常数（有时也将写成）。按白化导数定义有差分形式的微分方程，即显然，当时间密化值定义为1，即当时，上式可记为记为离散形式这显然表明是一次累计生成，因此上述方程可改写为这实际也表明，模型是以生成数（是以的一次累加）为基础的。当足够小时，到不会发生突变，因此可取与的平均值作为时的背景值，因此，背景值便可记为或于是白化的微分方程可改写为或即因此，上述方程可以改写为矩阵方程形式，即引入下列符号，设于是便有令则解得将求解得到的代入微分方程的解式（也称时间响应函数），则由于，因此求导还原得上述两式便为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）年份2001200220032004200520063.9364.5754.9685.0635.9685.507解：第一步，设原始数据为第二步，对原始数据进行累加生成，即因此累加生成数据为第三步，构造矩阵第四步，计算。先求，即根据逆矩阵的求解方法，得再求的值，即进而求得的值为计算GM1_1的程序如下function10toliti01(X0)[m,n]=size(X0);X1=cumsum(X0);X2=[];fori=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)-b4b6=-a*b5第五步，将的值代入微分方程的时间响应函数，令，得第六步，求导还原得第七步，对上述模型进行精度检验。常用的方法是回代检验，即分别用模型求出各时刻值，然后求相对误差。先利用时间响应函数模型求各时刻值（），并计算相对误差，结果如表10.2所示.表10.2精度检验实测值、残差值表GM计算值实测值残差相对残差8.605913.534418.735924.225430.01908.511013.479018.542024.510030.0170-0.0949-0.0554-0.19390.2846-0.0020-0.0112-0.0041-0.01050.0116-0.0001再利用时间响应函数模型求各时刻值（），并计算相对误差，结果如表10.3所示.表10.3计算值与实验原始数据值对照表GM计算值实测值残差相对残差4.79605.06165.34195.63775.94994.57504.96805.06305.96805.5070-0.2210-0.0936-0.27890.3303-0.4429-0.0483-0.0188-0.05510.0553-0.0804从残差检验结果看，累计生成数列曲线拟合较好，相对误差在0.01即1%左右；而还原数列的相对误差较大，其原因是累加生成数据将原始数据的随机性弱化，正负误差有抵消的，当数据再被还原回来时便表现出来。10.2.2GM（1，1）残差模型的应用——油菜发病率问题当ＧＭ（１，１）模型的精度不符合要求时，可用残差序列建立ＧＭ（１，１）模型，对原来的模型进行修正，以提高精度，即建立残差ＧＭ（１，１）模型，步骤如下第一步，利用原始数据建立GM（1，1）模型，得时间响应式其中第二个式子也成为导数还原值。鉴于导数还原值与原始数据（累减还原值）不一致，为减少往返运算造成的误差，往往用原始数据与导数还原值的残差修正的模拟值。第二步，利用残差数列建立新的GM（1，1）模型。建立残差模型的过程和计算方法同于GM（1，1）建模过程，只不过建立残差模型所用的原始数列采用的是残差数据。令为残差，则即或利用残差序列建立新的GM（1，1）模型，求解得时间响应式第三步，结合上两步的GM（1，1）模型，建立残差GM（1，1）模型结合上两步的GM（1，1）模型，则相应的残差修正时间响应式为称为导数还原式的残差修正模型。例10.2某县油菜发病率数据如表10.4所示，试建立残差GM模型并进行求解。表10.4某县油菜发病率数据（%）序号12345678910111213620402540453521141815.51715解：第一步，建立原始数据的GM（1，1）模型设原始数据为建立GM（1，1）模型，利用GM（1，1）的求解程序得时间响应式为第二步，误差检验利用时间响应函数模型计算各时刻值（），并计算相对误差，程序如下function10toliti02(X0)%formatlong;%X0=0.01*[620402540453521141815.51715];[m,n]=size(X0);s(1)=1;fori=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;endy'X0'z'w'z*z'sum(abs(w))/12计算结果如表10.5所示表10.5计算值与实验原始数据值对照表GM计算值实测值残差相对残差0.34490.32320.30290.28390.26610.24940.23370.21900.20530.19240.18030.16900.20.40.250.40.450.350.210.140.180.1550.170.15-0.14490.0768-0.05290.11610.18390.1006-0.0237-0.0790-0.0253-0.0374-0.0103-0.0190-0.72440.1919-0.21170.29020.40870.2875-0.1129-0.5645-0.1404-0.2412-0.0606-0.1265由表可以看出，最大误差高达72.44%，最低的也达到6.06%，模拟误差较大，进一步计算平均相对误差平均相对误差很较大，相对精度约70%。因此为了提高远原点（即现时）精度，即将最后一个误差减小，需采用残差模型进行修正。第三步，以部分残差数据为原始数据建立新的GM（1，1）模型取得残差尾端，即取最后5个数据的残差：-0.0790,-0.0253,-0.0374,-0.0103,-0.0190,用此尾段可建立残差尾段模型，取绝对值，得残差数列以上述的残差数列为原始数据建立新的GM（1，1）模型，得残差的时间响应式第四步，将原始数据和部分残差数据的两个GM（1，1）模型即和结合，得到修正后的残差GM（1，1）模型第五步，用修正后的模型对的模拟值进行修正，结果为：第六步，精度检验建立如下程序：function10toliti021(X0)%formatlong;%X0=0.01*[620402540453521141815.51715];[m,n]=size(X0);s(1)=1;fori=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;endy'X0'z'w'z*z'sum(abs(w))/5计算结果如表10.6所示表10.6修正后计算值与实验原始数据值检验结果GM计算值实测值残差相对残差0.21180.19930.18740.17620.16560.140.180.1550.170.15-0.0718-0.0193-0.0324-0.0062-0.0156-0.5130-0.1073-0.2093-0.0366-0.1040按此模型，可对五个模拟值进行修正，修正后的平均相对误差，精度有明显的提高。尤其对于原点附近的两个数据：0.17，0.15相对误差分别降低为3.66%和10.4%，低于允许误差要求。这说明，对原点数据GM（1，1）模型修正是有必要的。10.2.3GM模型在复杂问题中的应用——SARS疫情问题例10.3SARS疫情问题2003年的SARS疫情对我国的发展产生了一定影响，尤其在经济发展方面产生了很大的影响，下面就SARS疫情对我国经济的影响问题建立GM模型并求解。10.2.3.1问题的提出2003年的SARS疫情对中国部分行业的经济发展产生了一定影响，特别是对部分疫情较严重的省市的相关行业所造成的影响是显著的，经济影响主要分为直接经济影响和间接影响。直接经济影响涉及商品零售业、旅游业、综合服务等行业。很多方面难以进行定量地评估，现仅就SARS疫情较重的某市商品零售业、旅游业和综合服务业的影响进行定量的评估分析。究竟SARS疫情对商品零售业、旅游业和综合服务业的影响有多大，已知某市从1997年1月到2003年12月的商品零售额、接待旅游人数和综合服务收入的统计数据如表10.7—10.9所示表10.7商品的零售额（亿元）月份年份1月2月3月4月5月6月7月8月9月10月11月12月199719981999200020012002200383.079.878.185.186.688.290.386.793.392.590.996.9101.785.187.891.693.494.597.499.5104.2102.3101.0123.592.2114.093.3101.0103.5105.2109.5109.2109.6111.2121.7131.3105.0125.7106.6116.0117.6118.0121.7118.7120.2127.8121.8121.9139.3129.5122.5124.5135.7130.8138.7133.7136.8138.9129.6133.7137.5135.3133.0133.4142.8141.6142.9147.3159.6162.1153.5155.9163.2159.7158.4145.2124.0144.1157.0162.6171.8180.7173.5176.5表10.8接待海外旅游人数（万人）月份年份1月2月3月4月5月6月7月8月9月10月11月12月19971998199920002001200220039.411.316.819.820.318.820.924.924.724.319.418.69.611.715.819.919.517.817.823.321.424.520.115.910.112.917.721.021.020.421.925.829.329.823.616.511.426.019.625.927.624.323.027.827.328.532.818.511.526.420.426.128.928.025.230.828.728.122.220.713.729.723.128.929.027.426.032.231.432.629.222.915.417.123.511.61.782.618.816.220.124.926.521.8表10.9综合服务业累计数额（亿元）月份年份2月3月4月5月6月7月8月9月10月11月12月199719981999200020012002200396144194276383466554652747832972111169235400459565695805881101111391512383354255416417398669751087123816426337653160071191310381173129614971823184455767088561000114512921435166721636150464281897911421305147916441920241404584741923111412981492168418852218试根据这些历史数据建立预测评估模型，评估2003年SARS疫情给该市的商品零售业、旅游业和综合服务业所造成的影响。10.2.3.2模型的假设与分析模型假设：（1）假设该市的统计数据都是可靠准确的；（2）假设该市在SARS疫情流行期间和结束之后，数据的变化只与SARS疫情的影响有关，不考虑其它随机因素的影响。模型分析：根据所掌握的历史统计数据可以看出，在正常情况下，全年的平均值较好地反映了相关指标的变化规律，这样可以把预测评估分成两部分：（1）利用灰色理论建立GM(1,1)模型，由1997－2002年的平均值预测2003年平均值；（2）通过历史数据计算每个月的指标值与全年总值的关系，从而可预测出正常情况下2003年每个月的指标值，再与实际值比较可以估算出SARS疫情实际造成的影响。10.2.3.3建立灰色预测模型GM(1,1)第一步,数据处理（1）原始数据：根据表中的已知数据，计算1997－2002年某项指标的年平均值，作为原始数据，记为并要求级比（2）数据的累加生成：对原始数据进行一次累加生成，，因此累加生成数据记为（2）背景值的选择：取累加生成数据的加权平均值为背景值，即其中为确定参数。背景值记为第二步，GM（1，1）模型的建立（1）建立GM（1，1）的白化微分方程模型其中是a发展灰度，u是内生控制灰度。（2）转化为灰微分方程或即矩阵形式为其中(3)转化为时间响应函数利用最小二乘法得到参数的估计值，进而得到灰微分方程的解，对求导还原得。即参数的估计值为微分方程的解式（也称时间响应函数）为其中，称为还原值。第三步，利用模型预测指标值根据时间响应函数可以预测出正常情况下2003年的平均值，则，则预测2003年的总值为。根据历史数据，可以计算出2003年第i个月的指标值占全年总值的比例为，即记为，于是可以可到2003年每一个月的指标值10.2.4.4模型求解及结果分析（1）商品零售额根据商品零售额的数据表，计算得年平均值（即原始数据）和一次累加生成值，分别为显然的所有级比都在可容区域内，这里取，计算可得背景值计算得参数的估计值为，进而得到时间响应函数再根据时间响应函数预测可得，2003年的月平均值为亿元；年总值为亿元。又根据比例的表达式计算出每月的比例为因此2003年1~12月的预测值(单位：亿元)为将预测值与实际值进行比较，结果如表10.10所示表10.102003年商品的零售额比较表（亿元）月份1月2月3月4月5月6月7月8月9月10月11月12月实际值预测值163.2159.7158.4145.2124.0144.1157.0162.6171.8180.7173.5176.5152.9155.3144.2151.2157.7157.4162.6161.3168.0170.5166.7177.1图形如图10-1：（蓝线为实际值，红线为预测值）图10-12003年商品的零售额实际值与预测值比较图通过图形可以直观的看出：（1）预测值波动比较小，真实值波动比较剧烈；（2）5月份左右真实值远远低于预测值，年初和年末都高于预测值。这是由实际情况造成的，年初当SARS疫情刚刚开始的时候，人们储备保健药品和保健食物等，拉动了零售额的增长；5月份左右，SARS疫情比较猖獗，此时好多学校和单位等实行封闭管理，大大限制了人们的消费，因此零售额明显降低；年末SARS疫情慢慢远去，此前被限制的消费得以充分实现，又促进了零售额的增长。当然可以根据模型所得数据，对SARS疫情给该市的商品零售业造成的影响进行定量分析，这里不再详述。计算的MATLAB程序如下：function10toliti03clcclc,clearloadshuju1.txt%把原始数据保存在纯文本文件shuju1.txt中han1(end,:)=[];x0=mean(han1,2)m=size(han1,2);n=size(x0,1);z1=[];x1=cumsum(x0)alpha=0.5;fori=1:n-1z1(i,:)=(1-alpha)*x1(i)+alpha*x1(i+1);endz1Y=x0(2:n);B=[-z1,ones(n-1,1)];A=inv(B'*B)*B'*Y;a=A(1)u=A(2)b4=u/ab5=x1(1)-b4b6=-a*b5k=6;x7hat=(x0(1)-u/a)*(exp(-a*k)-exp(-a*(k-1)))z=m*x7hatu=sum(han1)/sum(sum(han1))v=z*u（2）接待海外旅游人数根据接待海外旅游人数的数据表，计算得年平均值（即原始数据）和一次累加生成值，分别为显然的所有级比都在可容区域内，这里取，计算可得背景值计算得参数的估计值为，进而得到时间响应函数再根据时间响应函数预测可得，2003年的月平均值为万人；年总值为万人。又根据比例的表达式计算出每月的比例为因此2003年1~12月的预测值(单位：万人)为将预测值与实际值进行比较，结果如表10.11所示表10.112003年接待海外旅游人数（万人）月份1月2月3月4月5月6月7月8月9月10月11月12月实际值预测值15.417.123.511.61.782.618.816.220.124.926.521.814.826.625.531.933.030.830.437.136.737.833.225.5图形如图10-2：（蓝线为实际值，红线为预测值）图10-22003年接待海外旅游人数实际值与预测值比较图通过图形可以直观的看出：（1）预测值波动比较小，真实值波动比较剧烈；（2）真实值低于预测值，尤其5月份左右真实值远远低于预测值，年初和年末相差不太大。这是由实际情况造成的，5月份左右SARS疫情比较猖獗，此时好多学校和单位等实行封闭管理，大大限制了人们的出行，同时人们也基于自身安全的因素，能不出门就不出门，因此旅游人数大大降低，旅游业处于低谷；年初和年末SARS疫情对人们出行的影响不大，因此年初和年末年末海外旅游人数的实值略低于预测值。（3）综合服务业累计数据根据综合服务业累计数据的数据表，计算得年平均值（即原始数据）和一次累加生成值，分别为显然的所有级比都在可容区域内，这里取，计算可得背景值计算得参数的估计值为。进而得到时间响应函数再根据时间响应函数预测可得，2003年的月平均值为亿元；年总值为亿元。又根据比例的表达式计算出每月的比例为因此2003年1~12月的预测值(单位：亿元)为将预测值与实际值进行比较，结果如表10.12所示表10.122003年综合服务业累计数据的比较表（亿元）月份2月3月4月5月6月7月8月9月10月11月12月实际值预测值241404584741923111412981492168418852218240390545744916110113161517170919072201图形如图10-3：（蓝线为实际值，红线为预测值）图10-32003年综合服务业累计数据实际值与预测值比较图通过图形可以直观的看出：（1）预测值与真实值相差不大；（2）5月份左右真实值与预测值最接近。这是由实际情况造成的，SARS疫情对于综合服务业中的部分行业影响较大，如航空交通运输、宾馆餐饮等，但有些行业影响不大，如电信、通讯等，因此总平均来看，影响还不算太大。（4）模型的评价从三方面的结果分析，可以看出模型的结论与实际情况相符，这说明了模型的正确性和可靠性。虽然该模型是就某经济指标的发展规律进行评估预测而建立的，但类似地也适用于其它方面的一些数据规律的评估预测问题，即该模型具有很广泛的应用性。10.2.4GM（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）模型表示一阶的含有两个变量灰色模型，其相应的白化微分方程模型为时间相应函数(即解)为还原值为例10.4某系统中两因素和相互之间存在关系，其中因素为系统特征因素，因素为相关因素，两因素的关系如表10.13所示表10.13两因素的关系序号12345因素因素10.712.315.419.221.540.350.560.365.575.2解：第一步，数据处理取原始数据为对原始数据累加生成，得第二步，建立GM（1，2）模型白化的微分方程模型为转化为灰微分方程或其中。即矩阵形式为其中(3)时间响应函数利用最小二乘法得到参数的估计值，即进而得到微分方程的解即时间响应函数及还原值，即第三步，模型的求解建立M文件10toliti04.m如下：clc,clearx10=[10.7,12.3,15.4,19.2,21.5];x20=[40.3,50.5,60.3,65.5,75.2];n=length(x10);x11=cumsum(x10)x21=cumsum(x20)fori=2:nz11(i)=0.5*(x11(i)+x11(i-1));endB=[-z11(2:n)',x21(2:n)'];Y=x10(2:n)';A=inv(B'*B)*B'*Yx=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-xdelta1=abs(epsilon1./x11)epsilon0=x10-xhatdelta0=abs(epsilon0./x10)利用程序10toliti04求解,得，进而得到时间响应函数及还原值第四步，精度检验先对累加生成数据进行的误差检验，如果进度较高可直接求解还原值，否则要进行模型的修正，