基于振荡序列的灰色GM（1，1|sin）幂模型及其应用

曾亮

（广东理工学院基础部，广东肇庆526100）

灰色系统理论[1]是研究小样本、贫信息和不确定问题的一种新方法，灰色预测理论是其核心体系之一，因其不要求序列具备典型的分布规律，且计算量较小，所以，在小样本序列预测上具有一定优势。传统的灰色预测模型通常对具有近似指数律的单调递增或递减序列模拟预测精度较高，在各领域已得到广泛应用[2]，而对于振荡序列的模拟预测往往难以达到预期的效果。因此，关于如何提高振荡序列的模拟预测精度是值得深入研究的方向，近年来，已取得一些有价值的研究成果，集中体现在两方面：一是对传统GM(1,1)模型进行拓展以适应振荡序列的预测；二是对原始振荡序列变换处理后建立传统GM(1,1)模型。关于模型拓展方面，DENG[3]提出了摆动型灰色GM(1,1|tan(k-π)p,sin(k-π)p)模型，由于模型中参数取值缺乏严格的理论依据，具有较强的随意性，所以对该模型的后续研究极少；MAO等[4]借鉴文献[3]中三角函数能反映振荡性的思路，提出了GM(1,1|sin)模型，并将该模型应用于交通流预测，取得了较高的模拟精度；ZHOU等[5]提出了一种改进的灰色GM(0,1|sin)振荡序列模型，并运用灰色理论，引入变权强化缓冲算子，削弱扰动项对原始数据序列的干扰；许泽东等[6]在文献[4]的基础上，将一阶累加生成拓展为分数阶累加，提出了分数阶累加GM(1,1|sin+cos)模型，并通过实例说明该模型具有较强的适用和拟合性能。由于GM(1,1)幂模型[7]是一类非线性灰色模型，通过幂指数取值的变化能适应具有不同波动特性的原始序列，近年来，关于利用GM(1,1)幂模型对振荡序列进行预测的研究成果较为丰富。比如，王正新等[8]在传统GM(1,1)幂模型的基础上通过优化幂指数的方法提高模型对小样本振荡序列的预测精度；王正新等[9]对GM(1,1)幂模型的病态性进行了研究，指出部分情形下可能出现数据矩阵求逆的病态性；王正新[10]在GM(1,1)幂模型的定义型和白化型基础上，推导了5种派生型GM(1,1)幂模型，构建了GM(1,1)幂模型群，并将成果推广到灰色多变量模型，提出了灰色多变量GM(1,N)幂模型及其派生GM(1,N,x(1))幂模型[11]；王正新[12]提出了含有系统延迟和时变参数的振荡GM(1,1)幂模型，并将该模型应用于应急资源需求预测，结果表明，其优于传统GM(1,1)幂模型和ARIMA方法；王正新[13]通过对原始序列建立GM(1,1)幂模型，并利用傅立叶级数提取模型残差中的周期性振荡规律，有效提高了对小样本振荡序列的预测精度；王俊芳等[14]借助分数阶累加算子，构建了分数阶离散GM(1,1)幂模型，并提出了一种新的累加阶数及幂参数的确定方法，提高了模型的预测精度；WANG等[15]利用傅立叶残差修正的GM(1,1)幂模型对我国高新技术产品国际贸易量进行了预测；PEI等[16]提出了一种基于非线性最小二乘法求解GM(1,1)幂模型参数的新方法，并以我国高新技术企业员工需求预测为例，对该模型进行了验证。关于对原始振荡序列的变换处理，钱吴永等[17]提出了先通过加速平移变换将振荡序列转化为单调序列，然后再进行加权均值生成变换，对变换后的序列建立GM(1,1)模型。算例表明，该方法能提高GM(1,1)模型的预测精度；ZENG等[18]研究了一种可压缩振荡序列幅值的平滑算法，并在此基础上推导出一种基于振荡序列的灰色预测模型，通过改善振荡序列的平滑性来提高灰色模型的模拟精度。

