[人工智能]灰色预测算法改进

😊😨😱😡:我的心路历程，妈耶，妥妥的翻车了，论文【1】里改进后效果明显变好了，采用了新的数据集，效果反而变差了，但仁者见仁吧，可能真的在别的数据集上效果就会好呢。后续会阅读别的论文，持续更新新的改进方法。
29/03/2022 10:44

灰色系统理论及其应用系列博文：
一、灰色关联度分析法(GRA)_python
二、灰色预测模型GM(1,1)
三、灰色预测模型GM(1,n)
四、灰色预测算法改进

参考文献：

[1] 改进灰色预测模型在电力负荷预测中的应用[J], 陈磊，张青云，向晓，刘红云，刘闯.

传统灰色预测模型可从两个方面入手优化，第一个是对原始数据的改造，即改善其平滑性，使其更接近指数发展规律；第二个是对模型参数的优化，即对求解灰参数时所用背景值进行优化改造。

一、修改背景值的参数设定

1、算法

GM（1，1）计算时，背景值由下式计算

$$z(k)=\alpha x^{(1)}(k?1)+(1?\alpha )x^{(1)}(k)(1)$$

在传统的计算中，一般将$\alpha$值取为0.5。ＧＭ（１，1）微分方程中的发展灰数ａ 与参数$\alpha$的关系为

$$\alpha =\frac{1}{a}?\frac{1}{\exp (a)?1}(2)$$

因此有ａ 与参数$\alpha$的关系如下：

由表１可得，当｜a｜较小时，$\alpha$非常接近０．５；当｜a｜较大时，$\alpha$偏离0.5较大。所以，将$\alpha$的值取0.5，显然是不合理的，需要通过a 的值来确定$\alpha$的值。
本文提出了 ＧＭ（１，ｎ）模型的$\alpha$参数修正预测方法，主要步骤如下。

• 1）初始化$\alpha$=0.5

• 2）对原始序列$x^{(0)}$计算累加值$x^{(1)}$

• 3）根据(1)生成背景值序列

• 4）构建A,B,利用最小二乘求解a,u

• 5）将a代入(2)，计算${\alpha }_{new}$，若$|{\alpha }_{new}?\alpha |>\epsilon$，则令$|\alpha =|{\alpha }_{new}$；循环3) 到 5）的步骤，若$|{\alpha }_{new}?\alpha |\le \epsilon$，则进入步骤6)

• 6)求解微分方程，得到$\widehat{x^{(1)}}$，累减还原，得到$\widehat{x^{(0)}}$

  ---

2、 案例及代码

2.1 数据

同样以江苏省无锡市锡北镇电力负荷预测为例，令$\epsilon =1E?10$给出灰度预测的结果

2.2 代码

