注明：本文根据数学建模BOOM网课简单整理，自用

❑ 模型简介

GM(1,1)灰色预测模型

❑ 名词释义

• G：Grey（灰色），M：Model（模型），(1,1)：只含有一个变量的一阶微分方程模型

❑ 白色系统

• 典型例子：电阻器件

• 特点：内部特征是完全已知的，电压、电流和电阻之间的关系(欧姆定律)是已知的

• 知道电阻大小后，输入电压值，就能算出电流值

• 白色系统可“计算”出想要的结果

❑ 黑色系统

• 典型例子：一辆车

• 特点：内部特征是完全未知的，开车并不需要懂发动机设计和工作原理等

• 非专业人士虽然可以控制汽车，然而内部出故障时并不会修

• 黑色系统具有“不可知性”

❑ 灰色系统

• 典型例子：GDP就是灰色系统

• 特点：介于黑色和白色之间，部分已知，部分未知，具有小样本数据的不确定性系统

• 我们有往年的数据和一定的理论基础（白色）

• 但无法精确计算得出下一年的值（黑色）

• 灰色无法“计算”，但并不是完全“不可知”，可以进行“预测”

❑ 如何进行灰色预测

• 根据原始数据，通过累加等方式削弱随机性、获得有规律的新序列

• 建立相应的微分方程模型，得到离散点处的解

• 再通过累减求得的原始数据的估计值，从而对原始数据预测

❑ 适用赛题

❑ 数列预测

• 特点：定时求量，已知xx年到xx年的数据，请预测下一年的数值

• 常见GDP、人口数量、耕地面积、粮食产量等问题

• 针对的问题往往短期波动小、可预测，但长期可能变化大、难以准确预测

❑ 灾变预测

• 特点：定量求时，已知xx年到xx年的数据和某灾变的阈值，预测下一次灾变发生的时间

• 常见洪涝灾害、虫灾等问题

• 模型中需要把超出阈值的数据（异常数据）对应的时间组成新序列

❑ 拓扑预测

• 特点：对数据波形进行预测，求的是多个模型构成的模型群，等于求解多个灾变预测

• 与灾变预测类似，不过有较详细的分级，例如虫灾“轻微”“中度”“重度”

❑ 注意事项

• 需要的数据量少，而且数据量太多了没意义，例如用近100年的GDP去预测下一年毫无意义

• 只能短期预测，究竟多短没有严格限制

❑ 典型例题

❑ 例题：某城市1986到1992年道路噪声平均声级数据见下表，请预测下一年的数据

• 某城市交通噪声数据/dB(A)

 

• 特点：数据少，看不出明显规律，适合用灰色预测

❑ 这一组数看不出规律怎么办？制造规律！

• 如何制造？累加：

 

• 累加生成序列：

 

• 弱化其随机性，显现其规律性

❑问题转化

❑ 如果一个东西长得像鸭子，叫声像鸭子，走路也像鸭子，那它就是一只鸭子

❑ 生成的新序列𝑥(1) ，看起来像一个指数曲线（直线）

• 因此可用一个指数曲线的表达式来逼近这个新序列

• 相应可构建一阶常微分方程来求解拟合指数曲线的函数表达式

 

❑ 关于一阶常微分方程，假设𝑥(1)满足形式：

• 要预测下一年数值，就需要新序列𝑥(1)的表达式，那就要解出微分方程(高数基本知识)

• 要解微分方程，就要先知道参数a和u

❑ 当前问题转化为：如何求微分方程中的参数a和u

❑表达式处理

❑模型求解

❑ 最小二乘法求解

❑ 拟合值与预测值

• 因为第1到7年的噪声数据是已知的，那么通过微分方程的函数解求出的值就是拟合值

• 拟合值与实际值的偏差，代表了模型的优劣

• 偏差越小，模型越好。偏差过大说明有问题，所以下一步需要拟合值进行模型检验

❑ 模型靠不靠谱，需要进行模型检验（k从1到7都要求）

❑原始数据的级比检验

❑ 在建模最开始，需要进行数据的级比检验

• 为了确定原始数据使用GM(1,1)模型的可行性，避免白忙活，需要对原始数据进行级比检验：

• 计算

• 如果𝜆(𝑘)在区间

内，说明可用GM(1,1)模型；（数值没有明确推导，套用即可）

• 如果在区间外，可尝试平移变换

• 也就是给每个数据都加上任意常数c后看是否在区间内；求解后再减去c

• 如果尝试多次平移变换后始终无法在区间内，说明题目不适合灰色预测

❑整体思路

❑ 流程

• 核心就是最小二乘法求参数，再求微分方程

• 在最开始要对原始数据进行检验，在最后要对结果进行检验

❑代码求解

❑ 结果分析

• 求解出预测下一年的噪声值为71.4 dB

• 模型检验结果见下表，可见该模型的精确较高，可用于本题预测。

原始数据级比检验

计算

，如果所有在区间

内，说明可用GM(1,1)模型

clc,clear
x0=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]'; % 注意这里为列向量
n=length(x0);
lamda=x0(1:n-1)./x0(2:n);  % 计算级比，为列向量
range=minmax(lamda');      % 计算级比的范围
lamrange_min=exp(-2/(n+1));     % 允许的范围下界
lamrange_max=exp(2/(n+1));      % 允许的范围上界
% 原始数据的级比检验
if (range(1)>lamrange_min && range(2)<lamrange_max)
    fprintf("原始数据通过级比检验")
else
    fprintf("原始数据未能通过级比检验！！！")
    return
end

最小二乘法求微分方程参数

离散形式下，且均值处理后的方程：

其中

以方程矩阵形式Y=BU求出U（对应最小二乘法部分）

x1=cumsum(x0)  %累加序列
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
% 注意此处是反斜杠，即左除
u=B\Y       % 拟合参数u(1)=a,u(2)=b

求解微分方程

以符号解的形式求出，就是求出表达式

注意：代码中的x是模型中的累加序列，x0是模型中的原始序列

此处代码中累加序列函数表达式的自变量取值与模型中的下标不一致

代码中的累加序列x是以函数的形式表示的，x(0)实际意义是累加序列第1年的数值！！！所以对应地，x(1)代表第2年，一直到x(n-1)年代表第n

即代码中的x(k)，对应模型中的

syms x(t)   %定义符号变量
%dsolve求微分方程的符号解，即累加序列的函数表达式，而不是数值；solve是数值
x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1));     % diff是一阶导或差分，初始值x0(1)原始矩阵的第一个值
xt=vpa(x,6) %以保留6位小数（有效小数）的格式显示微分方程的解,vpa设置精度
% 函数subs(x,t,[0:n-1])意味着把t=0到t=n-1依次代入x的表达式求出值
nihe1=subs(x,t,[0:n-1]); %求已知数据对应的拟合值
nihe1=double(nihe1);     %符号数转换成数值类型，否则无法作差分运算
nihe=[x0(1),diff(nihe1)] %差分运算，得到原始序列拟合值

预测下一年的道路噪声平均声级，模型中：

对应代码中：t=7时的x，减去t=6时的x

因为x是函数表达式，需要用subs代入自变量取值，再用double转为数值

yuce_result=double(subs(x,t,7))-double(subs(x,t,6))  % 原始序列的预测值，即第8年噪声

结果检验

delta=abs((x0'-nihe )./x0')  %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'  %计算级比偏差值