灰色系统gm(1,1)的使用

离落之猪

已于 2024-06-02 16:03:54 修改

阅读量504

点赞数 8

文章标签： matlab

于 2024-05-23 12:30:10 首次发布

本文链接：https://blog.csdn.net/weixin_62763029/article/details/139142256

版权

灰色系统简介

灰色系统：系统内部特性部分已知。

1982 年由 中国科学家邓聚龙 教授提出

特点：

1，基于“小样本、贫信息”不确定性数据挖掘有价值的信息

2 ，从系统内部结构及参数去研究系统

3，强调现实规律的认知

原理：利用离散随机数经过生成变为随机性被显著削弱而且较为

规律的生成数，建立起的微分方程形式的模型

GM(1,1)

已知参考数列： $x0=[X0(1),X0(2),...,X0(n)]$

累加生成数列（1-AGO):

$x1=[X1(1),X1(2),...,X1(n)]=[X1(1),X1(1)+X0(2),...,X(1)(n-1)+X0(n)]$

其中， $X1(k)=\sum_{1}^{k}X0(i)$

x0	1	2	3	4	5	6
x1	1	3	6	10	15	16

均值数列

$Z1(k)=0.5X1(k)+0.5X1(k-1),k=2,3,...,n$

背景值

$z1=[Z1(2),Z1(3),...,Z1(n)]$

灰色微分方程

$d(k)=X0(k)=X1(k)-X1(k-1)$

$X0(k)+aZ1(k)=b,k=2,3,...,n$

a:发展系数

b:灰作用量

最小二乘法求解灰微分方程

$X0(k)+aZ1(k)=b$

$Y=\begin{bmatrix} X0(2)\\ X0(3)\\ ...\\ X0(n) \end{bmatrix}$

$B=\begin{bmatrix} -Z1(2) &1 \\ -Z1(3) &1 \\ ...&... \\ -Z1(n)&1 \end{bmatrix}$

$u=(a,b)^T$

$u=(a,b)^T=(B^TB)^{-1}B^TY$

由此可以解的a,b的值

白微分方程

$dx1/dt+aX1(t)=b$

$X1(k+1)=(X0(1)-b/a)e^{-ak}+b/a,k=1,2,...,n-1$

GM(1,1)预测的步骤

已知参考数列： $x0=[X0(1),X0(2),...,X0(n)]$

计算数列级比： $\lambda (k)=X0(k-1)/X0(k),k=2,3,...,n$

要求级比都可覆盖在 $(e^{-2/n+1},e^{2/n+1})$ 内

否则，需要对x(0)变换处理

$Y0(k)=X0(k)+c,k=1,2,...,n$

$y0=[Y0(1),Y0(2),...,Y0(n)]$

$\lambda (k)=Y0(k-1)/Y0(k),k=2,3,...,n$

建模

预测值