本文在以上研究的基础上，进一步研究振荡序列的模拟预测问题。为更好地描述由多个不同周期运动叠加下系统的振荡特征，将文献[4]中GM(1,1|sin)模型的灰作用量sinpk项更改为(sinpk)γ，构建了GM(1,1|sin)幂模型。通过调整正弦函数项的幂指数，使其更适应振荡序列的变化，从而达到提高模拟预测精度的目的。由于增加了幂指数，模型的白化微分方程形式较为复杂，因此，难以得到模型的精确时间响应函数。为此，本文利用高精度的复合3/8辛普森积分公式得到模型的近似时间响应式。为得到最优的模拟预测效果，以平均模拟相对误差为适应度函数，利用粒子群优化算法，求得函数值最小化下的模型参数。最后，分别以城市交通流和高新技术产品出口额的模拟预测为例，并与传统GM(1,1)模型、GM(1,1)幂模型和GM(1,1|sin)模型进行了比较，以检验本文模型的有效性和实用性。

1 灰色GM（1，1|sin）模型简介

定义1设非负系统行为序列X={x(1),x(2),…,x(n)}，若

（1）对∀k=1,2,…,n-1，有x(k+1)-x(k)＞ 0，则称X为单调递增序列；

（2）对∀k=1,2,…,n-1，有x(k+1)-x(k)＜0，则称X为单调递减序列；

（3）若∃k1,k2∈ {1,2,…,n-1}，使得x(k1+1)-x(k1)＞ 0和x(k2+1)-x(k2)＜ 0同时成立，则称X为振荡序列。

针对振荡序列的模拟预测，文献[4]给出了一种灰作用量为bsinpk+c的GM(1,1|sin)模型。

定义2若X(0)={x(0)(1),x(0)(2),…,x(0)(n)}为非负原始数据序列，X(1)={x(1)(1),x(1)(2),…,x(1)(n)}为X(0)的一阶累加生成序列，为X(1)的紧邻均值生成序列，即z(1)(k)=0.5(x(1)(k-1)+x(1)(k))，则称

为GM(1,1|sin)模型。微分方程

为GM(1,1|sin)模型的白化形式。

由最小二乘辨识方法可得GM(1,1|sin)模型参数=(a,b,c)T的估计值，若令

则=(a,b,c)T=(BTB)-1BTY。

定理1假设B、Y和如定义2所述，即=(a,b,c)T=(BTB)-1BTY，则GM(1,1|sin)模型的时间响应序列为

还原式为

证明过程参见文献[4]，此略。

2 灰色GM（1，1|sin）幂模型

定义3若非负原始序列X(0)={x(0)(1),x(0)(2),…,x(0)(n)}，X(1)为X(0)的一阶累加生成序列，Z(1)为X(1)的紧邻均值生成序列，则称

为GM(1,1|sin)幂模型。微分方程

为GM(1,1|sin)幂模型的白化形式。

从GM(1,1|sin)幂模型的表达式中可看出，当γ=0时，GM(1,1|sin)幂模型为传统 GM(1,1)模型；当γ=1时，GM(1,1|sin)幂模型为文献[4]中的GM(1,1|sin)模型。

在GM(1,1|sin)幂模型中，当参数p和幂指数γ给定时，可利用最小二乘法估计模型的其他参数。

定理2若令

当|BTB|≠0时，由最小二乘法可得模型参数=(a,b,c)T的辨识值为

证明由矩阵B的表达式可得

则BTB的行列式值为

BTB的伴随矩阵为

当|BTB|≠0时，可得BTB的逆矩阵为

又由于BTY=

则由最小二乘法可得模型参数=(a,b,c)T的辨识值为

定理3设B，Y，如定理2所述，m为正整数，则步长为的复合3/8辛普森积分公式下的GM(1,1|sin)幂模型的近似时间响应序列为

还原值为证明由微分方程

可解得

则

由初始条件(1)(1)=x(1)(1)=x(0)(1)，得

将式（12）代入式（11），得

则近似时间响应序列为

由一阶累减还原公式可得还原值：

(0)(k)=(1)(k)-(1)(k-1)，k=2,3,…。

3 GM（1，1|sin）幂模型参数p和γ的确定

平均模拟或预测相对误差（ARPE）常被用来检验模型的模拟预测效果，其公式为

其中n为模拟或预测数据的个数。

由模型（6）可知，参数p和γ为模型的非线性部分，确定其精确值比较困难。本文利用平均模拟相对误差最小化方法确定p和γ的值，即建立以下优化模型：

为求解该优化模型，本文采用粒子群优化算法(PSO)[19]确定最优参数。粒子群优化算法是一种智能进化计算技术，首先对被优化的函数初始化一组随机解，每一个个体解称为一个“粒子”。所有的粒子都有一个由被优化的函数决定的适应值，每个粒子还有一个决定其搜索距离和方向的速度，然后就追随当前的最优粒子在解空间搜索，通过不断迭代找到最优解[20-22]。在迭代过程中，每个粒子都根据2个极值来更新自己的位置，一个是粒子本身目前搜索到的最优解，称为个体极值(pbest)，另一个就是群中所有粒子目前搜索到的最优解，称为全局极值(gbest)。每个粒子按以下公式不断更新自己的速度和位置：

