灰色预测系统基于GM(1,1)的matlab程序

灰色系统G(1,1)预测步骤【模板带代码】

=3499.075e -0.1062t

3641.075
编写程序
u=alpha(2)/alpha(1) v=X0(1)-u v=3499.075 u=—3641.075
(5)进行参差检验
1）根据预测公式，计算
v=3499.075 u=—3641.075
Xˆ
1
k
1

X
0
1
for n=0:10
X2(n+1)=v*exp(-alpha(1)*n)+u
end
X2
2.0690
2）累减生成序列
Xˆ X3 =1.0e+003 * (0) 0.1420 0.4079 0.4536
0.5044
0.7713 0.8577 0.9537 1.0605
源程序：X3(1)=X2(1)
for m=1:10
kesi =
4.4388 339.0664 176.2445 203.6132
0 0.1998 1.2682 0.2130 0.8524 0.1330
0.0089
0.3767
0.2203
0.4155
{0%，19.98%，126.82%，0.89%，37.67% ，22.03% ，41.55% ，21.30%，85.24%，13.30%}
e=
179.4592 111.5134 74.1747 175.0204 159.6072 29.2461 215.2168 33.1910
3.2147 24.1540
源程序：S0=0.6745*X0std e=abs(daita0-daita0mean) 对所有的 e 都小于 S0 ，故小参差概率 P(k S0) 1 0.95

灰色预测系统GM(11)模型及其Matlab实现

然后利用生成后的数列进行建模，预测时再情况下，对方于水旱灾害的要求将越来越高。展望增数列，还原成原数列，以恢复事物的原貌。２０１７年以后相当长的时期内，河南防御旱灾面临的
中图分类号：Ｓ１６２．３文献标识码：Ｂ
ＧＭ（１，１）ＭｏｄｅｌａｎｄＭａｔｌａｂＡｐｐｌｉｃａｔｉｏｎｏｆＧｒａｙＰｒｅｄｉｃｔｉｏｎＳｙｓｔｅｍ
文章编号：１００７— ７５９６（２０１７）０７— ００１６一Ｏ３
灰色预测系统ＧＭ（１，１）模型及其Ｍａｔｌａｂ实现
殷鹏远
（辽宁省锦州水文局，辽宁锦州１２１０００）
摘
要：灰色模型有严格的理论基础，最大的优点是实用，用灰色模型预测的结果比较稳定，不
ＹＩＮＰｅｎｇ — ｙｕａｎ
（ＬｉａｏｎｉｎｇＰｒｏｖｉｎｃｉａｌＪｉｎｚｈｏｕＨｙｄｒｏｌｏｇｉｃａｌＢｕｒｅａｕ，Ｊｉｎｚｈｏｕ１２１０００，Ｃｈｉｎａ）
Ｋｅｙｗｏｒｄｓ：ｇｒａｙｓｙｓｔｅｍ；ＧＭ（１，１）ｍｏｄｅｌ；Ｍａｔｌａｂ；ｐｒｅｄｉｃｔｉｏｎｏｆｄｉｓａｓｔｅｒｓ

灰色模型预测GM（1,1）MATLAB程序代码

灰⾊模型预测GM（1,1）MATLAB程序代码版权所有引⽤请注明出处function gmcal=gm1(x)%% ⼆次拟合预测GM(1,1)模型%x = [5999,5903,5848,5700,7884];sizexd2 = size(x,2);%求数组长度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加⽣成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，⽤于计算Byn1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k,yn1sizez1=size(z1,2);%size(yn1);z2 = z1';z3 = ones(1,sizez1)';YN = yn1'; %转置B=[z2 z3];au0=inv(B'*B)*B'*YN;au = au0';afor = au(1);ufor = au(2);ua = au(2)./au(1);constant1 = x(1)-ua;afor1 = -afor;x1t1 = 'x1(t+1)';estr = 'exp';tstr = 't';leftbra = '(';rightbra = ')';strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应⽅程k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(-(k2-1)*afor);endendsizeze1 = size(ze1,2);z4 = ones(1,sizeze1)';G=[ze1' z4];X1 = x1';au20=inv(G'*G)*G'*X1;au2 = au20';Aval = au2(1);Bval = au2(2);strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(Bval),rightbra) %输出时间响应⽅程nfinal = sizexd2-1 + 1; %决定预测的步骤数5 这个步骤可以通过函数传⼊%nfinal = sizexd2 - 1 + 1;%预测的步骤数 1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%⼀次拟合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendendx31fcast%⼀次拟合预测值for k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendendx41fcast,x%⼆次拟合预测值%***精度检验p C************//////////////////////////////////k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexd2elseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexd2elses1total = s1total + (x(k5) - xavg)^2;endends1suqare = s1total ./ sizexd2;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 残差数列x的⽅差 s1sqrt 为x⽅差的平⽅根S1k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexd2elses2total = s2total + (err1(k5) - err1avg)^2;endends2suqare = s2total ./ sizexd2;%s2suqare 残差数列err1的⽅差S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C检验值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )elseendendpval = pnum ./ sizexd2;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围预测步长和数据长度可调整程序参数进⾏改进。

GM(1,1)模型的Matlab实现(范例2)

