数学建模学习-灰色系统理论(Grey System Theory)教程(15)
写在最前
注意本文的相关代码及例子为同学们提供参考,借鉴相关结构,在这里举一些通俗易懂的例子,方便同学们根据实际情况修改代码,很多同学私信反映能否添加一些可视化,这里每篇教程都尽可能增加一些可视化方便同学理解,但具体使用时,同学们要根据实际情况选择是否在论文中添加可视化图片。
系列教程计划持续更新,同学们可以免费订阅专栏,内容充足后专栏可能付费,提前订阅的同学可以免费阅读,同时相关代码获取可以关注博主评论或私信。
目录
算法简介
灰色系统理论是邓聚龙教授于1982年首次提出的一种处理不确定性系统的方法。该理论主要用于处理部分信息已知、部分信息未知的"小样本"、"贫信息"不确定性系统。灰色系统理论中最常用的模型是GM(1,1)模型,其中第一个"1"表示模型的阶数,第二个"1"表示变量的个数。
GM(1,1)模型的基本思想是通过对原始数据序列进行累加生成(AGO)来减弱数据的随机性,使其具有较强的规律性,然后建立微分方程模型进行预测,最后通过累减还原得到预测结果。
算法特点
- 适用于小样本预测:只需要4个以上的数据点即可建模。
- 处理不确定性:能够处理部分信息未知的系统。
- 建模简单:不需要考虑数据的概率分布。
- 计算量小:相比其他预测方法计算效率高。
- 适用范围广:可用于时间序列预测、系统分析等多个领域。
理论基础
1. 灰色系统的基本概念
灰色系统理论将信息完全已知的系统称为白色系统,信息完全未知的系统称为黑色系统,而介于两者之间的系统称为灰色系统。在现实世界中,大多数系统都是灰色系统。
2. 灰色生成
灰色生成是灰色系统理论的核心概念之一,主要包括:
- 累加生成(AGO):通过累加减弱数据的随机性
- 累减生成(IAGO):AGO的逆运算,用于还原预测结果
- 紧邻均值生成:用于构建灰色微分方程
3. 灰色预测
灰色预测的基本思想是通过已知信息序列的累加生成来寻找系统发展的规律,主要特点包括:
- 不需要大量样本
- 不要求数据服从特定分布
- 计算简单,精度较高
数学原理
1. GM(1,1)模型的数学推导
设原始数据序列为:
 
     
      
       
        
         
          X
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         =
        
        
         (
        
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         1
        
        
         )
        
        
         ,
        
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         2
        
        
         )
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         n
        
        
         )
        
        
         )
        
       
       
        X^{(0)} = (x^{(0)}(1), x^{(0)}(2), ..., x^{(0)}(n))
       
      
     X(0)=(x(0)(1),x(0)(2),...,x(0)(n))
1.1 一次累加生成(AGO)
对原始序列进行一次累加:
 
     
      
       
        
         
          X
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         =
        
        
         (
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         1
        
        
         )
        
        
         ,
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         2
        
        
         )
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         n
        
        
         )
        
        
         )
        
       
       
        X^{(1)} = (x^{(1)}(1), x^{(1)}(2), ..., x^{(1)}(n))
       
      
     X(1)=(x(1)(1),x(1)(2),...,x(1)(n))
 其中:
 
     
      
       
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         =
        
        
         
          ∑
         
         
          
           i
          
          
           =
          
          
           1
          
         
         
          k
         
        
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         i
        
        
         )
        
        
         ,
        
        
         k
        
        
         =
        
        
         1
        
        
         ,
        
        
         2
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         n
        
       
       
        x^{(1)}(k) = \sum_{i=1}^k x^{(0)}(i), k = 1,2,...,n
       
      
     x(1)(k)=i=1∑kx(0)(i),k=1,2,...,n
1.2 紧邻均值生成
计算
    
     
      
       
        
         X
        
        
         
          (
         
         
          1
         
         
          )
         
        
       
      
      
       X^{(1)}
      
     
    X(1)的紧邻均值序列:
 
     
      
       
        
         
          Z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         =
        
        
         (
        
        
         
          z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         2
        
        
         )
        
        
         ,
        
        
         
          z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         3
        
        
         )
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         
          z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         n
        
        
         )
        
        
         )
        
       
       
        Z^{(1)} = (z^{(1)}(2), z^{(1)}(3), ..., z^{(1)}(n))
       
      
     Z(1)=(z(1)(2),z(1)(3),...,z(1)(n))
 其中:
 
     
      
       
        
         
          z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         =
        
        
         0.5
        
        
         (
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         +
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         −
        
        
         1
        
        
         )
        
        
         )
        
        
         ,
        
        
         k
        
        
         =
        
        
         2
        
        
         ,
        
        
         3
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         n
        
       
       
        z^{(1)}(k) = 0.5(x^{(1)}(k) + x^{(1)}(k-1)), k = 2,3,...,n
       
      
     z(1)(k)=0.5(x(1)(k)+x(1)(k−1)),k=2,3,...,n