其中，vk是粒子的速度向量；xk是粒子的当前位置；pbestk是粒子本身搜索到的最优解，即个体极值；gbestk是粒子群搜索到的最优解，即全局极值；ω为惯性权重，用来保持全局搜索和局部搜索的平衡，其值取0.1～0.9；c1和c2是加速系数，值为0～4，通常取2。

利用粒子群优化算法计算GM(1,1|sin)幂模型参数的步骤如下：

步骤1在可行域内初始化粒子种群和粒子速度；

步骤2以平均模拟相对误差为适应度函数，计算每个粒子的适应值，并选择个体极值和全局极值；

步骤3若计算结果达到一定的精度或达到最大迭代次数，转步骤5，否则转步骤4；

步骤4分别按照式（16）和（17）更新每个粒子的位置和速度，然后转步骤2；

步骤5输出全局极值及相应的解，退出循环，算法结束。

4 应用实例

4.1 城市道路交通流预测

文献[4]给出了某天6:00～18:00(纽约时间，以1 h为间隔)纽约市的道路交通流量数据，见表1。

表1 某天6：00～18：00纽约市道路交通流量Table1 City traffic flow from 6 ：00 to 18：00 in New York

以6:00～18:00的交通流量为原始数据，分别建立 GM(1,1)模型、GM(1,1)幂模型[8]、GM(1,1|sin)模型[4]和本文提出的GM(1,1|sin)幂模型，4种模型的模拟值及相对误差见图1和表2。

图1 4种模型对道路交通流的模拟比较图Fig.1 Comparison of simulation results of four models for road traffic flow

从图1中可看出，原始数据为振荡序列，本文模型模拟序列的变化趋势与原始序列的变化趋势较为吻合，模拟效果要好于另外3种模型。从表2中可看出，本文模型的平均模拟相对误差为2.140 2%，显著低于GM(1,1)模型的10.565 2%、GM(1,1)幂模型的7.730 2%和GM(1,1|sin)模型的4.465 1%。综合分析可知，本文提出的模型较3种经典模型具有更高的模拟预测精度，更适合振荡序列的模拟预测。

表2 4种模型对城市交通流的模拟效果比较Table2 Comparison of the simulation effects of four models for city traffic flow

4.2 高新技术产品出口额预测

表3为2007―2015年广东省高新技术产品出口额数据(来源于广东科技统计网：http://www.sts.gd.cn/)。

以2007―2015年广东省高新技术产品出口额为原始数据，分别建立GM(1,1)模型、GM(1,1)幂模型、GM(1,1|sin)模型和GM(1,1|sin)幂模型，4种模型的模拟值及相对误差结果见图2和表4。

表3 2007―2015年广东省高新技术产品出口额Table3 The gross of exports of high and new technology products in Guangdong province form 2007 to 2015

表4 4种模型对高新技术产品出口额的模拟效果比较Table4 Comparison of the simulation effects among four models for the gross of exports of high and new technology products

图2 4种模型对高新技术产品出口额的模拟比较图Fig.2 Comparison of the simulation results among four models for the gross of exports of high and new technology products

从图2中可以看出，原始数据为振荡序列时，本文模型的模拟效果要好于其他3种模型。从表4可中看出，GM(1,1)模型、GM(1,1)幂模型和GM(1,1|sin)模型的平均模拟相对误差分别为6.601 0%，5.281 2%和4.960 9%，而本文模型的平均模拟相对误差仅为2.114 6%，其精度显著高于前三者。综合分析可知，本文模型较3种经典模型更适合对广东省高新技术产品出口额的模拟预测，进一步证实了本文模型对振荡序列的适应性。

5 结 论

提出了一种针对振荡序列模拟预测的灰色GM(1,1|sin)幂模型，可以根据原始序列的振荡特征，灵活调整幂指数的取值，即调整预测曲线的形状。由于无法给出模型的解析解，用精度较高的复合3/8辛普森积分公式代替模型的近似时间响应式。为了得到较好的模拟预测效果，采用粒子群算法确定模型参数的最优值。通过2个应用实例证明了本文模型对振荡序列具有较高的预测精度，拓展了灰色预测模型的应用范围。