本发明涉及一种复杂时间序列预测模型，特别涉及一种基于改进的多目标差分优化回声状态网络的时间序列预测模型。
背景技术：
时间序列广泛存在于社会生活的各个方面。它在经济领域中广泛存在，例如商品价格，股票指数。它也广泛存在于工业领域，例如轴承健康监测。在水文领域，年径流量和星系形态，等等。因此，时间序列预测的准确性通常是许多研究人员追求的目标。在过去的几十年中，已经提出了各种时间序列预测模型，例如自回归模型，神经网络，支持向量回归和模糊系统。尤其是神经网络在处理非线性方面表现出极大的优势。在许多学者的不断努力下，神经网络已经得到了深刻的发展。例如，jaeger和haas提出了一种用于训练神经网络的随机算法，称为回声状态网络(esn)。发现esn在时间序列预测中具有一定优势。它训练算法简单，训练时间短，并且保证了解的全局最优性。尽管回声状态网络具有上述优点，但也存在一些问题，例如储备池的适应性和稳定性问题，病态解的出现以及对不同时间序列esn参数值设定的问题。
近些年来，许多学者利用群智能优化算法优化神经网络，得到了很好的研究成果。例如乔玉晶等的《基于蜂群算法优化bp神经网络三目视觉标定》(公开号：cn110009696a)利用人工蜂群算法为bp神经网络选取最佳的初始权重与偏置。但是人工蜂群算法对具体问题和应用环境的依赖性较大，且计算复杂度高。bp神经网络的训练时间较长，整体来看蜂群算法优化bp神经网络这个算法的适用性不高。玄东吉等的《基于蚁群算法优化后的bp神经网络来估算soc的方法》(公开号：cn109738807a)选取蚁群算法来优化bp神经网络来估算现阶段动力电池的soc。尽管蚁群算法具有很强的鲁棒性，但它具有明显的经验性，收敛速度慢和适用性低的特点，因此它不适合直接用于实际问题。易辉等的《基于改进的萤火虫算法优化概率神经网络的电力变压器故障诊断方法及系统》(公开号：cn110363277a)改进的萤火虫算法来优化概率神经网络进行电力变压器故障诊断，达到诊断精度高，误差小的效果。但是萤火虫算法发现率低、求解精读不高、求解速度慢，导致萤火虫算法优化概率神经网络的开销很大。
目前熟知的群智能优化算法有一定的不足。遗传算法的局部搜索能力不强，无法找到全局最优解，搜索效率低。粒子群算法的鲁棒性很差，而且粒子群算法容易陷入局部最优解，参数的初始化对算法的性能影响很大。人工鱼群算法尽管不易陷入局部最优，但不适用于同时优化具有不同范围的参数，其本身算法结构复杂且效率低下。鸽群优化算法收敛速度慢且稳定性相对较差。教与学优化算法具有过早收敛的现象。研究发现，差分进化算法结构很简单而且具有强大的全局搜索能力，并且可以平衡局部和全局信息以进行搜索。因此本发明选择差分进化来优化回声状态网络的储备池参数，为了进一步提高预测的准确定和稳定性，本发明选择多目标差分来进行两个目标函数的优化。
技术实现要素：
本发明的主要目的是解决回声状态网络的储备池参数应对不同时间序列时难以确定的问题，提供一种利用改进的多目标差分进化算法来优化回声状态网络的储备池参数的预测模型，不必人工手动调节参数以此节省时间，而且适用于不同的时间序列预测，提高预测精度。
为了达到上述目的，本发明采用的技术方案为：
回声状态网络由输入层、储备池、输出层组成，储备池参数的设定直接决定着回声状态网络的预测性能。人工调节既费时间又不能选出最合适的参数值，因此本发明利用改进的多目标差分优化回声状态网络的储备池参数，构造一种基于改进多目标差分算法优化回声状态网络的时间序列预测模型。待优化的参数是回声状态网络的储备池参数，分别是储备池规模、谱半径、稀疏度、输入变换因子。本发明使用的差分进化结构简单、易于实现、鲁棒性高、全局搜索能力强。差分进化算法主要分为四个操作：初始化、变异操作、交叉操作和选择操作。单目标优化的目标函数是预测值的绝对误差最小，缺点是不能保证预测效果的稳定性，本发明对改进的多目标差分进化算法采用两个目标函数，分别为f1(x)和f2(x)，其中，第一个函数表示绝对误差的最小和，以确保预测准确性。另一个函数表示误差的方差最小，以保证预测的稳定性。目标函数如公式(1)所示：
在这里y(t)为数据的真实值，y′(t)是回声状态网络的输出值即预测值。n为样本数据的数量。
其步骤如下：
步骤1：随机初始化种群，np为种群个体数，g为当前迭代次数。种群个体的每一维代表回声状态网络中储备池的一个参数,每一个参数都有选择范围，因此对个体的每一维进行约束。
步骤2：种群中的个体用xi表示，其中i＝1,2,...,np。将个体xi的每一维分别赋值给回声状态网络中储备池对应的参数：储备池规模、稀疏度、谱半径和输入变换因子。实验数据被划分为训练集和测试集。训练集数据输入到回声状态网络中，得出预测值并进行该个体的适应度评价即计算f(x)。依次对种群中的个体进行适应度评价。
步骤3：设置种群最大的迭代次数为max-iteration，如果迭代次数大于设定的最大迭代次数，则停止迭代，执行步骤8；反之，执行步骤4。
步骤4：采用差分进化的变异操作与交叉操作为种群中的每个个体生成变异个体和试验个体，将生成的试验个体组成试验种群，记作trial-pop。此时将当代种群和试验种群组成一个混合种群，记作mix-pop。所述的变异操作和交叉操作中用到的缩放因子f及交叉概率cr均为0～1之间的任意值，且不能取0、1。
步骤5：利用一组随机均匀分布的参考向量l＝{l1,l2,...,lk},k＝1,2,...,k将混合种群分解成k个子种群。种群分解原则是为混合种群中每个个体找到离其本身最近的参考向量，并将每个个体划入离它最近的向量所对应的子种群中。种群个体xi在向量lk上的投影距离记作d1，垂直距离记作d2。d1和d2分别代表种群个体的收敛性和多样性。
种群分解原则具体为：分别计算混合种群中每个种群个体与k个参考向量之间的夹角，找到最小夹角，最小夹角表示为αi＝min{arg(cos(xi,lk))}，并将该种群个体xi划分到向量lk对应的子种群subpop{k}中。
步骤6：将每个子种群根据
子种群更新规则如下：
①如果|pcon|>np/k，对子集pcon按照ψ(x)＝d1(x) θ×d2(x)进行升序排列，选取前np/k个解，组成下一代子种群。其中，θ为惩罚函数，|pcon|代表子集的大小。np表示种群的大小，k表示子种群个数。
②如果|pcon|＝np/k，本发明将子集pcon中的解作为下一代子种群。
③如果|pcon|<np/k，此时子集pcon中的解全进入到下一代子种群，并将子集pdiv按照d1值进行升序排列，之后在子集pdiv中选取前(np/k-|pcon|)个解进入下一代子种群。
步骤7：所有更新完的子种群即为下一代新的种群，此时g＝g 1，执行步骤3-步骤6。
步骤8：当迭代次数达到最大，得到的一个最终种群，从种群中随机选出一个个体作为最优的输出个体，该输出个体的每一维赋值给回声状态网络的储备池参数，并将测试集数据输入到该回声状态网络，输出预测值。根据测试集数据的真实值和网络的预测值绘制预测曲线和误差曲线。
综上，步骤1至步骤8为本发明提出的基于改进的多目标差分优化回声状态网络的时间序列预测模型的步骤描述。
本发明的特点与益处在于：多目标差分算法结构简单、易于实现、鲁棒性高、全局搜索能力强，同时多目标优化使得预测效果同时达到准确性和稳定性。改进的多目标差分算法得到的最终种群具有很好的收敛性和多样性，也为不同的时间序列选择合适的回声状态网络储备池参数，提高预测精度。
附图说明
图1为本文提出模型的步骤流程图。
图2为多目标差分进化算法的种群分解示意图。
图3为多目标差分进化算法的子种群分解示意图。
图4为实施例一的大连月平均气温:imode-esn的预测曲线及误差曲线仿真结果图；图4(a)为预测曲线图，图4(b)为误差曲线图，图4(c)为局部的预测曲线，图4(d)为局部的误差曲线。
图5为实施例二的北京pm2.5时间序列:imode-esn的预测曲线及误差曲线仿真结果图；图5(a)为预测曲线图，图5(b)为误差曲线图，图5(c)为局部的预测曲线，图5(d)为局部的误差曲线。
具体实施方式
下面结合附图对本发明的技术方案进行详细说明：
回声状态网络(echostatenetwork,esn)优点在于网络学习过程中仅需调节储备池到输出层之间的输出连接权值，而其他连接权值一般都随机赋值后保持不变。储备池参数设置直接影响esn的预测性能，人工调节参数既费时又不能选择出最适合的参数值，本发明基于改进的多目标差分进化算法优化esn，简记imode-esn模型。该模型实现esn参数的自适应，要优化的参数是储备池规模、储备池的谱半径、稀疏度、输入变换因子。
差分进化算法(differentialevolution,de)有结构简单、易于实现、鲁棒性高、全局搜索能力强等特点。差分进化算法主要分为四个操作：初始化、变异操作、交叉操作和选择操作。同时为了平衡预测的准确性和稳定性，本发明使用改进的多目标差分来优化储备池参数，以弥补单目标优化算法的缺点。因此，本发明使用两个目标函数，分别为f1(x)和f2(x)，其中，第一个函数表示绝对误差的最小和，以确保预测准确性。另一个函数表示误差的方差最小，以保证预测的稳定性。
本发明基于改进的多目标差分优化回声状态网络的时间序列预测模型，如图1所示，按照以下步骤进行：
步骤1：随机初始化种群，种群大小np＝100，g为当前迭代次数。种群个体的每一维代表回声状态网络中储备池的一个参数,每一个参数都有选择范围，因此对个体的每一维进行约束。
步骤2：种群中的个体用xi表示，其中i＝1,2,...,50。将个体xi的每一维分别赋值给回声状态网络中储备池对应的参数：储备池规模、稀疏度、谱半径和输入变换因子。实验数据被划分为训练集和测试集。训练集数据输入到回声状态网络中，得出预测值并进行该个体的适应度评价即计算f(x)。依次对种群中的个体进行适应度评价。
步骤3：设置种群最大的迭代次数max-iteration＝30，如果迭代次数大于设定的最大迭代次数，则停止迭代，执行步骤8；反之，执行步骤4。
步骤4：采用差分进化的变异操作与交叉操作为种群中的每个个体生成变异个体和试验个体，将生成的试验个体组成试验种群，记作trial-pop。此时将当代种群和试验种群组成一个混合种群，记作mix-pop。所述的变异操作和交叉操作中用到的缩放因子f＝0.9交叉概率cr＝0.3。
步骤5：如图2所示，利用一组随机均匀分布的参考向量l＝{l1,l2,...,l20},k＝1,2,...,20将混合种群分解成20个子种群，。种群分解原则是为混合种群中每个个体找到离其本身最近的参考向量，并将每个个体划入离它最近的向量所对应的子种群中。种群个体xi在向量lk上的投影距离记作d1，垂直距离记作d2，公式如(2)、(3)所示。d1和d2分别代表种群个体的收敛性和多样性。
d1(x)＝ltf*(x)(2)
其中，l表示随机均匀分布的参考向量，f(x)min＝{min(f1(x)),min(f2(x))}表示一个理想点。f*(x)＝f(x)-f(x)min。
种群分解原则具体为：分别计算混合种群中每个种群个体与20个参考向量之间的夹角，找到最小夹角，最小夹角表示为αi＝min{arg(cos(xi,lk))}，并将该种群个体xi划分到向量lk对应的子种群subpop{k}中。
步骤6：如图3所示，将每个子种群根据
子种群更新规则如下：
①如果|pcon|>10，对子集pcon按照ψ(x)＝d1(x) θ×d2(x)进行升序排列，选取前10个解，组成下一代子种群。其中，θ为惩罚函数，|pcon|代表子集的大小。np表示种群的大小，np＝100，k表示子种群个数，k＝20。
②如果|pcon|＝10，本发明将子集pcon中的解作为下一代子种群。
③如果|pcon|<10，此时子集pcon中的解全进入到下一代子种群，并将子集pdiv按照d1值进行升序排列，之后在子集pdiv中选取前(10-|pcon|)个解进入下一代子种群。
步骤7：所有更新完的子种群即为下一代新的种群，此时g＝g 1，执行步骤3-步骤6。
步骤8：当迭代次数达到最大，得到的一个最终种群，从种群中随机选出一个个体作为最优的输出个体，该输出个体的每一维赋值给回声状态网络的储备池参数，并将测试集数据输入到该回声状态网络，输出预测值。根据测试集数据的真实值和网络的预测值绘制预测曲线和误差曲线。
本发明选择了大连月平均气温时间序列和北京pm2.5时间序列进行仿真实验。同时，本发明在相同条件下选择四种优化算法优化esn来作为对比模型。它们分别是人工鱼群算法优化esn(afsa-esn)、粒子群算法优化esn(pso-esn)、教与学优化算法优化esn(tlbo-esn)、基于指标的进化算法(ibea)优化esn(ibea-esn)。预测性能的评估指标是均方根误差(rmse)，对称均方根绝对百分比误差(smape)和归一化均方根误差(nrmse)。rmse，smape和nrmse的计算公式定义如下：
在这里y(t)为数据的真实值，y′(t)是网络预测值。
本发明给出了两组实验的待优化参数值的范围。然后通过优化算法找到给定范围内最合适的参数值。储备池参数范围设置如表1所示。
表1储备池参数范围
实施例一：
实施例一通过大连市月平均气温和月降雨量数据集来说明本发明的有益效果。选取大连地区1951年1月至2017年12月的月平均温度降水数据集，共792个样本。采样间隔为月。大连月降水量用作实验数据的第二维，以辅助大连月平均气温的预测。75％的数据用于训练，25％的数据用于测试。在训练样本中，前50个样本的状态被丢弃，以清除初始暂态的影响。为了证明本发明的有效性，对于pso-esn，afsa-esn，tlbo-esn，ibea-esn，本发明使esn的储备池参数的范围与表1所示相同。种群大小和最大迭代次数也相同。表2给出了大连的月平均气温时间序列选择最合适的储备池参数值。
表2大连月平均气温:imode-esn模型参数
表3给出了不同模型对大连市月平均气温时间序列进行单步预测的结果和运行时间。在rmse，smape和nrmse方面，imode-esn优于其他对比模型。从表3可以看出，pso-esn的smape很大，说明pso-esn无法预测大连月平均气温的序列趋势。尽管afsa-esn比pso-esn更好的预测序列趋势，但rmse和nrmse的误差很大，这表明在某些点上的预测存在较大偏差。在rmse，smape和nrmse方面，imode-esn小于其他模型。表3也给出了大连月平均气温时间序列在不同模型上的运行时间。可以看出，afsa-esn的运行时间最长。imode-esn的运行时间比pso-esn和ibea-esn略长，tlbo-esn的运行时间是imode-esn的四倍。尽管ibea-esn的运行时间最短，但预测精度却比tlbo-esn和imode-esn差。可以看出，imode-esn模型在rmse上的优势是显而易见的。图4(a)给出了imode-esn对大连月平均气温的预测曲线，图4(b)给出了imode-esn对大连月平均气温的误差曲线。从图4(a)中可以看出，imode-esn能很好地拟合月平均气温的真实值曲线，从图4(b)中可以看出真实值与预测值的误差很小。因为图4(a)中的真实值曲线与预测值曲线拟合的很好，为增加可观性，图4(c)给出了图4(a)中黑框标出的曲线的放大图。图4(d)为对应的放大的误差曲线。
表3大连月平均气温:测试集仿真结果
实施例二：
较佳地，实施例二采用来自北京空气污染指数(aqi)的北京pm2.5时间序列。本发明使用收敛交叉映射(ccm)为北京pm2.5时间序列选择辅助变量，发现pm10，co和so2对pm2.5的影响大。因此，本发明选择pm10，co和so2作为pm2.5时间序列预测的辅助变量。本发明选取从2016年1月2日到2016年12月31日，共有8759个样本用于仿真。总样本的70％用于训练，总样本的30％用于测试。在训练样本中，前100个样本的状态被丢弃以消除初始暂态的影响。
从表1中可以知道储备池参数的范围设置。为了证明本发明提出模型的有效性，本发明对afsa-esn，pso-esn，tlbo-esn，ibea-esn模型进行了相同的设置。表4给出了本发明为北京pm2.5时间序列选择的最合适的储备池参数。从表4可知，储备池规模为22，略大于储备池范围的下限。对于esn，当前的储备池结构简单，运行速度更快。
表4北京pm2.5时间序列:imode-esn模型参数
表5给出了不同模型对北京pm2.5时间序列的单步预测结果。可以看出，就rmse，smape和nrmse而言，imode-esn仍然优于其他比较模型。可以看出，对于afsa-esn，tlbo-esn，ibea-esn，rmse很大，这表明它们不能很好地预测北京pm2.5时间序列。图5(a)和图5(a)分别显示了imode-esn根据测试数据得出的北京pm2.5时间序列的预测曲线和误差曲线。imode-esn可以很好地预测北京pm2.5时间序列的变化趋势。因为图5(a)中的真实值曲线与预测值曲线拟合的很好，为增加可观性，图5(c)给出了图5(a)中1000步至1020步的曲线的放大图。图5(d)为对应的放大的误差曲线。
表5显示了不同模型在20次迭代后的运行时间。可以看出，afsa-esn仍具有最长的运行时间。imode-esn的运行时间比ibea-esn的运行时间稍长。尽管ibea-esn的运行时间最短，但预测精度却比imode-esn差。
表5北京pm2.5时间序列:测试集仿真结果
以上仿真结果表明，imode-esn模型在大连月平均气温时间序列和北京pm2.5时间序列上具有出色的表现。可以得出本发明可以选择合适的储备池参数设置以适用于不同的时间序列，从而提高了模型的泛化性能。从预测精度方面看，该模型预测出来的rmse最小；从运行时间的角度来看，该模型的计算复杂度相对较低。因此本发明提出的模型在预测精度和计算复杂度方面具有良好的前景。
上述说明书虽然结合附图详细描述了本发明的具体实施方式，但本发明并不局限于上述实施例的描述范围。本领域技术人员应该知道，在本专利的实质范围内，进行的各种修改和替换，都应属于本专利的保护范围。
1.一种基于改进多目标差分算法优化回声状态网络的时间序列预测模型，其特征在于，所述改进的多目标差分进化算法采用两个目标函数，分别为f1(x)和f2(x)，其中，第一个函数表示绝对误差的最小和，以确保预测准确性；另一个函数表示误差的方差最小，以保证预测的稳定性；目标函数如公式(1)所示：
在这里y(t)为数据的真实值，y′(t)是回声状态网络的输出值即预测值；n为样本数据的数量；
其步骤如下：
步骤1：随机初始化种群，np为种群个体数，g为当前迭代次数；种群个体的每一维代表回声状态网络中储备池的一个参数，每一个参数都有选择范围，因此对个体的每一维进行约束；
步骤2：种群中的个体用xi表示，其中i＝1,2,...,np；将个体xi的每一维分别赋值给回声状态网络中储备池对应的参数：储备池规模、稀疏度、谱半径和输入变换因子；将实验数据划分为训练集和测试集；训练集数据输入到回声状态网络中，得出预测值并进行该个体的适应度评价即计算f(x)；依次对种群中的个体进行适应度评价；
步骤3：设置种群最大的迭代次数为max-iteration，如果迭代次数大于设定的最大迭代次数，则停止迭代，执行步骤8；反之，执行步骤4；
步骤4：采用差分进化的变异操作与交叉操作为种群中的每个个体生成变异个体和试验个体，将生成的试验个体组成试验种群，记作trial-pop；此时将当代种群和试验种群组成一个混合种群，记作mix-pop；
步骤5：利用一组随机均匀分布的参考向量l＝{l1,l2,...,lk},k＝1,2,...,k将混合种群分解成k个子种群；种群分解原则是为混合种群中每个个体找到离其本身最近的参考向量，并将每个个体划入离它最近的向量所对应的子种群中；种群个体xi在向量lk上的投影距离记作d1，垂直距离记作d2；
步骤6：将每个子种群根据
子种群更新规则如下：
①如果|pcon|>np/k，对子集pcon按照ψ(x)＝d1(x) θ×d2(x)进行升序排列，选取前np/k个解，组成下一代子种群；其中，θ为惩罚函数，|pcon|代表子集的大小；np表示种群的大小，k表示子种群个数；
②如果|pcon|＝np/k，本发明将子集pcon中的解作为下一代子种群；
③如果|pcon|<np/k，此时子集pcon中的解全进入到下一代子种群，并将子集pdiv按照d1值进行升序排列，之后在子集pdiv中选取前(np/k-|pcon|)个解进入下一代子种群；
步骤7：所有更新完的子种群即为下一代新的种群，此时g＝g 1，执行步骤3-步骤6；
步骤8：当迭代次数达到最大，得到的一个最终种群，从种群中随机选出一个个体作为最优的输出个体，该输出个体的每一维赋值给回声状态网络的储备池参数，并将测试集数据输入到该回声状态网络，输出预测值；根据测试集数据的真实值和网络的预测值绘制预测曲线和误差曲线。
2.根据权利要求1所述的一种基于改进多目标差分算法优化回声状态网络的时间序列预测模型，其特征在于，步骤4中所述的变异操作和交叉操作中用到的缩放因子f及交叉概率cr均为0～1之间的任意值，且不能取0、1。
3.根据权利要求1所述的一种基于改进多目标差分算法优化回声状态网络的时间序列预测模型，其特征在于，步骤5所述的种群分解原则具体为：分别计算混合种群中每个种群个体与k个参考向量之间的夹角，找到最小夹角，最小夹角表示为αi＝min{arg(cos(xi,lk))}，并将该种群个体xi划分到向量lk对应的子种群subpop{k}中。
技术总结