GM （1，1）模型的MATLAB 实现摘要本文介绍了灰色预测()1,1GM 模型的基本原理，并在其基础上利用MATLAB 软件进行实现，给出相应的MATLAB 算法。

并通过实例检验MATLAB 算法的正确性与通用性。

关键词模型；灰色系统； MATLAB ；1 引言1.1灰色预测理论基本介绍灰色系统是指系统论中一部分信息是已知的，但另一部分信息是未知的，系统内各因素间有不确定的关系。

而灰色预测模型是灰色系统理论的重要内容之一，其中灰色模型是邓聚龙教授在1982年提出的一种新的模型计算方法，被广泛应用于农业、工业、医学、军事、经济、交通、生态等许多科学领域，解决了许多领域中的复杂实践问题[13]-。

且对于杂乱无章表征数据，灰色预测能发现其内在必然的联系与规律，即通过鉴别系统各因素之间发展趋势的相异程度，进行关联分析，并对原始数据进行一定的处理，来寻找系统变动的规律，生成有较强规律性的数据序列，建立相应的微分方程模型，从而预见系统的未来的发展趋势及状况。

灰色预测法是用等时距观测到的反映预测对象特征的一系列数量值，构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

灰色预测有其四种常见的类型，分别为：灰色时间序列预测、畸变预测、系统预测、拓扑预测。

灰色模型（Grey Models ）简称GM 模型，是灰色系统理论（Grey System Theory ）的一个基本模型。

GM 模型的一般模型为，当中时，GM 模型不能做预测，只能用于分析因子之间的相互作用。

做预测用的一般为模型，其中，用的最多的则是用灰色微分拟合法建立的模型。

由于它预测要求所需样本量少（至少四组），不需要计算统计特征量、不用考虑样本变化的趋势、运算简便、在短期内精度高、预测效果好、容易检查等优点。

本文通过运用MATLAB 软件对模型进行程序实现。

增强模型通用性与实用性，方便今后对其使用。

假设我们有一段原始时间序列：为了弱化原始时间序列的随机性，在建立灰色预测模型之前，我们需要先对原始时间序列进行数据处理，经过数据处理后的时间序列称为生成列。

灰色预测模型matlab程序精确版

%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型function gmcal=gm1(x)if nargin==0x=[1019,1088,1324,1408,1601]endformat long gsizex=length(x);%求数组长度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，用于计算Byn1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k,yn1sizez1=length(z1);%size(yn1);z2 = z1';z3 = ones(1,sizez1)';YN = yn1'; %转置%YNB=[z2 z3];au0=inv(B'*B)*B'*YN;au = au0';%B,au0,auafor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,ua%输出预测的 a u 和 u/a的值constant1 = x(1)-ua;afor1 = -afor;x1t1 = 'x1(t+1)';estr = 'exp';tstr = 't';leftbra = '(';rightbra = ')';%constant1,afor1,x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra)%输出时间响应方程%******************************************************%二次拟合k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(-(k2-1)*afor);endend%ze1sizeze1=length(ze1);z4 = ones(1,sizeze1)';G=[ze1' z4];X1 = x1';au20=inv(G'*G)*G'*X1;au2 = au20';%z4,X1,G,au20Aval = au2(1);Bval = au2(2);%Aval,Bval%输出预测的 A,B的值strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra)%输出时间响应方程nfinal = sizex-1 + 1;（其中+1可改为+5等其他数字，即可预测更多的数字）%决定预测的步骤数5 这个步骤可以通过函数传入%nfinal = sizexd2 - 1 + 1;%预测的步骤数 1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fcast%一次拟合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendendx31fcast%一次拟合预测值for k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendendx41fcast,x%二次拟合预测值%***精度检验p C************////////////////////////////////// k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexelseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses1total = s1total + (x(k5) - xavg)^2;endends1suqare = s1total ./ sizex;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 残差数列x的方差 s1sqrt 为x方差的平方根S1 k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses2total = s2total + (err1(k5) - err1avg)^2; endends2suqare = s2total ./ sizex;%s2suqare 残差数列err1的方差S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C检验值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )elseendendpval = pnum ./ sizex;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围预测步长和数据长度可调整程序参数进行改进运行结果x =1019 1088 1324 1408 1601ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929)ans =x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685)x31fcast =Columns 1 through 31019 1122.89347857097 1264.43142178303 Columns 4 through 61423.80987235488 1603.27758207442 1805.36675232556 x41fcast =Columns 1 through 31019 1118.05685435129 1269.65470492098 Columns 4 through 61429.69153740195 1609.90061644041 1812.82460377782 x =1019 1088 1324 1408 1601Cval =0.139501578334155 pval =1。

灰色预测MATLAB程序

灰色预测专设工⑼他QA—叫吋)为原始数列.其1次累❖加生成数列为恥=妙①曲⑵，…卅何),其中X° 仇)二工* ° (0.址=1=2= -:n5-1卷定义卫的决导数为d(k) = *町(上)=x 叫咼-x cl)(Jt-l).令为数列工①的邻值生成数列.即却(去)=^(*) + (1- a)x0)(t-lX于是定义GM (L 1)的灰微分方程模型为d(k)-血⑴住)=K即或严＞(£) + “尹⑻=人⑴在式(1)中』。

