一、簡介

1 名稱由來

灰色模型（Gray Model）,鄧聚龍教授1982年提出。
常見系統分類：
白色系統是指一個系統的內部特徵是完全 已知的，即系統的資訊是完全充分的。
黑色系統是指一個系統的內部資訊對外界來說是一無所知的，只能通過它與外界的聯絡來加以觀測研究。
灰色系統內的一部分資訊是已知的，另一部分資訊是未知的，系統內各因素間有不確定的關係。

2 演算法原理

二、原始碼

clc;                                                                  % 清屏
clear;                                                                % 清理記憶體
close all;                                                            % 關閉所有視窗
I_2=load('武漢地區PM2.5資料.txt');                                                     % 航空燃油泵電流（單位：A）
% Liuliang=[0 2000 4000 6000 8000 10000 12000 14000 16000 18000];
% LL=0:600:18000;
% I1=interp1(Liuliang,I0,LL,'pchip');
% I=I1';
% I2=I1(1:14,3);
% I3=I1(1:14,4);
% I4=[I2',I3'];
% I_1=I4';
I_1=I_2+172;
[a,b]=size(I_1);                                                        % 計算I的行列大小
for i=1:b                                                             % 求取原始資料的一次累加序列
    I_1_AGO(:,i)=cumsum(I_1(:,i));
end
I_1_AGO;                                                              % 累加序列，ok
[c,d]=size(I_1_AGO);                                                  % 生成序列Y
y=I_1(2:a,:);
Y=y;
beta=0.5;                                                             % 一般取為0.5
for k=1:d                                                             % 生成序列L
    for j=2:c
        l(j-1,k)=(beta)*I_1_AGO(j,k)+(1-beta)*I_1_AGO(j-1,k);         %構建緊鄰均值
    end
end
l;
[e,f]=size(l);
l1=ones(e,1);                                                         % 生成e*1的全1矩陣
L=[-l,l1];                                                            %生成矩陣B
A1=inv(L'*L)*L'*Y;                                                    % 由最小二乘法生成n元一階微分方程的兩個待定引數矩陣A和B，從而對原灰色模型進行白化
A1=A1';
[g,h]=size(A1);
A=A1(:,1:h-1);                                                        % 求發展係數（反映了預測值序列的發展態勢。為負值時（當然指的是最小二乘求取法），表示指數衰減；為正值時，表示指數上升）。
B=A1(:,end);                                                          % 求灰色作用量（反映資料變化的關係且是內涵外延化的具體體現，是區別灰色建模和一般建模過程的根本途徑）。
disp(['發展係數:',num2str(A)])
disp(['灰色作用量:',num2str(B)])

for k=0:a+1600                                                             % 累加序列預測值
    II(k+1)=(I_1(1)-B/A)*exp(-A*k)+B/A;                                 % 通過時間響應序列求預測後累加序列In_1(k-1,:)=(expm(A*(k-1)))*(L'+inv(A)*B)-inv(A)*B;
end
    II;
for kk=1:a+1600 
    yuce(kk)=II(kk+1)-II(kk);                                         % 通過累減操作後即可得到所求預測值
end
yuce';                                                               % 模擬結果
%wucha(:,1)=I_1'-yuce;                                                   % 誤差
%Jieguo=[I_1,yuce',wucha(:,1)];                                          % [原始資料，模擬結果，誤差]
%ave=(sum(abs(wucha))/a)*100                                           % 平均相對誤差
x=1:a+1600;
hold on
box on
%grid on
plot(x(1:a),I_1,'k.-')
hold on

plot(x(1:a+1600),yuce,'b.-')
plot(35*ones(2500,1),'r-')
% plot(x,I_1,'k',x,yuce,'b')                                              % 繪圖
legend('原始資料','預測資料','最終年終平均濃度')
% title('基於GM(1,1)模型的複合材料結構損傷趨勢分析')
xlabel('治理時長/天');
ylabel('PM2.5');
xlim([1,2000]);
ylim([0,650]);
axis fill;
hold off

yucend=[280,0,0,0,0,0]';
for i=1:1:5
    yucend(i+1,1)=yuce(365*i);
    data1=data(:,1)+171;
qt=230;%  控制建模資料個數
for at=qt:length(data)-1
    at
I=data1(w:at);
[a,b]=size(I);  % 計算I的行列大小
for i=1:b                                                             % 求取原始資料的一次累加序列
    I_1_AGO=cumsum(I);       %cumsum(a) 求陣列的累加和
end
I_1_AGO;                                                            % 累加序列，ok
[c,d]=size(I_1_AGO);                                                  % 生成序列Y
y=I(2:a,:);    % 從第2個開始取出原始資料構成Y矩陣
Y=y;
beta=0.5;                                                             % 一般取為0.5
for k=1:d                                                             % 生成序列L
    for j=2:c
        l(j-1,k)=(beta)*I_1_AGO(j,k)+(1-beta)*I_1_AGO(j-1,k);
    end
end
l;   % 計算緊鄰均值矩陣  30行1列
[e,f]=size(l);    
l1=ones(e,1);                                                         % 生成e*1的全1矩陣
L=[-l,l1];        % 構建 B矩陣
A1=inv(L'*L)*L'*Y;                                                    % 由最小二乘法生成n元一階微分方程的兩個待定引數矩陣A和B，從而對原灰色模型進行白化
A1=A1';
[g,h]=size(A1)
A=A1(:,1:h-1) %將a b 的值提取出來                            % 求發展係數（反映了預測值序列的發展態勢。為負值時（當然指的是最小二乘求取法），表示指數衰減；為正值時，表示指數上升）。
B=A1(:,end)                                                         % 求灰色作用量（反映資料變化的關係且是內涵外延化的具體體現，是區別灰色建模和一般建模過程的根本途徑）。
disp(['發展係數:',num2str(A)])
%Fazhanxishu=A
disp(['灰色作用量:',num2str(B)])
Huisezuoyongliang=B
for k=0:a                                                             % 累加序列預測值
    II(k+1)=(I(1)-B/A)*exp(-A*(k))+B/A;                                 % 通過時間響應序列求預測後累加序列In_1(k-1,:)=(expm(A*(k-1)))*(L'+inv(A)*B)-inv(A)*B;
end
 II;      % 算出X（1）的
for kk=1:a
    yuce(kk+1)=II(kk+1)-II(kk);                                         % 通過累減操作後即可得到所求預測值
end
yuce1(w,:)=yuce;
w=w+1;
end
yuce1;
lyuce=[data1(1,1),yuce1(1,2:qt),yuce1(1:length(data)-qt,qt+1)'];
[o,p]=size(lyuce);
wucha=data1-lyuce';                                                   % 誤差
%Jieguo=[I,lyuce,wucha(:,1)];                % [原始資料，模擬結果，誤差]
 

三、執行結果

四、備註

版本：2014a

完整程式碼或代寫加1564658423