• 代码

  # -*- coding: utf-8 -*- 
  
  class GM11():
      def __init__(self):
          self.f = None
          self.alpha = 1 / 2
          self.alpha_b = 1
  
      def isUsable(self, X0):
          '''判断是否通过光滑检验'''
          X1 = X0.cumsum()
          rho = [X0[i] / X1[i - 1] for i in range(1, len(X0))]
          rho_ratio = [rho[i + 1] / rho[i] for i in range(len(rho) - 1)]
          print(rho, rho_ratio)
          flag = True
          for i in range(2, len(rho) - 1):
              if rho[i] > 0.5 or rho[i + 1] / rho[i] >= 1:
                  flag = False
          if rho[-1] > 0.5:
              flag = False
          if flag:
              print("数据通过光滑校验")
          else:
              print("该数据未通过光滑校验")
  
          '''判断是否通过级比检验'''
          lambds = [X0[i - 1] / X0[i] for i in range(1, len(X0))]
          X_min = np.e ** (-2 / (len(X0) + 1))
          X_max = np.e ** (2 / (len(X0) + 1))
          for lambd in lambds:
              if lambd < X_min or lambd > X_max:
                  print('该数据未通过级比检验')
                  return
          print('该数据通过级比检验')
  
      def train(self, X0):
          X1 = X0.cumsum()
          Z = (np.array([-((1 - self.alpha) * X1[k - 1] + self.alpha * X1[k]) for k in range(1, len(X1))])).reshape(
              len(X1) - 1, 1)
          # 数据矩阵A、B
          A = (X0[1:]).reshape(len(Z), 1)
          B = np.hstack((Z, np.ones(len(Z)).reshape(len(Z), 1)))
          # 求灰参数
          a, u = np.linalg.inv(np.matmul(B.T, B)).dot(B.T).dot(A)
          # u = Decimal(u[0])
          # a = Decimal(a[0])
          self.alpha_b = self.alpha  # 上一次的alpha
          self.alpha = 1 / a - 1 / (np.exp(a) - 1)  # 这一次的alpha
          print("上一次的alpha：", self.alpha_b, "本次的alpha：", self.alpha)
          print("灰参数a：", a, "，灰参数u：", u)
          self.f = lambda k: (X0[0] - u / a) * np.exp(-a * k) + u / a
  
      def predict(self, k):
          X1_hat = [float(self.f(k)) for k in range(k)]
          X0_hat = np.diff(X1_hat)
          X0_hat = np.hstack((X1_hat[0], X0_hat))
          return X0_hat
  
      def evaluate(self, X0_hat, X0):
          '''
          根据后验差比及小误差概率判断预测结果
          :param X0_hat: 预测结果
          :return:
          '''
          S1 = np.std(X0, ddof=1)  # 原始数据样本标准差
          S2 = np.std(X0 - X0_hat, ddof=1)  # 残差数据样本标准差
          C = S2 / S1  # 后验差比
          Pe = np.mean(X0 - X0_hat)
          temp = np.abs((X0 - X0_hat - Pe)) < 0.6745 * S1
          p = np.count_nonzero(temp) / len(X0)  # 计算小误差概率
          print("原数据样本标准差：", S1)
          print("残差样本标准差：", S2)
          print("后验差：", C)
          print("小误差概率p：", p)
  
  def MAPE(y_true, y_pred):
      """计算MAPE"""
      n = len(y_true)
      mape = sum(np.abs((y_true - y_pred) / y_true)) / n * 100
      return mape
  
  
  if __name__ == '__main__':
      import matplotlib.pyplot as plt
      import numpy as np
      import pandas as pd
  
      plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
      plt.rcParams['axes.unicode_minus'] = False  # 步骤二（解决坐标轴负数的负号显示问题）
      data = pd.read_csv("data.csv")
      X = np.array(
          [21.2, 22.7, 24.36, 26.22, 28.18, 30.16, 32.34, 34.72, 37.3, 40.34, 44.08, 47.92, 51.96, 56.02, 60.14,
           64.58,
           68.92, 73.36, 78.98, 86.6])
      # 训练集
      X_train = X[::]
      # 测试集
      X_test = []
  
      model = GM11()
      model.isUsable(X_train)  # 判断模型可行性
      thresh = 1E-10
  
      while abs(model.alpha - model.alpha_b) > thresh :
          model.train(X_train)  # 训练
      Y_pred = model.predict(len(X))  # 预测
      Y_train_pred = Y_pred[:len(X_train)]
      Y_test_pred = Y_pred[len(X_train):]
      score_test = model.evaluate(Y_train_pred, X_train)  # 评估
  
      # 可视化
      plt.grid()
      plt.plot(np.arange(len(X_train)), X_train, '->')
      plt.plot(np.arange(len(X_train)), Y_train_pred, '-o')
      plt.legend(['负荷实际值', '灰色预测模型预测值'])
      print("mape:",MAPE(Y_pred,X),"%")
      plt.title('训练集')
      plt.show()
  

2.3 效果对比

1）原始的GM(1,1)，$\alpha =0.5$（可以令上述算法的$\epsilon =10$，就等同于原始的GM(1,1)了）

数据通过光滑校验
该数据通过级比检验
上一次的alpha： 0.5 本次的alpha： [0.50620611]
灰参数a： [-0.07448023] ，灰参数u： [20.13471058]
原数据样本标准差： 20.132184863258658
残差样本标准差： 0.5563201606178899
后验差： 0.0276333723535977
小误差概率p： 1.0
mape: 0.8581559747654708 %

2）令$\epsilon =1E?10$给出灰度预测的结果：

数据通过光滑校验
该数据通过级比检验
上一次的alpha： 0.5 本次的alpha： [0.50620611]
灰参数a： [-0.07448023] ，灰参数u： [20.13471058]
上一次的alpha： [0.50620611] 本次的alpha： [0.50620325]
灰参数a： [-0.07444585] ，灰参数u： [20.12539706]
上一次的alpha： [0.50620325] 本次的alpha： [0.50620325]
灰参数a： [-0.07444586] ，灰参数u： [20.12540136]
上一次的alpha： [0.50620325] 本次的alpha： [0.50620325]
灰参数a： [-0.07444586] ，灰参数u： [20.12540136]
原数据样本标准差： 20.132184863258658
残差样本标准差： 0.5580281824309876
后验差： 0.027718212713682754
小误差概率p： 1.0
mape: 0.8761683399209106 %