1.3 灰色微分方程
GM(1,1)模型的基本形式为:
 
     
      
       
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         +
        
        
         a
        
        
         
          z
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         =
        
        
         b
        
       
       
        x^{(0)}(k) + az^{(1)}(k) = b
       
      
     x(0)(k)+az(1)(k)=b
 其中,
    
     
      
       
        a
       
      
      
       a
      
     
    a为发展系数,
    
     
      
       
        b
       
      
      
       b
      
     
    b为灰作用量。
对应的白化微分方程为:
 
     
      
       
        
         
          
           d
          
          
           
            x
           
           
            
             (
            
            
             1
            
            
             )
            
           
          
         
         
          
           d
          
          
           t
          
         
        
        
         +
        
        
         a
        
        
         
          x
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         =
        
        
         b
        
       
       
        \frac{dx^{(1)}}{dt} + ax^{(1)} = b
       
      
     dtdx(1)+ax(1)=b
1.4 参数估计
使用最小二乘法求解参数:
 
     
      
       
        
         
          [
         
         
          
           
            
             
              a
             
            
           
          
          
           
            
             
              b
             
            
           
          
         
         
          ]
         
        
        
         =
        
        
         (
        
        
         
          B
         
         
          T
         
        
        
         B
        
        
         
          )
         
         
          
           −
          
          
           1
          
         
        
        
         
          B
         
         
          T
         
        
        
         Y
        
       
       
        \begin{bmatrix} a \\ b \end{bmatrix} = (B^TB)^{-1}B^TY
       
      
     [ab]=(BTB)−1BTY
其中:
 
     
      
       
        
         B
        
        
         =
        
        
         
          [
         
         
          
           
            
             
              
               −
              
              
               
                z
               
               
                
                 (
                
                
                 1
                
                
                 )
                
               
              
              
               (
              
              
               2
              
              
               )
              
             
            
           
           
            
             
              1
             
            
           
          
          
           
            
             
              
               −
              
              
               
                z
               
               
                
                 (
                
                
                 1
                
                
                 )
                
               
              
              
               (
              
              
               3
              
              
               )
              
             
            
           
           
            
             
              1
             
            
           
          
          
           
            
             
              
               ⋮
              
              
               
              
             
            
           
           
            
             
              
               ⋮
              
              
               
              
             
            
           
          
          
           
            
             
              
               −
              
              
               
                z
               
               
                
                 (
                
                
                 1
                
                
                 )
                
               
              
              
               (
              
              
               n
              
              
               )
              
             
            
           
           
            
             
              1
             
            
           
          
         
         
          ]
         
        
       
       
        B = \begin{bmatrix} -z^{(1)}(2) & 1 \\ -z^{(1)}(3) & 1 \\ \vdots & \vdots \\ -z^{(1)}(n) & 1 \end{bmatrix}
       
      
     B=
              −z(1)(2)−z(1)(3)⋮−z(1)(n)11⋮1
              
Y = [ x ( 0 ) ( 2 ) x ( 0 ) ( 3 ) ⋮ x ( 0 ) ( n ) ] Y = \begin{bmatrix} x^{(0)}(2) \\ x^{(0)}(3) \\ \vdots \\ x^{(0)}(n) \end{bmatrix} Y= x(0)(2)x(0)(3)⋮x(0)(n) 
1.5 时间响应函数
求解白化微分方程得到时间响应函数:
 
     
      
       
        
         
          
           x
          
          
           ^
          
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         =
        
        
         (
        
        
         
          x
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         1
        
        
         )
        
        
         −
        
        
         
          b
         
         
          a
         
        
        
         )
        
        
         
          e
         
         
          
           −
          
          
           a
          
          
           (
          
          
           k
          
          
           −
          
          
           1
          
          
           )
          
         
        
        
         +
        
        
         
          b
         
         
          a
         
        
       
       
        \hat{x}^{(1)}(k) = (x^{(0)}(1) - \frac{b}{a})e^{-a(k-1)} + \frac{b}{a}
       
      
     x^(1)(k)=(x(0)(1)−ab)e−a(k−1)+ab
1.6 预测值还原
通过累减还原得到预测值:
 
     
      
       
        
         
          
           x
          
          
           ^
          
         
         
          
           (
          
          
           0
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         =
        
        
         
          
           x
          
          
           ^
          
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         )
        
        
         −
        
        
         
          
           x
          
          
           ^
          
         
         
          
           (
          
          
           1
          
          
           )
          
         
        
        
         (
        
        
         k
        
        
         −
        
        
         1
        
        
         )
        
       
       
        \hat{x}^{(0)}(k) = \hat{x}^{(1)}(k) - \hat{x}^{(1)}(k-1)
       
      
     x^(0)(k)=x^(1)(k)−x^(1)(k−1)
