灰色系统预测GM(1,1)模型及其Matlab实现三天三夜72小时：读懂题目-》查找文献资料-》选择题目-》重查找文献资料-》精读其中几篇-》查找资料的资料。
在数学建模中常常会遇到数据的预测问题，有些赛题中，预测占主导地位，例如：2003年A题 SARS的传播问题；2005年A题长江水质的评价和预测问题；2006年B题艾滋病疗法的评价及疗效的预测问题；2007年A题中国人口增长预测问题。
有些问题则是需要在求解的过程中进行预测，如2009年D题“会议筹备”对与会人数的确定等。
参考资料：《灰色系统理论及其应用第五版》作者：刘思峰，党耀国等著出版时间：2010.05 校超星数字图书馆可阅读。
灰色模型（Gray Model）有严格的理论基础，最大优点是实用。
用灰色模型预测的结果比较稳定，不仅适用于大数据量的预测，在数据量较少时(>3)预测结果依然较准确。
预备知识(1)灰色系统白色系统是指系统内部特征是完全已知的,即人们不仅知道该系统的输入——输出关系，而且知道实现输入——输出关系的结构与过程；黑色系统是指系统内部信息完全未知的,即人们只知道该系统输入——输出关系，但不知道实现输入——输出关系的结构与过程；而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知，另一部分信息未知或不确定。
例如，一个加有电压的电阻，也是一个系统，根据欧姆定律，I=U/R,当电阻的大小知道后，便可由多大电压算出能得到多大电流。
电压与电流之间有明确的关系或函数，这便是白色系统。
因此，这样的系统要求有明确的作用原理，一个有明确作用原理的系统必定是具有确定结构的，必定是有物理原型的。
然而许多社会经济系统都没有物理原型，虽然知道影响系统的某些因素，但很难明确全部因素，更不可能确定因素之间的映射关系。
这种没有确定的映射关系(函数关系)的系统是灰色系统。
(2)灰色预测灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测，也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。
尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。
灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。
目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。
它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。
经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。
因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。
目前，灰色模型GM(1,1)己广泛应用于工程技术、社会、经济、农业、生态、环境等各种系统的预测中。
1 灰色预测基础知识灰色系统理论认为:系统的行为现象尽管是朦胧的,数据是复杂的,但它毕竟是有序的,是有整体功能的。
在建立灰色预测模型之前,需先对原始时间序列进行数据处理,经过数据预处理后的数据序列称为生成列。
对原始数据进行预处理,不是寻求它的统计规律和概率分布,而是将杂乱无章的原始数据列通过一定的方法处理,变成有规律的时间序列数据,即以数找数的规律,再建立动态模型。
灰色系统常用的数据处理方式有累加和累减两种,通常用累加方法。
灰色预测通过鉴别系统因素之间发展趋势的相异程度,并对原始数据进行生成处理来寻找系统变动的规律,生成有较强规律性的数据序列,然后建立相应的微分方程模型,从而预测事物的未来发展趋势。
灰色预测的数据是通过生成数据的模型所得到的预测值的逆处理结果。
灰色预测是以灰色模型为基础的,在诸多的灰色模型中,以灰色系统中单序列一阶线性微分方程模型 GM(1,1)模型最为常用。
下面简要地介绍 GM(1,1)模型。
设有原始数据列x(0)=(x(0)(1),x(0)(2),…,x(0)(n)),n 为数据个数。
如果根据x(0)数据列建立 GM(1,1)来实现预测功能,则基本步骤如下:(1)原始数据累加以便弱化随机序列的波动性和随机性,得到新数据序列:x (1) = {x (1)(1),x (1)(2),…,x (1)(n)}其中,x (1)(t)中各数据表示对应前几项数据的累加。
x (1)(t)= ∑=t1k )0()k (x ， t=1,2,…,n(2)对x (1)(t)建立下述一阶线性微分方程: 即GM(1,1)模型。
u ax dtdx )1()1(=+ 其中,a,u 为待定系数,分别称为发展系数和灰色作用量,a 的有效区间是 (-2,2),并记a,u 构成的矩阵为灰参数∧a =⎪⎪⎭⎫⎝⎛u a 。
只要求出参数a,u,就能求出x (1)(t),进而求出x (0)的未来预测值。
(3)对累加生成数据做均值生成B 与常数项向量Yn,即B ＝⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡+++- 1 (n))X 1)-(n (X 21 ... 1 (3))X (2)X (211 (2))X (1)X (21(1)1(1)(1)(1)(1)）（－－Yn=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛)n (x )3(x )2(x )0()0()0((4)用最小二乘法求解灰参数∧a ,则∧a ＝(B T B)-1B T Y n(5)将灰参数∧a 代入u ax dt dx )1()1(=+ ,并对u ax dtdx )1()1(=+进行求解,得 ∧X (1)(t+1)＝(X (0)(1)－a u )at e -＋au由于∧a 是通过最小二乘法求出的近似值,所 以xˆ(1)(t+1)是一个 近似表达式,为了与原序列x (1)(t+1)区分开来,故记为xˆ(1)(t+1)。
式中t 为时间序列，可取年、季或月。
(6)对函数表达式xˆ(1)(t+1)及x ˆ(1)(t)进行离散,并将二者做差以便还原x (0)原序列,得到近似数据序列xˆ(0)(t+1)如下: xˆ(0)(t+1)=x ˆ(1)(t+1)-x ˆ(1)(t) (7)对建立的灰色模型进行检验,步骤如下:计算x (0)(t)与xˆ(0)(t)之间的残差e (0)(t)和相对误差q (0)(t): e (0)(t)=x (0)(t)-xˆ(0)(t) q (0)(t)=e (0)(t)/x (0)(t)等等。
(8)利用模型进行预测:x ˆ(0) ={x ˆ(0)(1),x ˆ(0)(2),…,x ˆ(0)(n ），xˆ(0)(n+1),…,x ˆ(0)(n+m)} 原数列的模拟 未来数列的预测应用举例取某高校2000年～2005年的某专业招生数据建模,见表1。
表以表1中的数据构造原始数据列X (0),即X (0)＝{X (0)(1),X (0)(2),X (0)(3),X (0)(4),X (0)(5),X (0)(6)}＝{132,92,118,130,187,207}对X(0)进行一次累加(1—AGO)，生成数列：X (1)(k)＝∑=ki 1X (0)(i)即X (1)＝{X (1)(1),X (1)(2),X (1)(3),X (1)(4),X (1)(5),X (1)(6)}＝{132,224,342,472,659,866} 和数据阵B 、数据列YnB ＝⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡----- 1 5.762 1 5.5651 4071 2831178，Yn ＝(92,118,130,187,207)T经计算可得∧a ＝[a,u]T＝⎥⎦⎤⎢⎣⎡-7878184.56205.0 进一步得到灰色预测模型GM(1,1)为∧X (1)(k)＝(X (0)(1)－a u ))1(--k a e ＋a u＝(132＋277.))1(205.0-k e －277.＝409.)1(205.0-k e－277.预测值及预测精度见表2。
表2 某高校专业招生预测值及预测精度表由表2知预测精度较高。
2006年某专业招生人数预测值为259人。
由于人数为整数，所以结果取整数部分。
GM(1,1)也是一种长期预测模型，在没有大的市场波动及政策性变化的前提下，该预测值应是可信的。
众所周知，影响招生人数的因素很多且难以预测。
因此，在采用灰色系统理论进行定量预测时，如果存在对预测对象影响较大的因素，就要在定性分析的基础上，寻找原始数据信息的突变点的量化值，然后再对预测值进行必要的修正，使预测值更接近实际情况，提高预测值的可信度，为科学决策提供可靠的数据。
另外，若作长期预测，要考虑对上限值的约束条件。
2 灰色预测的 MATLAB程序2.1 典型程序结构灰色预测中有很多关于矩阵的运算,这可是 MATLAB的特长,所以用 MATLAB是实现灰色预测过程的首选。
用 MATLAB编写灰色预测程序时,可以完全按照预测模型的求解步骤,即(1)对原始数据进行累加。
(2)构造累加矩阵B 与常数向量Yn。
(3)求解灰参数。
(4)将参数带入预测模型进行数据预测。
下面以某公司收入预测问题为例介绍灰色预测的 MATLAB实现过程。
已知某公司1999—2008年的利润为 (单位:元/年): [89677,99215,109655,120333,135823,159878,182321,209407,2466 19,300670],现在要预测该公司未来几年的利润情况。
具体的 MATLAB程序如下：clearsyms a u;c=[a u]'; %灰参数cA=[89677,99215,109655,120333,135823,159878,182321,209407,24 6619,300670]; %原始数据Ago=cumsum(A); %原始数据一次累加n=length(A); %原始数据个数for i=1:(n-1)C(i)=(Ago(i)+Ago(i+1))/2; %生成累加矩阵end% 计算待定参数的值Yn=A;Yn(1)=[];Yn=Yn';E=[-C;ones(1,n-1)];c=inv(E*E')*E*Yn;c=c';a=c(1);u=c(2);% 预测后续数据F=[];F(1)=A(1);for i=2:(n+10)F(i)=(A(1)-u/a)/exp(a*(i-1))+u/a ;endG=[];G(1)=A(1);for i=2:(n+10)G(i)=F(i)-F(i-1); %得到预测出来的数据endt1=1999:2008;t2=1999:2018;Gplot(t1,A,'o',t2,G) %原始数据与预测数据的比较运行该程序,得到的预测数据如下:G =1.0e+006 * Columns1through140.0897 0.0893 0.1034 0.1196 0.1385 0.1602 0.1854 0.2146 0.2483 0.2873 0.3325 0.3847 0.4452 0.5152 Columns15through200.5962 0.6899 0.7984 0.9239 1.0691 1.2371 该程序还显示了预测数据与原始数据的比较图。