＞(灼称为灰导数，我称为发展系数, 弧称为白化背景值，b称为灰作用量乜将时刻表殳二2「3「/代入(1)式有V!1「—ay=代⑶ B =Ib*- :X闵0)-Z，：］(K)1于是G\I <1»1）複至可表示为Y = Bu.現在问题归结为求sb 在值。

用一元线性回归・即最小二秦法求它们的活计值为注二实陌上回归分析中求估计值是用软件计尊的・有标准程序求解，iOmaClab 等。

GM <1» 1>的白化晏対于G\I <1> 1）的灰微分方程（1） >如果将灰导数打（Q 的时刻视为连绫变里"则x°）视为时问（函数卅⑺,于是*〉（Q 対血于导数里级心2 >白化背臬值申的对应于导数卅⑴。

于是G\I （1，1）的坝徽分方樂対应于的白微分方程为内・则数堀列X©可以塗互G\I <19 1）且可以进行页色预测。

否朋，対数摄做适当的克换处理■如平移叢换:取C 使得鞍据列严伙）=工⑴伙）+ G 上=1,2,…,的级比都華住可吝禎盖内。

心⑴⑴ + o?i> (r)二◎ dr<2)GM mi ）质色预测的步骤1 •教摇的枪绘与处連为了ftilGAl （1,1）建複方法的可行性，亲要为已知期S 做必要的检蛉处理。

设原始教据列为了逛=（乂°（1）*6（2）严炉00； >计算数列的级比如果所有的级比都落在可容覆盖区间 • fc =A-2,3"・如果対所有的|p 伙)|<0・1 -则认为达到较高的要求,否则若旳所有的|。

灰色预测MATLAB程序

灰色预测心设尹曲⑴#为原始数列，其1次累＜加生成数列为炉=(孝①宀2\S,其中©=2^°:⑺卫=12…止i-1尋定文沙的灰导数为d(Jt)=玄㈣(Jt)=尤⑴的-工⑴(*-1).令尹为数列壬⑴的邻值生成数列，即尹)(町=加小(町十(1—a)x山(k-1).于是定文GM(1T1)的灰微分方程模型为d(k)+az①(上)=&_即或.严⑹+盘⑴懐)=乩⑴在式(1)中口①的称为灰导数’熬称为发展系数'弧称为白化背景值，b称为灰作用量。

将时刻表庄=23…用代入(O式有j<0)(2)-az⑴(2)=工®⑶—俺叫巧=»于是GMIL)樫型可表示为r=现在问题归结为求巧h在值。

用一元绒性回归，即最小二垂進求它们的估计值住=[]卜护跖护F奕厢上回归分析中求诂计值是用软件计算的，有标淮程博求解，如山訥甜等。

GM(1.1)的白化型对于的(1-1)的获微分方程⑴，如果将解导教矿悶的时報=%…屮观対连续叢里"则工⑴衩为时间i函敕卅®,于是-<'W耐应于导敕重级必％),白化背杲值刃(時对应于导數申⑴。

于是GM(1，1)的换微分方嗨对应于的白微分方程为写®4曲％「)=也⑵GAI(1>1)换色预刪的步叢1-數堀的椅噓弓处理为了保证©M(B1)屋複方达的可行性・需要対已却皴堀锁必要的检峻处Ho 设療皓数攥列为了-计算埶列的级比如果所有的级比都落在可容覆盖区间盂-內・则數摒列X糾可咲建立G*ICL-1)複型且可以避行页色预测。

否则，丙軌据懺适当的叢换处理，如平移銮换：取C使得敕培列严⑹二工蚀盘)+匚用二12…”的级比都落在可啓禎盖内。