$X^{'}1(k+1)=(X0(1)-b/a)e^{-ak}+b/a,k=1,2,...,n-1$

$X^{'}0(k+1)=X^{'}1(k+1)-X^{'}1(k)$

检测预测值

残差检验： $\varepsilon (k)=\frac{X0(k)-X^{'}0(k)}{X0(k)},k=1,2,...,n$

$\varepsilon (k)<0.2$ ，达到一般要求

$\varepsilon (k)<0.1$ ，达到较高要求

运用matlab实现Gm(1,1)模型

例子

先将数据存入变量中

运行代码即可实现

由此可见，预测精度较差

可选取最近五年的数据预测

代码如下

%年份与车数量的坐标图
x=2021:-1:2017;
i=1:length(car);
y=car(i);
plot(x,y,'b*-');
xlabel('year');
ylabel('car');
hold on
%级比检验
a=[];
b=[];
for i=length(car):-1:1
       a=[a;car(i)];
end

for i1=2:length(a)
       re=a(i1-1)/a(i1);
       b=[b;re];
end
tmin=exp(-2/(length(car)+1));
m=['理论最小级比为',num2str(tmin)];
disp(m)
tmax=exp(2/(length(car)+1));
rmin=min(b);
rmax=max(b);
n=['理论最大级比为',num2str(tmax)];
disp(n)
s=['实际最小级比为',num2str(rmin)];
disp(s)
pp=['理论最小级比为',num2str(rmax)];
disp(pp)
%能直接使用gm（1,1)模型
if rmin>tmin && rmax<tmax
    disp('能直接使用gm（1,1)模型');
    X11=0;
    X1=[];
for i8=1:length(a)
    X11=a(i8)+X11;
    X1=[X1;X11];
end
z11=[];
for i9=2:length(X1)
    er=0.5*X1(i9-1)+0.5*X1(i9);
    z11=[z11;-er 1];
end
B11=z11;
Y111=[];
for i10=2:length(a)
    Y111=[Y111;a(i10)];
end
J=((B11.'*B11)^-1)*(B11.')*Y111;
a11=J(1)
b11=J(2)
F1=[a(1)];
for i11=2:(length(car))
    reel=(a(1)-b11/a11)*exp(-a11*(i11-1))+b11/a11-(a(1)-b11/a11)*exp(-a11*(i11-2))-b11/a11;
    F1=[F1;reel];
end
x=2017:2021;
i13=1:length(F1);
y=F1(i13);
plot(x,y,'r*-');
xlabel('year');
ylabel('car');
hold off
%残差
E=[];
for i12=1:length(a)
    lo=(a(i12)-F1(i12))/a(i12);
    E=[E;lo];
end
e1=max(E);
pi=['最大残差为',num2str(e1)];
disp(pi)
if e1<0.1
    disp('达到较高要求')
end
if e1<0.2 && e1>0.1
    disp('达到一般要求')
end
%预测
forcast1=@(k) (a(1)-b11/a11)*exp(-a11*(k-1))-(a(1)-b11/a11)*exp(-a11*(k-2));



%不能直接使用gm（1,1)模型
else
    
    disp('不能直接使用gm（1,1)模型')
syms x
c=a+x;
d=[]; 
for i2=2:length(c)
       rec=c(i2-1)/c(i2);
       d=[d;rec];
end
k=[];
for i3=1:length(d)
    cond1=d(i3)>tmin;
    cond2=d(i3)<tmax;
    cond3=x>0;
    conds=[cond1 cond2 cond3];
    [s q c]=solve(conds,x,'ReturnConditions',true,'IgnoreAnalyticConstraints',true);
    k=[k;c];
end
[S Q C]=solve(k,x,'ReturnConditions',true);
e=char(C);
m=' < x';
f=regexprep(e,m,'');
need_to_add=str2num(f);
ADD=need_to_add+300;
kk=['需要加上',num2str(ADD)];
disp(kk)
Y=a+ADD;
Y1=[];
Y11=0;
for i4=1:length(Y)
    Y11=Y(i4)+Y11;
    Y1=[Y1;Y11];
end
Z1=[];
for i5=1:(length(Y1)-1)
    RESU=0.5*Y1(i5)+0.5*Y1(i5+1);
    Z1=[Z1;-RESU 1];
end
B1=Z1;
A=[];
for i6=2:length(Y)
    A=[A;Y(i6)];
end
u=((B1.'*B1)^-1)*B1.'*A;
a1=u(1)
b1=u(2)
F=[a(1)];
for i7=2:(length(car))
    ree=(Y(1)-b1/a1)*exp(-a1*(i7-1))+b1/a1-(Y(1)-b1/a1)*exp(-a1*(i7-2))-b1/a1;
    F=[F;ree-ADD];
end
x=2017:2021;
i=1:length(F);
y=F(i);
plot(x,y,'r*-');
xlabel('year');
ylabel('car');
hold off
%残差
E1=[];
for i13=1:length(a)
    loY=abs((a(i13)-F(i13)))/F(i13);
    E1=[E1;loY];
end
e12=max(E1);
p=['最大残差为',num2str(e12)];
disp(p)
if e12<0.1
    disp('残差＜0.1,达到较高要求')
end
if e12<0.2 && e12>0.1
    disp('残差＜0.2,达到一般要求')
end
end
%预测
forcast=@(k) (Y(1)-b1/a1)*exp(-a1*(k-1))-(Y(1)-b1/a1)*exp(-a1*(k-2))-ADD;