环境准备
本教程需要以下Python库:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
代码实现
1. GM(1,1)模型类定义
class GM11:
    def __init__(self, data):
        """
        初始化GM(1,1)模型
        :param data: 原始数据序列
        """
        self.raw_data = np.array(data)
        self.n = len(data)
        self.ago_data = None  # 一次累加生成序列
        self.z_data = None    # 紧邻均值生成序列
        self.a = None         # 发展系数
        self.b = None         # 灰作用量
        self.fitted = None    # 拟合值
        self.residuals = None # 残差
        self.relative_errors = None  # 相对误差
2. 累加生成和均值生成
def ago(self):
    """一次累加生成"""
    self.ago_data = np.cumsum(self.raw_data)
    return self.ago_data
def z_mean(self):
    """计算紧邻均值生成序列"""
    self.z_data = 0.5 * (self.ago_data[1:] + self.ago_data[:-1])
    return self.z_data
3. 模型拟合
def fit(self):
    """拟合GM(1,1)模型"""
    self.ago()
    self.z_mean()
    
    # 构建B矩阵和Y矩阵
    B = np.vstack([-self.z_data, np.ones(self.n-1)]).T
    Y = self.raw_data[1:]
    
    # 使用最小二乘法求解参数
    params = np.linalg.inv(B.T @ B) @ B.T @ Y
    self.a, self.b = params
    
    # 计算预测值
    self.fitted = np.zeros(self.n)
    self.fitted[0] = self.raw_data[0]
    for k in range(1, self.n):
        self.fitted[k] = (self.raw_data[0] - self.b/self.a) * (1 - np.exp(self.a)) * np.exp(-self.a * (k))
    
    # 计算残差和相对误差
    self.residuals = self.raw_data - self.fitted
    self.relative_errors = np.abs(self.residuals / self.raw_data)
    
    return self.fitted
4. 预测未来值
def predict(self, steps=1):
    """预测未来值"""
    k = np.arange(self.n, self.n + steps)
    predictions = (self.raw_data[0] - self.b/self.a) * (1 - np.exp(self.a)) * np.exp(-self.a * k)
    return predictions
实例分析
我们以某地区2015-2020年的GDP数据为例,使用GM(1,1)模型进行拟合和预测:
# 示例数据:某地区2015-2020年GDP数据(单位:亿元)
gdp_data = [1200, 1350, 1480, 1650, 1850, 2100]
# 创建并拟合GM(1,1)模型
model = GM11(gdp_data)
fitted_values = model.fit()
# 预测未来两年的GDP
predictions = model.predict(2)
结果分析
1. 拟合效果

从上图可以看出:
- 模型拟合值与实际值非常接近,说明模型拟合效果良好。
- 相对误差都控制在2%以内,最大相对误差为1.62%。
- 预测结果显示2021年和2022年的GDP预测值分别为2327.61亿元和2604.04亿元,保持了稳定增长趋势。
2. 模型评价
- 拟合精度:从相对误差来看,模型的拟合精度很高,所有点的相对误差都小于2%。
- 预测合理性:预测值保持了数据的增长趋势,且增长率与历史数据相符。
- 模型稳定性:相对误差的分布较为均匀,没有出现异常波动,说明模型稳定性好。
应用场景
GM(1,1)模型在实际中有广泛的应用,主要包括:
1. 经济领域
- GDP增长预测
- 产业结构演变分析
- 经济发展趋势预测
- 投资回报预测
2. 工程领域
- 设备故障预测
- 系统可靠性分析
- 工程项目进度预测
- 材料性能预测
3. 环境领域
- 污染物浓度预测
- 气象数据分析
- 生态系统演变预测
- 环境质量评估
4. 社会领域
- 人口增长预测
- 交通流量预测
- 能源消耗预测
- 社会发展指标预测
注意事项
- 使用GM(1,1)模型时,原始数据序列应该具有明显的指数增长或指数衰减趋势。
- 数据量不宜过大,一般4-10个数据点较为合适。
- 预测步数不宜过长,建议不超过原始数据长度的1/3。
- 如果数据波动较大,可以考虑对数据进行预处理。
- 模型预测精度与数据的光滑度有关,数据越光滑,预测精度越高。
- 建议在使用模型前进行数据平滑处理。
- 对于季节性或周期性明显的数据,需要先进行季节性调整。
- 在实际应用中,建议结合其他预测方法进行交叉验证。
同学们如果有疑问可以私信答疑,如果有讲的不好的地方或可以改善的地方可以一起交流,谢谢大家。
 
                   
                   
                   
                   
                             
                     
       
           
                 
                 
                 
                 
                 
                
               
                 
                 
                 
                 
                
               
                 
                 扫一扫
扫一扫
                     
              
             
                   2188
					2188
					
 被折叠的  条评论
		 为什么被折叠?
被折叠的  条评论
		 为什么被折叠?
		 
		  到【灌水乐园】发言
到【灌水乐园】发言                                
		 
		 
    
   
    
   
             
            


 
            