(1)残差檢验：计算相对薙差Z 建立GM (L T 1)複型不妬设少弋以m 叫唠霸足上面的要求，以它芮議堀列建立GM(1>1)型蛊(仍(i)+血C1\A)=b ・用回归分祈求得目上的估计值"于是相应的白化模型为气^十小卄工解为工叱)=0)①—勺中1-色-⑶ 应Q于是停到预测值壬⑴(上+1)=0叫1)一勺>加+仝血二12…卫一1=aa伙而相应地得到预«=x co \t +1)=x 0)(t+l)-x a)(i)3i =1,2,-?n-l ?如果对所有的^<0.1・则认为达到鞭嵩的要求：否则，若耐所有的|^)1<0^，则认対达到一般要求©(2)级比偏差値桧验：计算能)=1-呂学©如果对所有的|，则认为达列较高的要求孑吾则若对斫有的，则认为达到一般要求O灰色预测计算实例^…;=:=-■■■■昏例北方某城市1986—1992年道路交通噪声平均声级数据见表6序号年吶寺表拆市近年来交通噪声数据[眶(应)]二諾;二319S872.4第—爭:级比检验建立丢通噪屛均声级数锯时间序列如下：4198972.1j 1990?1.4 619?17201199771.6艸=(•严①卫购(2)厂卅⑺) =(711,72.4.71.4,72.1.71.4,7UQ.71.6)些(1)求级比k(k)忠防护住T)2=(几⑵山⑶.…也⑺)g=(0.982JJ.0042J.0098-0.9917J.0056)(2)级比判断由于所有的X.(10e[0.982J.009S],k=2,3.6故可以用双0)作满意的GM(1,1)建模’第二步：GM(1,1)建模(1)对原始数据X®作一次累加，即卞⑴=(71.L143.5215.9.288359.4.431.4,503)(2)构造数据矩阵B及数据向量Y-2)—H 弋3/>1⑶讦算1T心求解得F'⑴=(工倒〔1〉_-)e 弋Q f+-1*^+1)=0＜l,U)--)£-t +-=-3092^--^+31000⑶求生咸数列值歸型齊看:n令“is 那血由上面的碉醯数可甲得，其中取菱由龙⑴(i}=恥壮曲5加得丁I —"炉閃=进悶-进德-尊(71儿72.4.72.2:72.1:71.9:71.7,71.6)^}=(s"a ＞亍⑴⑵，…，网⑺A＜第三步;模型检验•＞模型的各种检验指标值的计算结果见表工 •t*表7GM(1检验表＜序号年俯原始值模型值残差相对误差级比偏差•>1 19S6 71.1 71.1<219S7 72.4 72.4 -0.0057 0.01%0.0023 <3 19S S 72.4 72.2 0.163S 0.23%0.0203 •>4 19S9 72.1 72.1 0.0329 0.05%-O.(K H8 •>5199071.4 71.9 -0-49S4 0.7%-0.0074 <61991 72.0 71.7 0.21599 037%0.0107<71992 71.6 71.6 0.037S0.05%-0.0032于是得到目=山的餡,立欖型7-B)'1B TY=(dt0.0023 72.6573dt+0.002ix (1>=72.657^心经验证・该模型的精度较高.可进行预测和预报计算的Matlab 程序如下:仃坝测和预报n=length(x); z=0;%取输入数据的样本量for i=1:nz=z+x(i,:)be(i,:)=z; %计算累加值，并将值赋予矩阵beend for i=2:n %对y(i-1,:)=x(i,:)%对原始数列平行移位 endfor i=1：n-1%计算数据矩阵B 的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:)); clCjdearxO=[71H 72.472A 72J71477m c n.lengthtxO)；*'b%注意这里为列帖lamda =xD(l ：n-1),A0(2：n)%计算级比range =minmaxflamda f )%计算级比的范阖 X1=cumsum(xO)；%累加运算B=['0,5*(xl(l ；n ^l)+xl(2:n))t ones(n -1,1)]TY 二甸(2：町；口=B\Y%拟合参数u(l>=a .u(2)=bx=dsolve (+a 'x =b\f x(0)-xO^J ；%求徴分方程的特号解x =subs(xJ*a\,b r /xO ，Mu(l)P u(2)t xO(l)|)i%代入荷计痹擞值和初蜡值yucel =subs %求巳知数擁的扳测位y-vpa(x,6)奄其中的石表示显不白位数字yuce=[x0(l)T diff(yucel)]%羔分运算，还原数据 epsiIon=-yuce%计算战羞作用：求累加数列、求ab 的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.172.472.472.171.472.071.6]; format long ;%设置计算精度if length(x(:,1))==1%对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x endM.I-JTVorhlllst 模型endfor j=1:n-1%计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1%构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y;%计算参数矩阵即ab的值for i=1:n+1%计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha( 2,:)/alpha(1,:);%显示输出预测值的累加数列endvar(1,:)=ago(1,: )for i=1:n%显示输出预测值%如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:);%估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:);%计算残差endc=std(error)/std(x);%调用统计工具箱的标准差函数计算后验差的比值cago alpha var%显示输出预测值的累加数列%显示输出参数数列%显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求ab的值、求预测方程clc,clearx0=[71.172.472.472.171.472.071.6]';%注意这里为列向量n=length(x0);lamda=x0(1:n-1)./x0(2:n)%计算级比range=minmax(lamda')%计算级比的范围x1=cumsum(x0)%累加运算B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];Y=x0(2:n);u=B\Y%拟合参数u(1)=a,u(2)=bx=dsolve('Dx+a*x=b','x(0)=x0');%求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)})%代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]);%求已知数据的预测值y=vpa(x,6)%其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)]%差分运算，还原数据。

灰色系统预测GM(1-1)模型及其Matlab实现

灰色系统预测GM(1,1)模型及其Matlab 实现预备知识（1）灰色系统白色系统是指系统内部特征是完全已知的;黑色系统是指系统内部信息完全未知的;而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知,另一部分信息未知或不确定。

(2)灰色预测灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。

尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。

灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。

目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的ＧＭ（1，1)模型。

它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。

经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。

因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。

