一种基于量子引力算法的时间序列预测方法技术领域
本发明属于信号处理及信号预测等领域，更为具体地讲，涉及一种基于量
子引力算法的时间序列预测方法。
背景技术
非线性时间序列预测方法自20世纪80年代中期开始发展,近年来得到了
更深入的研究和更广泛的应用。其中，Volterra级数模型分析方法因其输出是其
滤波器核的线性函数,可用现有的线性工具来分析其滤波性能,已成为广泛使用
的预测模型之一。Volterra级数是一种泛函数,大多数非线性动态系统都可以用
Volterra级数逼近到任意准确的程度。换句话说,一大类的非线性动态系统都有它
固有的Volterra核函数。Volterra核函数可以完全表征非线性系统的性质。因此,
只要知道核就足以确定系统对任何输入激励的响应。但是，建立非线性系统
Volterra级数模型,不是一件容易的事情，其高阶的核函数估计是其面临的一个
最大的困难，这在很大程度上制约了Volterra泛函级数模型的有效应用。而现
有的Volterra核函数辨识主要采用传统的最小二乘法，该法采用梯度信息进行搜
索，易陷入局部最小，得不到满意结果，尤其在干扰严重情况下，辨识的稳定
性并不理想。其结果就导致了建模的不精确,影响模型的预测性能，尤其是当训
练集中包含的与预测轨道相关的信息过少或者夹杂许多无关的信息时极易影响
训练效果，造成模型训练不充分或过拟合,导致多步预测性能不佳。
而量子引力算法(QuantumInspiredGravitationalSearchAlgorithm,简称
QGSA)作为一种新型的智能优化算法，其从量子力学角度提出的改进的万有引
力搜索算法(GravitationalSearchAlgorithm，简称GSA)模型，将粒子定义在一
个量子空间内，该量子空间由概率密度函数决定(处于量子束缚态的粒子能够
以一定的概率出现在任何空间点)，认为粒子具有量子的行为，能在整个可行解
空间中进行搜索，对需求解的具体问题本身没有特别的要求，因而全局搜索能
力较强，克服了上述Volterra核函数辨识方法的一些不足，另外，还具有容易实
现、收敛快且参数少等优点。因此，研究基于量子引力算法的Volterra核函数辨
识及其时间序列预测方法具有重要的理论价值和现实意义。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于量子引力算法的时
间序列预测方法，将量子引力算法引入非线性系统的Volterra核函数辨识中，以
提高电子产品退化状态预测的准确度。
为实现上述发明目的，本发明基于量子引力算法的时间序列预测方法，其
特征在于，包括以下步骤：
(1)、数据预处理
(1.1)、设预选的数据集DS＝{x(1),x(2),…,x(D)}，其中D为数据个数；
(1.2)、根据微分熵率的方法确定数据集DS的最佳嵌入维数d和延迟时间
τ；
(1.3)、利用Matlab中的函数windowize对数据集DS进行分帧重构，并将
各帧数据分别映射到d维特征空间，且以Volterra级数的方式p(p≥1)阶截断，从
而得到N0(N0＝D-d×τ)个数据帧{(Xt,Yt)}(t＝1,2,…,N0)，并作为预测模型的输入和目标输出
数据，其中：输入数据为Xt＝{1,x(t),x(t+τ),…,x(t+(d-1)×τ),x2(t),x(t)×x(t+τ),…,x2(t+(d-1)×τ),…}目
标输出数据为Yt＝x(t+d×τ)(t＝1,2,…,N0)；
(2)、初始粒子的局部最佳位置
设有NS个粒子，粒子位置用M维向量Wi表示，粒子群位置可用矩阵表
示，i＝1,2,…,NS；随机初始化NS个粒子的初始位置Wi，，并令各粒子局部最佳
位置为Lbesti＝Wi；
(3)、获取NS个粒子构成的群体的初始全局最佳位置
第i(i＝1,2,…,NS)个粒子在寻优过程中的适应值函数为：
fit i ( W ) = 1 N 0 Σ t = 1 N 0 [ Y t - W T X t ] 2 - - - ( 1 ) ]]>
选取适应值最小的粒子的位置作为NS个粒子构成的群体的初始全局最佳
位置Fbest；
(4)、根据万有引力算法计算出每一维空间中粒子i(i＝1,2,…,NS)的加速度
因此M维空间中全体粒子的加速度可以用NS×M维矩阵a表示；
(5)、更新粒子的位置和适应值：
m b e s t = 1 M Σ i = 1 M F b e s t ( i ) ]]>
p ( K ) = a 0 · F b e s t + a · f i t ( F b e s t ) a 1 + a ]]>
b(K)＝1-0.5K/MAXITER
r(K)＝[r1(K),…,rd(K),…,rM(K)],其中 r d ( K ) = - 1 , i f δ ( K ) > 0.5 - 1 , i f δ ( K ) ≤ 0.5 ]]>
W ( K ) = p ( K ) + r ( K ) b ( K ) | m b e s t - W ( K - 1 ) | l n ( 1 u ( K ) ) ]]>
其中，δ(K),u(K)为0至1之间的随机数，a0,a1分别为0至1之间的NS维
和NS×M维随机向量，a为粒子的加速度，fit(Fbest)为全局最佳位置所对应的
适应值，即最优适应值，K为当前迭代步数，MAXITER为最大迭代步数；
随着粒子位置W(K)的更新，从而利用公式(1)可以计算得到粒子
i(i＝1,2,…,NS)的新适应值fit(Wi(K))；
(6)、更新每个粒子的局部最佳位置，如果每个粒子fit(Wi(K))＜fit(Lbesti)，
则Lbesti＝Wi(K)，否则Lbesti保持不变(i＝1,2,…,NS)；
更新全局最佳位置， F b e s t = arg m i n 1 ≤ i ≤ N S f i t ( Lbest i ) ; ]]>
(7)、迭代操作：判断迭代步数是否达到最大迭代步数MAXITER或者最优
适应值fit(Fbest)是否达到稳定；如果是，则终止迭代，得到一组最优的M维核
向量执行步骤(8)；否则，返回到步骤(4)继续执行；
(8)、计算预测值：当预测第D+1个数据，即预测数据x(D+1)时，根据步骤
(1)可知此时预测模型的输入数据帧为 X t ^ = [ x t ^ , 1 , x t ^ , 2 , ... , x t ^ , M ] T , ( t ^ = D + 1 - d × τ ) ; ]]>
将步骤(7)中的核向量代入，从而得到时刻，预测模型的一步预
测值：
y t ^ = f ( X t ^ ) = f ( x t ^ , 1 , x t ^ , 2 , ... , x t ^ , M ) = W ^ T X t ^ . ]]>
进一步的，本发明还提供了计算每一维空间中粒子的加速度的方法为：
第K次迭代时，每个粒子的惯性质量mi(K):
c i ( K ) = [ fit i ( W ) ] K - w o r s t ( K ) b e s t ( K ) - w o r s t ( K ) ]]>
m i ( K ) = c i ( K ) Σ j = 1 N S c i ( K ) ]]>
其中，ci(K)为计算粒子惯性质量的中间变量，[fiti(W)]K是第K次迭代时粒
子i的适应值，best(K)和worst(K)代表最优适应度和最差适应度；
对于Volterra核辨识而言，目的为寻求一组M维核向量并使最小，
所以在搜索函数最小值时，best(K)和worst(K)分别为：
根据万有引力公式：F＝Gm1m2/R2,得到第K次迭代时粒子i和粒子j之间
的引力为：
F i j M ( K ) = G ( K ) m i ( K ) × m j ( K ) R i j ( K ) + ϵ ( W j ( K ) - W i ( K ) ) ]]>
其中，Wi(K),Wj(K)是M维空间中粒子i和粒子j的位置，ε是一个很小的常
量，Rij(K)是第K次迭代时，粒子i和粒子j的欧式距离，即Rij(K)＝||Wi(K),Wj(K)||2；
G(K)是万有引力系数，可表示为：G0＝100,β＝20,K为当前
迭代步数，MAXITER为最大迭代步数；
所以得到每维空间中粒子i的加速度
a i M ( K ) = Σ j = 1 , j ≠ i N S h j F i j M ( K ) m i ( K ) ]]>
其中，hj(1≤j≤NS,j≠i)代表第j个随机数。
更进一步的，本发明还描述了计算预测值时，还可以一次预测多步数据，
即先预测时刻的数据，再以时刻的预测数据为基础，预测时刻的数据，并
依次类推，从而得到q(q≥1)步超前预测值为：
y t ^ + q = f ( x t ^ , q + 1 , ... , x t ^ , m , y t ^ , ... , y t ^ + q - 1 ) q ∈ { 1 , 2 , ... , M - 1 } f ( y t ^ + q - M , ... , y t ^ + q - 1 ) q ∈ { M , M + 1 , ... } . ]]>
本发明的发明目的是这样实现的：
本发明一种基于量子引力算法的时间序列预测方法，通过采用Volterra级数
展开式来构造时间序列预测模型。首先，在引用相空间分帧重构技术的基础上，
构建以Volterra级数P阶截断的输入信号矢量；其次，引入加速度作为某个变量
参数到MohadesehSoleimanpour等人提出的量子引力算法中，利用此算法有效
地训练出预测模型的Volterra核函数；最后，通过预测点时刻的输入信号矢量与
Volterra核函数的线性组合，从而求得该时刻的预测值。通过实验验证，将量子
引力算法引入到非线性系统的Volterra核函数辨识中，可以有效地提高时间序列
预测的准确度。
附图说明
图1是本发明基于量子引力算法的时间序列预测方法的一种具体实施方式
流程图；
图2是MG、RDD、Laser三种序列在寻优过程中粒子的位置；
图3是本发明在9种不同数据集下的学习曲线示意图；
图4是本发明对Laser序列进行单步预测的预测结果和实际数据的对比图；
图5是本发明对Laser序列进行50步预测的预测结果和实际数据的对比图；
图6是不同预测模型预测Laser序列时，其均方误差的三维圆柱对比图。
具体实施方式
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更
好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设
计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
实施例
图1是本发明基于量子万有引力算法的时间序列预测方法的一种具体实施
方式流程图。
在本实施例中，如图1所示，本发明基于量子引力算法的时间序列预测方
法可以划分为三个阶段：数据预处理阶段、预测模型训练阶段和预测阶段。下
面分别对各个阶段进行详细说明：
(一)数据预处理阶段
S101：设预选的数据集DS＝{x(1),x(2),…,x(D)}，其中D为数据个数；
首先，根据微分熵率的方法确定数据集DS的最佳嵌入维数d和延迟时间τ；
其次，利用Matlab中的函数windowize对数据集DS进行分帧重构，并将各
帧数据分别映射到d维特征空间，且以Volterra级数的方式p(p≥1)阶截断，从而
得到N0(N0＝D-d×τ)个数据帧{(Xt,Yt)}(t＝1,2,…,N0)，并作为预测模型的输入和目标输出数
据，其中：输入数据为Xt＝{1,x(t),x(t+τ),…,x(t+(d-1)×τ),x2(t),x(t)×x(t+τ),…,x2(t+(d-1)×τ),…}目标
输出数据为Yt＝x(t+d×τ)(t＝1,2,…,N0)；
在本实施例中，为算法的简洁性，设输入数据的维数大小为M，则输入数
据Xt＝[xt,1,xt,2,…,xt,M]T。
(二)预测模型训练阶段
S102：设有NS个粒子，粒子位置用M维向量Wi表示，粒子群位置可用矩
阵表示，i＝1,2,…,NS；随机初始化NS个粒子的初始位置Wi，，并令各粒
子局部最佳位置为Lbesti＝Wi。
S103：分别计算NS个粒子的适应值，其中，第i(i＝1,2,…,NS)个粒子在寻
优过程中的适应值函数为：
fit i ( W ) = 1 N 0 Σ i = 1 N 0 [ Y t - W T X t ] 2 - - - ( 1 ) ]]>
选取适应值最小的粒子的位置作为NS个粒子构成的群体的初始全局最佳
位置Fbest；
S104：根据万有引力算法计算出每一维空间中粒子i(i＝1,2,…,NS)的加速度
因此M维空间中全体粒子的加速度可以用NS×M维矩阵a表示。
其中，计算每一维空间中粒子的加速度，具体过程为：
第K次迭代时，每个粒子的惯性质量mi(K):
c i ( K ) = [ fit i ( W ) ] K - w o r s t ( K ) b e s t ( K ) - w o r s t ( K ) ]]>
m i ( K ) = c i ( K ) Σ j = 1 N S c i ( K ) ]]>
其中，ci(K)为计算粒子惯性质量的中间变量，[fiti(W)]K是第K次迭代时粒
子i的适应值，best(K)和worst(K)代表最优适应度和最差适应度；
对于Volterra核辨识而言，目的为寻求一组M维核向量并使最小，
所以在搜索函数最小值时，best(K)和worst(K)分别为：
b e s t ( K ) = min j ∈ { 1 , ... , N S } [ fit j ( W ) ] K , w o r s t ( K ) = m a x j ∈ { 1 , ... , N S } [ fit j ( W ) ] K ]]>
根据万有引力公式：F＝Gm1m2/R2,得到第K次迭代时粒子i和粒子j之间
的引力为：
F i j M ( K ) = G ( K ) m i ( K ) × m j ( K ) R i j ( K ) + ϵ ( W j ( K ) - W i ( K ) ) ]]>
其中，Wi(K),Wj(K)是M维空间中粒子i和粒子j的位置，ε是一个很小的常
量，Rij(K)是第K次迭代时，粒子i和粒子j的欧式距离，即Rij(K)＝||Wi(K),Wj(K)||2；
G(K)是万有引力系数，可表示为：G0＝100,β＝20,K为当前
迭代步数，MAXITER为最大迭代步数；
所以得到每维空间中粒子i的加速度
a i M ( K ) = Σ j = 1 , j ≠ i N S h j F i j M ( K ) m i ( K ) ]]>
其中，hj(1≤j≤NS,j≠i)代表第j个随机数。
S105:更新粒子的位置和适应值：
m b e s t = 1 M Σ i = 1 M F b e s t ( i ) ]]>
p ( K ) = a 0 · F b e s t + a · f i t ( F b e s t ) a 1 + a ]]>
b(K)＝1-0.5K/MAXITER
r(K)＝[r1(K),…,rd(K),…,rM(K)],其中 r d ( K ) = - 1 , i f δ ( K ) > 0.5 - 1 , i f δ ( K ) ≤ 0.5 ]]>
W ( K ) = p ( K ) + r ( K ) b ( K ) | m b e s t - W ( K - 1 ) | l n ( 1 u ( K ) ) ]]>
其中，δ(K),u(K)为0至1之间的随机数，a0,a1分别为0至1之间的NS维
和NS×M维随机向量，a为粒子的加速度，fit(Fbest)为全局最佳位置所对应的
适应值，即最优适应值；K为当前迭代步数，MAXITER为最大迭代步数；
随着粒子位置W(K)的更新，从而利用公式(1)可以计算得到粒子
i(i＝1,2,…,NS)的新适应值fit(Wi(K))。
S106:更新每个粒子的局部最佳位置，如果每个粒子fit(Wi(K))＜fit(Lbesti)，
则Lbesti＝Wi(K)，否则Lbesti保持不变(i＝1,2,…,NS)；
更新全局最佳位置， F b e s t = arg m i n 1 ≤ i ≤ N S f i t ( Lbest i ) . ]]>
S107:迭代操作：判断迭代步数是否达到最大迭代步数MAXITER或者最优适
应值fit(Fbest)是否达到稳定；如果是，则终止迭代，得到一组最优的M维核向
量执行步骤(8)；否则，返回到步骤(4)继续执行。
(三)预测阶段
S108：计算预测值：当预测第D+1个数据，即预测数据x(D+1)时，根据步骤
S101可知此时预测模型的输入数据帧为 X t ^ = [ x t ^ , 1 , x t ^ , 2 , ... , x t ^ , M ] T , ( t ^ = D + 1 - d × τ ) ; ]]>
将步骤S107中的核向量代入，从而得到时刻，预测模型的一步预
测值：
y t ^ = f ( X t ^ ) = f ( x t ^ , 1 , x t ^ , 2 , ... , x t ^ , M ) = W ^ T X t ^ ]]>
在实际中，可以一次只预测一步的数据，即时刻的数据，也可以一次预测
多步数据，即时刻的数据。如果一次预测多步数据，则先预测时刻
的数据，再以时刻的预测数据为基础，预测时刻的数据，并依次类推，从
而得到q(q≥1)步超前预测值为：
y t ^ + q = f ( x t ^ , q + 1 , ... , x t ^ , m , y t ^ , ... , y t ^ + q - 1 ) q ∈ { 1 , 2 , ... , M - 1 } f ( y t ^ + q - M , ... , y t ^ + q - 1 ) q ∈ { M , M + 1 , ... } . ]]>
实例
为了说明本发明的技术效果，本发明将分别从仿真和实际的混沌时间序列
进行实验分析，以衡量预测模型的有效性。其中仿真序列为著名的Mackey-Glass
(MG)数据，实际序列为TimeSeriesPredictionGroup小组提供的几组数据
集:Laser、DailyMinimumTemperatures(DMT)、Electricity_Demand(ED)、
CATS_Benchmark(CATS_B)、SunspotNumber(SN)、PolandElectricity(PE)、Dslp。
同时，本发明还利用随机退化模型产生的一组退化测量数据(RDD)用于验证本发
明所提出的方法。
对于这9个不同的仿真和实际混沌时间序列，首先将其进行归一化，然后
将2/3的数据作为训练数据，1/3的数据作为测试数据。在预测模型的参数设置
中，采用微分熵率的方法选取最佳嵌入维数和延迟时间；设置模型的粒子群个
数为50，最大迭代步数MAXITER为1000，以及Volterra级数的阶数为3。
在上述参数条件下，图2(a)(b)(c)分别给出了MG、RDD、Laser三种序列在
寻优过程中粒子的位置，图中五角星代表了适应值达到稳定时所对应的粒子最
佳位置，圆圈代表每次迭代时粒子的最佳位置。虽然MG序列随机性较弱，而
RDD和Laser序列随机性较强，但从图2的三幅图可以看出，本发明所提出的
方法对这两种类型的随机序列均具有较强的搜索能力。进而，图3给出了该方
法在上述9种不同数据集下的学习曲线示意图。可见，在函数寻优过程中，针
对不同的数据集，该方法都能拥有较高的收敛速度以及收敛精度。
训练结束后利用预测模型对测试数据点进行预测。图4是采用本发明对
Laser序列进行单步预测的预测结果和实际数据的对比图。图5是采用本发明对
Laser序列进行50步预测的预测结果和实际数据的对比图。如图4和图5所示，
以实线表示实际数据，虚线表示预测结果。从图4和图5可以看出，不管是单
步预测还是多步预测，本发明所提出的方法都具有良好的预测性能，即使在50
步的超长预测情况，其预测输出曲线也非常逼近于真实序列。
为更好地验证本发明所提出的方法的预测性能，这里采用均方误差(MSE)
和规范化均方根误差(NRMSE)作为性能评价标准，利用本发明所提出的预测
模型分别对上述9个不同的仿真和实际混沌时间序列进行单步和多步预测。为
了确保测试的准确性，对不同步长的预测均重复30次实验，记录每次实验的
MSE和NRMSE，最后分别求其平均值。
表1则是不同数据集下该模型预测结果的误差分析。
表1
从表1的误差分析结果可以看出，不管是随机性较弱的MG序列，还是对
于随机性很强的DMT、SN、PE、RDD等序列，本发明提出的方法对于这两种
类型的随机序列均可取得较好的预测效果。
进一步，为了定量衡量和比较预测性能，我们还利用了自回归模型(AR)，
标准的径向基神经网络模型(RBF)，最优剪枝极限学习机模型(OPELM)，
Volterra,裁剪懒惰学习模型(LLpruned)，最小二乘支持向量机(LSSVM),以
及本发明这七种不同的预测模型分别对Laser序列进行预测，所有预测模型的训
练数据和测试数据在输入模型前都进行了归一化处理。为了对结果能有直观地
认识，图6画出了不同预测模型预测Laser序列时，其均方误差的三维圆柱对比
图。从对比图可以看出，本发明所提出的模型要比目前比较流行且常用的大多
数预测模型具有更高的预测精度。
上述实验结果表明，本发明一种基于量子引力算法的时间序列预测方法，
不仅在函数寻优过程中拥有较高的收敛速度及收敛精度，而且在单步或多步超
前预测情况下，均可对不同类型的随机序列取得较好的预测效果。
尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的
技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本
技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的
本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明
创造均在保护之列。