灰色预测算法改进_改进背景值Z_python

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

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

参考文献：

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

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

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

1、算法

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

z ( k ) = α x ( 1 ) ( k − 1 ) + ( 1 − α ) x ( 1 ) ( k )                        ( 1 ) z(k)=\\alpha x^{(1)}(k-1)+(1-\\alpha)x^{(1)}(k)~~~~~~~~~~~~~~~~~~~~~~(1) z(k)=αx(1)(k−1)+(1−α)x(1)(k)                      (1)

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

α = 1 a − 1 exp ⁡ ( a ) − 1                                        ( 2 ) \\alpha = \\frac{1}{a}-\\frac{1}{\\exp(a)-1}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(2) α=a1​−exp(a)−11​                                      (2)

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

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

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

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

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

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

• 5）将a代入(2)，计算 α n e w \\alpha_{new} αnew​，若 ∣ α n e w − α ∣ > ε |\\alpha_{new}-\\alpha|>\\varepsilon ∣αnew​−α∣>ε，则令 ∣ α = ∣ α n e w |\\alpha=|\\alpha_{new} ∣α=∣αnew​；循环3) 到 5）的步骤，若 ∣ α n e w − α ∣ ≤ ε |\\alpha_{new}-\\alpha|\\leq\\varepsilon ∣αnew​−α∣≤ε，则进入步骤6)

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

  ---

2、 案例及代码

2.1 数据

同样以江苏省无锡市锡北镇电力负荷预测为例，令 ε = 1 E − 10 \\varepsilon=1E-10 ε=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)， α = 0.5 \\alpha=0.5 α=0.5（可以令上述算法的 ε = 10 \\varepsilon=10 ε=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）令 ε = 1 E − 10 \\varepsilon=1E-10 ε=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 %