1 灰色系统的模型GM(1,１）１．1 GM(1，1)的一般形式设有变量X (0)={X (0)(i ）,i ＝1,2,．..，n}为某一预测对象的非负单调原始数据列，为建立灰色预测模型:首先对X（0)进行一次累加（１—ＡGO ， Acum ｕl ａt ｅｄ Ge ｎeｒa ｔｉng Opera ｔo ｒ)生成一次累加序列: X (1)={X（1）(k ),k =1，2,…,n}其中X (1）(k )＝∑=ki 1X (0）(i)＝X (1)(ｋ-１）+ X （0）(ｋ) (1)对X（1)可建立下述白化形式的微分方程：dtdX )1(十)1(aX =u (2)即G Ｍ（１，1)模型。

上述白化微分方程的解为(离散响应): ∧X（1)(k +１)＝(X (０)(1)－a u ）ake -+au (3）或∧X （1）（k )=(X (0）(1)-a u ))1(--k a e +au （４)式中：k 为时间序列，可取年、季或月。

Matlab灰色预测模型GM(1,1)代码

end
s=0;xx0(1)=X0(1);
for jj=2:1:m2;
xx0(jj)=xx1(jj)-xx1(jj-1);
end
disp('GM(1,1)对数列进行预测结果');xx0
disp('数列1原始观测数据');X0
disp('a');AU(1)
disp('u');AU(2)
z=z+x(i,:);
be(i,:)=z;
end
for i=2:n %对原始数列平行移位
y(i-1,:)=x(i,:);
end
for i=1:n-1 %计算数据矩阵B的第一列数据
c(i,:)=-0.5*(be(i,:)+be(i+1,:));
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
yy(i)=yy(i-1)+y(i);
disp(['再下个拟合值为',num2str(ys(n+2))]);
x(0)={x(0)(1),x(0)(2),……,x(0)(N)}=
{ 724.57, 746.62, 778.27, 800.8, 827.75,871.1, 912.37, 954.28, 995.01, 1037.2}
function c7fun73

MATLAB实现灰色预测程序

MATLAB实现灰色预测程序function [y,p,e]=huise_1_1(X,k) %灰色模型的malab程序%Example [y,p]=gm_1_1([200 250 300 350],2)%接口描述：X的预测的初始数列，|X|>4，K是指向后进行预测的个数%命令格式：程序保存的文件名，eg：huise.m 则命令是：huise([579.8 547.5 527.0 492.3 437.0],5) if nargout>3;error('Too many output argument.');endif nargin==1,k=1;x_orig=X;elseif nargin==0|nargin>2error('Wrong number of input arguments.');endx_orig=X;predict=k; %AGO 处理，即是对初始数列进行一阶累加x=cumsum(x_orig); %计算系数(a 和u)------------------------n=length(x_orig); %生成矩阵Bfor i=1:(n-1);B(i)=-(x(i)+x(i+1))/2;endB=[B' ones(n-1,1)]; %生成矩阵Yfor i=1:(n-1);y(i)=x_orig(i+1);endY=y'; %计算系数a=au(1) u=au(2)au=(inv(B'*B))*(B'*Y); %--------------------------------------------------------%把huise模型公式转换成符号coef1=au(2)/au(1);coef2=x_orig(1)-coef1;coef3=0-au(1);costr1=num2str(coef1);costr2=num2str(abs(coef2));costr3=num2str(coef3);eq=strcat(costr1,'+',costr2,'e^',costr3,'*(t-1))'); %计算每一个值for t=1:(n+predict)mcv(t)=coef1+coef2*exp(coef3*(t-1));endx_mcv0=diff(mcv);x_mcve=[x_orig(1) x_mcv0] %输出图形中的各点x_mcv=diff(mcv(1:end-predict));x_orig_n=x_orig(2:end);x_c_error=x_orig_n-x_mcv;x_error=mean(abs(x_c_error./x_orig_n));if x_error>0.2 %相对误差的均值disp('model disqualification!');elseif x_error>0.1disp('model check out');elsedisp('model is perfect!');endplot(1:n,x_orig,'o',1:n+predict,x_mcve);p=x_mcve(end-predict+1:end); %画出预测模型和初始数列的点xlabel('年份(从第一个数据年份起)');ylabel('产水量（万吨）');title('灰度模型GM(1,1)');grid ony=eq;e=x_error;p=x_mcve(end-predict+1:end);。

灰色模型GM(1,1)&matlab实现

单位：亿请根据1978-2012年的全国人口数据，用GM(1,1)来预测未来人口数量，能预测50年内的数据吗？100年呢？为什么？摘要：关键字：人口预测、GM(1,1)一、问题的重述二、问题的分析三、模型假设1、假设本问题所使用的数据均真实有效，具有统计分析价值2、不考虑战争、灾害、疾病对人口数目的影响；3、...4、...四、符号说明与名词解释4.1符号说明见模型假设。

4.2名词解释...五、模型的建立与求解5.1模型建立在灰色系统理论中，称抽象的逆过程为灰色模型，也称GM。

它是根据关联度、生成数灰导数，灰微分等观点和一系列数学方法建立起来的连续型的微分方程。

通常GM表示为GM(n,h)。

当n=h=1时即构成了单变量一阶灰色预测模型。

设原始时间序列为：(0)(0)(0)(0)(0)[(1),(2),(3),,()]X x x x x n=设(1)X为(0)X的一次累加序列：即：(1)(0)(1)(1)(1)(1)(0)()(-1)(),2,, x xx k x k x k k n ⎧=⎪⎨⎪=+=⎩得(1)(1)(1)(1)(1)[(1),(2),(3),,()] X x x x x n=利用(1)X 计算GM(1,1)模型参数a 、u ，令ˆ[,]T aa u =, 则有1ˆ()T T N aB B B Y -= 其中(1)(1)(1)(1)(1)(1)1((1)(2))121((2)(3))121((1)())12x x x x B x n x n ⎡⎤-+⎢⎥⎢⎥⎢⎥-+⎢⎥=⎢⎥⎢⎥⎢⎥--+⎢⎥⎣⎦(0)(0)(0)[(2),(3),,()]T N Y x x x n =由此获得：(1)(1)ˆ(1)((1))ak u uxk x e a a -+=-+ 于是：(0)(1)(1)ˆˆˆ(1)(1)()1xk x k x k +=+-（）或(0)(0)ˆ(1)-((1))ak uxk a x e a-+=-（2）注：（1）式是根据(0)X 和(1)X 的关系的到的（2）式是利用数学求导还原得到的至于用哪个，最好看相对误差5.2模型求解：求解得：六、模型检验6.1残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验。

改进的灰色预测GM(1,1)模型的MATLAB实现

一
６９一
２０１４年４月第７期
江苏科技信息・网络与通信
Ｎｏ．７
Ａｐｒｉｌ，２０１４
＾
＾
＾
（ｏ（＋１）：（１（＋１）一（（）＝Ｃ・ｅ－ａｋｋ＝ｌ，２， …，一１（６）
Ｘ（（ｋ＋１）＝ｃｅ－ａｋ＋，ｋ＝０，１， …，假定一个初始值。假定ｘ（（１）；０（１），
则有
图１模型改进流程图
作者简介：杨旭（１９９２一），男，河南新野，本科生；研究方向：水污染控制理论与技术。
数据处理问题提供了方便。
关键词：灰色预测；ＧＭ（１，１）模型；改进模型；ＭＡＴＩＡＢ算法程序
Ｏ引言
龛㈩（１）：ｃ＋：（０（１）
ｃ＝Ｘ（０（１）一
ｎ
灰色系统理论 … 是由我国学者邓聚龙教授于１９８２年在国际上首先提出来的，用于研究少数据、贫信息的不确定性问题的理论方法。该理论的主要内容之一就是以ＧＭ（１，１）模型为核心的预测模型体系。该模型在工业、农业、商业等经济领域以及环境、社会等领域中都有广泛应用。然而在使用ＧＭ（１，１）模型进行预测的过程中，也会出现预测模型精度较低的情况。许多学者提出了改进预测模型精度的方法［２－３］。其中，杨华龙４等学者在分析了以往学者的改进方法后认为虽然以往学者提出的模型改进方法对模型精度的提高有所帮助，但模型预测公式本身存在的缺陷并未得到有效的改进。因此在分析了ＧＭ（１，１）模型预测公式的形成过程后，提出并使用自动寻优定权对背景值进行了选择，使用最小二乘法原理对ＧＭ（１，１）模型的初始值进行了改进。且通过实例结果表明，提出的改进方法是有效和完

数学建模-灰色预测模型GM（1,1）_MATLAB

数学建模-灰⾊预测模型GM（1,1）_MATLAB %GM(1,1).m%建⽴符号变量a(发展系数)和b(灰作⽤量)syms a b;c = [a b]';%原始数列 AA = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];%填⼊已有的数据列！n = length(A);%对原始数列 A 做累加得到数列 BB = cumsum(A);%对数列 B 做紧邻均值⽣成for i = 2:nC(i) = (B(i) + B(i - 1))/2;endC(1) = [];%构造数据矩阵B = [-C;ones(1,n-1)];Y = A; Y(1) = []; Y = Y';%使⽤最⼩⼆乘法计算参数 a(发展系数)和b(灰作⽤量)c = inv(B*B')*B*Y;c = c';a = c(1);b = c(2);%预测后续数据F = []; F(1) = A(1);for i = 2:(n+10) %这⾥10代表向后预测的数⽬，如果只预测⼀个的话为1F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;end%对数列 F 累减还原,得到预测出的数据G = []; G(1) = A(1);for i = 2:(n+10) %10同上G(i) = F(i) - F(i-1); %得到预测出来的数据enddisp('预测数据为：');G%模型检验H = G(1:10); %这⾥的10是已有数据的个数%计算残差序列epsilon = A - H;%法⼀：相对残差Q检验%计算相对误差序列delta = abs(epsilon./A);%计算相对误差Qdisp('相对残差Q检验：')Q = mean(delta)%法⼆：⽅差⽐C检验disp('⽅差⽐C检验：')C = std(epsilon, 1)/std(A, 1)%法三：⼩误差概率P检验S1 = std(A, 1);tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);disp('⼩误差概率P检验：')P = length(tmp)/n%绘制曲线图t1 = 1995:2004;%⽤⾃⼰的，如1 2 3 4 5...t2 = 1995:2014;%⽤⾃⼰的，如1 2 3 4 5... plot(t1, A,'ro'); hold on;plot(t2, G, 'g-');xlabel('年份'); ylabel('污⽔量/亿吨');legend('实际污⽔排放量','预测污⽔排放量'); title('长江污⽔排放量增长曲线'); %都⽤⾃⼰的grid on;。

灰色预测的matlab实现,代码大全

disp('确定的模型为:x1(k+1)=(x(1)-u_a)*exp(-k*a)+u_a') 改进灰色模型代码： for i=1:inf if i==1
e(i)=input('输入数据:'); else if e(i-1)==0 break; else end
e(i)=input('输入数据:');
6
1
x2(k+1)=(x0(1)-u_a)*exp(-k*a)+u_a; end %求出新数列的模拟值 x2 x3=x2; x3(m2+k0)=[]; x4=[0 x3]; x5=x2-x4;%利用累减生成法求出原始数据的模拟值 x5 disp('3.一次累加（1-AGO）生成的数据的模拟值:') Y=''; for z=1:m2+k0 Y=strcat(Y,'(',num2str(x2(z)),')'); end disp(Y); disp('4.原始数据的模拟值:') Y=''; for z=1:m2+k0 Y=strcat(Y,'(',num2str(x5(z)),')'); end disp(Y); % 4. 模型检验（算出的值到等级参照表中检查其精度等级） %计算后验差比 C 和残差序列 Q x6=x5(1:m2); Q=x0-x6;%Q 为残差序列 s1=std(Q);%s1 为残差序列 Q 的标准差 s2=std(x0);%s2 为初始序列 x0 的标准差 C=s1/s2;%后验差比 C(越小越好) w1=1:m2; w1=[ones(m2,1) w1']; w2=Q'; [bb,bint,r1,rint,stats]=regress(w2,w1); rcoplot(r1,rint) C1=strcat('5.后验差比(均方差比值): C=',num2str(C));disp(C1); if C<=0.35 disp(' 由于 C<=0.35，则此模型精度等级为 1 级（好）。'); else if C<=0.5 disp(' 由于 0.35<C<=0.5，则此模型精度等级为 2 级（合格）。'); else if C<=0.65 disp(' 由于 0.5<C<=0.65，则此模型精度等级为 3 级（勉强）。'); else disp(' 由于 C>0.65，则此模型精度等级为 4 级（不合格）。'); end end end

灰色预测系统基于GM(1-1)的matlab程序

function GM1_1(X0)%format long ;X0=input('请输入实测数据');%实测值[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] ; % 求B矩阵YN=X0(2:end) ;Pt=YN./X1(1:(length(X0)-1)) %对原始数据序列X0进行准光滑性检验， %序列X0的光滑比P(t)=X0(t)/X1(t-1)A=inv(B.'*B)*B.'*YN.' ;a=A(1)u=A(2)c=u/a ;b=X0(1)-c ;X=[num2str(b),'exp','(',num2str(-a),'k',')',num2str(c)];strcat('X(k+1)=',X)%syms k;for t=1:length(X0)k(1,t)=t-1;endkY_k_1=b*exp(-a*k)+c;for j=1:length(k)-1Y(1,j)=Y_k_1(j+1)-Y_k_1(j);endXY=[Y_k_1(1),Y] %预测值CA=abs(XY-X0) ; %残差数列Theta=CA %残差检验绝对误差序列XD_Theta= CA ./ X0 %残差检验相对误差序列AV=mean(CA); % 残差数列平均值R_k=(min(Theta)+0.5*max(Theta))./(Theta+0.5*max(Theta)) ;% P=0.5R=sum(R_k)/length(R_k) %关联度Temp0=(CA-AV).^2 ;Temp1=sum(Temp0)/length(CA);S2=sqrt(Temp1) ; %绝对误差序列的标准差%----------AV_0=mean(X0); % 原始序列平均值Temp_0=(X0-AV_0).^2 ;Temp_1=sum(Temp_0)/length(CA);S1=sqrt(Temp_1) ; %原始序列的标准差TempC=S2/S1*100; %方差比C=strcat(num2str(TempC),'%') %后验差检验 %方差比%----------SS=0.675*S1 ;Delta=abs(CA-AV) ;TempN=find(Delta<=SS);N1=length(TempN);N2=length(CA);TempP=N1/N2*100;P=strcat(num2str(TempP),'%') %后验差检验 %计算小误差概率m=input('请输入预测期数:');for g=1:(length(X0)+m)v(1,g)=g-1;endvm=b*exp(-a*v)+c;for j=1:length(v)-1l(1,j)=m(j+1)-m(j);endxyz=[m(1),l]%预测值disp(['小误差概率为：',num2str(P)]);disp(['后验方差比为：',num2str(C)]);disp(['预测值为：',num2str(xyz)]);。

灰色预测模型GM(1,1)的matlab运行代码

灰色预测模型GM(1,1)的matlab 运行代码例由1990—2001年中国蔬菜产量，建立模型预测2002年中国蔬菜产量，并对预测结果作检验。

分析建模：给定原始时间1990—2001年资料序列X )0((k),对X )0((k)生成1-AGO(累加)序列X )1((k)及Y n 。

见下表X )1((k)=)(1i )0(i X k∑=; Y n =T X X X )]12(,),3(),2([)0()0()0( 对上述X )0((k)的GM(1,1),得到[][][][][][][][][][][]⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡+-+-+-+-+-+-+-+-+-+-+-==⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-----------=15236.2-1331173.5-1244348.01204848.5-1168369.5-1135943.5-1107892.5-186730.0-168581.5-148915.5-129308.0-112115.0111105.011095.01985.01875.01765.01655.01545.01435.01325.01215.01)12(1)11(1)10(1)9(1)8(1)7(1)6(1)5(1)4(1)3(1)2()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1(）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（X X X X X X X X X X X X X X X X X X X X X X z z z z z z z z z z z B 将 B 和Y n 代入辨识算式，有：10.1062105()13999.9T Tn a B B B Y b α--⎡⎤⎡⎤==•=⎢⎥⎢⎥⎣⎦⎣⎦得灰色GM(1,1)模型为（1）灰微分方程X )0((k)-0.1062105 Z )1((k)=13999.9(2)白化方程9.139991062105.0)1()1(=-X dtdX (3)白化方程的时间响应式5.1318135.151332)1()1(ˆ1062105.0)0()1(-=+⎥⎦⎤⎢⎣⎡-=+--t at e a b e a b X t X (4)还原为原始数据预测方程：)(ˆ)1(ˆ)1(ˆ)1()1()0(t X t X t X-+=+,即 t e t X1062105.0)0(15248.968)1(ˆ=+ (5)残差检验：残差error1=e1=)()ˆ)0()0(i X i X -（,这里残差有12个。

灰色模型GM11matlab实现

单位：亿请根据1978-2012年的全国人口数据，用GM(1,1)来预测未来人口数量，能预测50年内的数据吗？100年呢？为什么？摘要：关键字：人口预测、GM(1,1)一、问题的重述二、问题的分析三、模型假设1、假设本问题所使用的数据均真实有效，具有统计分析价值2、不考虑战争、灾害、疾病对人口数目的影响；3、...4、...四、符号说明与名词解释4.1符号说明见模型假设。

4.2名词解释...五、模型的建立与求解5.1模型建立在灰色系统理论中，称抽象的逆过程为灰色模型，也称GM。

它是根据关联度、生成数灰导数，灰微分等观点和一系列数学方法建立起来的连续型的微分方程。

通常GM表示为GM(n,h)。

当n=h=1时即构成了单变量一阶灰色预测模型。

设原始时间序列为：(0)(0)(0)(0)(0)[(1),(2),(3),,()]X x x x x n=L设(1)X为(0)X的一次累加序列：即：(1)(0)(1)(1)(1)(1)(0)()(-1)(),2,, x xx k x k x k k n ⎧=⎪⎨⎪=+=⎩L得(1)(1)(1)(1)(1)[(1),(2),(3),,()] X x x x x n=L利用(1)X 计算GM(1,1)模型参数a 、u ，令ˆ[,]T aa u =, 则有1ˆ()T T N aB B B Y -= 其中(1)(1)(1)(1)(1)(1)1((1)(2))121((2)(3))121((1)())12x x x x B x n x n ⎡⎤-+⎢⎥⎢⎥⎢⎥-+⎢⎥=⎢⎥⎢⎥⎢⎥--+⎢⎥⎣⎦L L (0)(0)(0)[(2),(3),,()]T N Y x x x n =L由此获得：(1)(1)ˆ(1)((1))ak u uxk x e a a -+=-+ 于是：(0)(1)(1)ˆˆˆ(1)(1)()1xk x k x k +=+-（）或(0)(0)ˆ(1)-((1))ak uxk a x e a-+=-（2）注：（1）式是根据(0)X 和(1)X 的关系的到的（2）式是利用数学求导还原得到的至于用哪个，最好看相对误差5.2模型求解：求解得：六、模型检验6.1残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验。

MATLAB实现灰色预测程序

MATLAB实现灰色预测程序function [y,p,e]=huise_1_1(X,k) %灰色模型的malab程序%Example [y,p]=gm_1_1([200 250 300 350],2)%接口描述：X的预测的初始数列，|X|>4，K是指向后进行预测的个数%命令格式：程序保存的文件名，eg：huise.m 则命令是：huise([579.8 547.5 527.0 492.3 437.0],5) if nargout>3;error('Too many output argument.');endif nargin==1,k=1;x_orig=X;elseif nargin==0|nargin>2error('Wrong number of input arguments.');endx_orig=X;predict=k; %AGO 处理，即是对初始数列进行一阶累加x=cumsum(x_orig); %计算系数(a 和u)------------------------n=length(x_orig); %生成矩阵Bfor i=1:(n-1);B(i)=-(x(i)+x(i+1))/2;endB=[B' ones(n-1,1)]; %生成矩阵Yfor i=1:(n-1);y(i)=x_orig(i+1);endY=y'; %计算系数a=au(1) u=au(2)au=(inv(B'*B))*(B'*Y); %--------------------------------------------------------%把huise模型公式转换成符号coef1=au(2)/au(1);coef2=x_orig(1)-coef1;coef3=0-au(1);costr1=num2str(coef1);costr2=num2str(abs(coef2));costr3=num2str(coef3);eq=strcat(costr1,'+',costr2,'e^',costr3,'*(t-1))'); %计算每一个值for t=1:(n+predict)mcv(t)=coef1+coef2*exp(coef3*(t-1));endx_mcv0=diff(mcv);x_mcve=[x_orig(1) x_mcv0] %输出图形中的各点x_mcv=diff(mcv(1:end-predict));x_orig_n=x_orig(2:end);x_c_error=x_orig_n-x_mcv;x_error=mean(abs(x_c_error./x_orig_n));if x_error>0.2 %相对误差的均值disp('model disqualification!');elseif x_error>0.1disp('model check out');elsedisp('model is perfect!');endplot(1:n,x_orig,'o',1:n+predict,x_mcve);p=x_mcve(end-predict+1:end); %画出预测模型和初始数列的点xlabel('年份(从第一个数据年份起)');ylabel('产水量（万吨）');title('灰度模型GM(1,1)');grid ony=eq;e=x_error;p=x_mcve(end-predict+1:end);。