灰色预测--matlab&python实现

function SGrey

X0 = input('请输入原始负荷数据：');   %输入原始数据
n = length(X0);  %原始n年数据

%累加生成
X1 = zeros(1,n);
for i = 1:n
    if i == 1
        X1(1,i) = X0(1,i);
    else
        X1(1,i) = X0(1,i) + X1(1,i-1);
    end
end
X1

%计算数据矩阵B和数据向量Y
B = zeros(n-1,2);
Y = zeros(n-1,1);
for i = 1:n-1
    B(i,1) = -0.5*(X1(1,i) + X1(1,i+1));
    B(i,2) = 1;
    Y(i,1) = X0(1,i+1);
end
B,Y

%计算GM(1,1)微分方程的参数a和u
A = zeros(2,1);
A = inv(B'*B)*B'*Y;
a = A(1,1);
u = A(2,1);
a,u

%建立灰色预测模型
XX0(1,1) = X0(1,1);
for i = 2:n
    XX0(1,i) = (X0(1,1) - u/a)*(1-exp(a))*exp(-a*(i-1));
end
XX0
%模型精度的后验差检验
e = 0;          %求残差平均值
for i =1:n
    e = e + (X0(1,i) - XX0(1,i));
end
e = e/n;
e
aver = 0;     %求历史数据平均值
for i = 1:n
    aver = aver + X0(1,i);
end
aver = aver / n;
aver
s12 = 0;     %求历史数据方差
for i = 1:n
    s12 = s12 + (X0(1,i)-aver)^2;
end
s12 = s12 / n;
s12
s22 = 0;       %求残差方差
for i = 1:n
    s22 = s22 + ((X0(1,i) - XX0(1,i)) - e)^2;
end
s22 = s22 / n;
s22
C = s22 / s12;    %求后验差比值
C
cout = 0;    %求小误差概率
for i = 1:n
    if abs((X0(1,i) - XX0(1,i)) - e) < 0.6754*sqrt(s12)
        cout = cout+1;
    else
        cout = cout;
    end
end
P = cout / n;
P
if (C < 0.35 & P > 0.95)
    disp('预测精度为一级');
    m = input('请输入需要预测的年数: m = ');   %预测往后各年的负荷
    disp('往后m各年负荷为：');
    f = zeros(1,m);
    for i = 1:m
        f(1,i) = (X0(1,1) - u/a)*(1-exp(a))*exp(-a*(i+n-1));
    end
    f
else
    disp('灰色预测法不适用');
end

matlab输出

输入:[724.57 746.62 778.27 800.8 827.75 871.1 912.37 954.28 995.01 1037.2]

输出:

>> SGrey
请输入原始负荷数据：[724.57 746.62 778.27 800.8 827.75 871.1 912.37 954.28 995.01 1037.2
]

X1 =

  1.0e+003 *

  Columns 1 through 8

    0.7246    1.4712    2.2495    3.0503    3.8780    4.7491    5.6615    6.6158

  Columns 9 through 10

    7.6108    8.6480


B =

  1.0e+003 *

   -1.0979    0.0010
   -1.8603    0.0010
   -2.6499    0.0010
   -3.4641    0.0010
   -4.3136    0.0010
   -5.2053    0.0010
   -6.1386    0.0010
   -7.1133    0.0010
   -8.1294    0.0010


Y =

  1.0e+003 *

    0.7466
    0.7783
    0.8008
    0.8277
    0.8711
    0.9124
    0.9543
    0.9950
    1.0372


a =

   -0.0420


u =

  693.9403


XX0 =

  1.0e+003 *

  Columns 1 through 8

    0.7246    0.7398    0.7715    0.8046    0.8391    0.8750    0.9125    0.9517

  Columns 9 through 10

    0.9925    1.0350


e =

    0.1818


aver =

  864.7970


s12 =

  1.0357e+004


s22 =

   26.8113


C =

    0.0026


P =

     1

预测精度为一级
请输入需要预测的年数: m = 10
往后m各年负荷为：

f =

  1.0e+003 *

  Columns 1 through 8

    1.0794    1.1257    1.1739    1.2242    1.2767    1.3315    1.3885    1.4481

  Columns 9 through 10

    1.5101    1.5749

Python实现

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
import numpy as np
import math

history_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]
n = len(history_data)
X0 = np.array(history_data)

#累加生成
history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
X1 = np.array(history_data_agg)

#计算数据矩阵B和数据向量Y
B = np.zeros([n-1,2])
Y = np.zeros([n-1,1])
for i in range(0,n-1):
    B[i][0] = -0.5*(X1[i] + X1[i+1])
    B[i][1] = 1
    Y[i][0] = X0[i+1]

#计算GM(1,1)微分方程的参数a和u
#A = np.zeros([2,1])
A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
a = A[0][0]
u = A[1][0]

#建立灰色预测模型
XX0 = np.zeros(n)
XX0[0] = X0[0]
for i in range(1,n):
    XX0[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i));


#模型精度的后验差检验
e = 0      #求残差平均值
for i in range(0,n):
    e += (X0[i] - XX0[i])
e /= n

#求历史数据平均值
aver = 0;     
for i in range(0,n):
    aver += X0[i]
aver /= n

#求历史数据方差
s12 = 0;     
for i in range(0,n):
    s12 += (X0[i]-aver)**2;
s12 /= n

#求残差方差
s22 = 0;       
for i in range(0,n):
    s22 += ((X0[i] - XX0[i]) - e)**2;
s22 /= n

#求后验差比值
C = s22 / s12   

#求小误差概率
cout = 0
for i in range(0,n):
    if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12):
        cout = cout+1
    else:
        cout = cout
P = cout / n

if (C < 0.35 and P > 0.95):
    #预测精度为一级
    m = 10   #请输入需要预测的年数
    #print('往后m各年负荷为：')
    f = np.zeros(m)
    for i in range(0,m):
        f[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i+n))    
else:
    print('灰色预测法不适用')