CN108898257A - 一种基于粒子群优化的灰色预测方法 - Google Patents

Abstract

一种基于粒子群优化的灰色预测方法，首先构造改进的灰色非线性模型，然后利用粒子群优化算法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in‑1,N_i,ξ_i]进行寻找最优解。本发明不仅提高了预测准确度，而且显著提升了预测效率。

Description

一种基于粒子群优化的灰色预测方法

技术领域

本发明涉及一种基于粒子群优化的灰色预测方法。

背景技术

灰色预测在已知较少信息的情况预测对象未来发展趋势方面上具有十分突出的效果。在泛滥的不确定随机系统信息面前，有效地提取出、筛选出更多有价值的信息，过滤价值少的信息部分，成为了人们日常活动的任务和目标之一。灰色理论充分有效地利用已知信息来探究系统内部变化规律，通过对部分少量信息的累加、累减等生成处理，能够获得有价值的部分信息，进而得知系统运行变化趋势等规律的有效预测和合理监控。灰色预测模型以GM(1,1)模型为核心，能根据已知的少量信息进行建模并预测，目前已经拓展到诸多的实际预测模型中。

经典灰色系统GM(1,1)模型的背景值构造公式为：

式中：x⁽⁰⁾(i)为原始序列。如图1所示，传统背景值等于梯形的面积，但是曲线x⁽¹⁾(t)在[k-1,k]与t轴围成的面积才是真正意义上的背景值，而两者之差，即图中曲线与直线围成的面积部分，就是传统背景值导致的误差。当最初的预测序列变动幅度不大时，取紧邻均值计算得到的背景值是符合的，预测模型的误差也不会太大。然而在实际情况中，尽管是在较短的时间隔间Δt里，但是预测序列的变化速度不是以常见的平均速率产生的，所以会无法避免的产生突变量。换言之，如果原始数据存在较大的不稳定性时应用公式(1)计算得到的背景值容易存在很大的误差。

发明内容

本发明提供一种基于粒子群优化的灰色预测方法，不仅提高了预测准确度，而且显著提升了预测效率。

为了达到上述目的，本发明提供一种基于粒子群优化的灰色预测方法，包含以下步骤：

步骤S1、构造改进的灰色非线性模型；

步骤S2、利用粒子群优化算法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解。

所述的步骤S1中，构造改进的灰色非线性模型的方法包含以下步骤：

步骤S1.1、灰色非线性预测模型的背景构造公式为：

z⁽¹⁾(k)＝(1-p)x⁽¹⁾(k)+px⁽¹⁾(k-1),k＝2,3,…,n (2)

设X⁽⁰⁾＝[x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)]为非负序列，X⁽¹⁾为累加序列若X⁽¹⁾可以用光滑的曲线F(t,a,u,N)拟合，其中N＞0为待定常数，并且F(t,a,u,N)满足非线性微分方程：

步骤S1.2、对一阶累加生成序列建立如下方程：

x⁰(k)+a(z¹(k))^N＝u,k＝2,3,…,n (4)

式中，

式(4)简记为GM(1,1，N)模型，其白化微分形式为：

步骤S1.3、求解灰参数估值为：

记参数B和Y分别为：

步骤S1.4、权值p和参数a存在如下关系：

在不同的时刻引入不同的p_i改进背景值的计算公式：

z⁽¹⁾(k)＝p_k-1x⁽¹⁾(k-1)+(1-p_k-1)x⁽¹⁾(k),k＝2,3,…,n (8)

灰色参数为：

其中，

步骤S1.5、对模型的初始值进行改正，其改进的白化形式为：

步骤S1.6、四阶龙格库塔法经典表述如下：

y′＝f(t,y)，y(t₀)＝y₀ (11)

其中：

k₁＝f(t_n,y_n)

k₄＝f(t_n+h,y_n+hk₃)

式中，f代表需逼近的函数，待求值y_n+1由y_n加上时间间隔h与一个估算的斜率乘积决定，该斜率由以下斜率加权平均：k₁是时间开始时的斜率；k₂是时间中点的斜率，通过采用斜率k₁决定y在的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y的值；k₄是时间终点的斜率，其y值由k₃决定，当四个斜率取平均时，给中点的斜率赋较大的权值。对给定修正的初始值x⁽⁰⁾(1)＝x⁽⁰⁾(1)+ξ代替y₀作为四阶龙格库塔法的初值，利用上式(12)递推预测。

所述的步骤S2中，利用粒子群优化法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解的方法包含以下步骤：

步骤S2.0、构造适应度函数；

适应度函数选择平均相对误差均值最小，即：

步骤S2.1、初始化粒子群；

种群取值大小M、粒子维数n+1、最大遍历次数i_max、粒子Y轴坐标y_i＝[p_i1,p_i2，…,p_in-1,N_i,ξ_i]和速度大小v_i＝[v_i1,v_i2，…,v_in-1,v_i,n,v_i,n+1]，0.1≤p_i≤1，N∈(0,2)；

步骤S2.2、根据适应度函数计算每个粒子的适应度f_i；

步骤S2.3、根据所有的粒子适应度大小选择各自的极值p_i和全局极值p_g，遍历寻找最优解，优化惯性权重，优化所有粒子的速度和坐标，同时控制其速度符合v∈(-2,2)；

设搜索的全局为D维空间，种群中共有M个粒子，粒子i在第t次遍历的X轴位置为粒子的速度为在每次迭代中粒子i跟踪个体极值点p_i＝[p_i1,p_i2,…,p_iD]^T和全局极值点p_g＝[p_g1,p_g2,…,p_gD]^T来优化自己的位置和速度，原则为：

为了满足粒子的速度在迭代时不会失控，令v∈(-2,2)，公式(13)的为粒子原来的速度；为认知方面，表示粒子原有的信息；为群体信息方面，表示各个粒子之间的彼此分享和信息传递；在式(13)中，t为目前遍历次数；c₁为粒子优化本体的最佳状态的学习因子；c₂为粒子优化所有粒子最佳状态的学习因子；r₁和r₂为[0,1]范围的待定值，ω为惯性权重，随迭代次数线性递减的ω值为：

ω_max＝0.95，ω_min＝0.4，i为目前遍历次数，i_max为最大遍历次数；

步骤S2.4、将各个粒子的适应度f_i同目前的最佳位置f(P_i)作对比，如果f_i＜f(P_i)，则让f_i作为当前粒子的最佳值，f(P_i)＝f_i，同时优化当前粒子的最佳位置P_i＝y_i；

步骤S2.5、把所有粒子的目前的最佳状态和全局最佳状态作对比，如果f(P_i)＜f(P_g)，则将f(P_i)当作目前整体的最佳状态，即f(P_i)＝f(P_g)，同时更新整体的最佳位置V_i＝V_i+c₁×rand()×(pbest_i-X_i)+c₂×rand()×(gbest_i-X_i)，X_i＝X_i+V_i，V_i为粒子的速度，pbest_i，gbest_i为粒子的个体最优值和整体最优值，rand()为介于0～1之间的随机数，X_i为粒子的当前速度，c₁和c₂为学习因子；

步骤S2.6、判断遍历次数是否大于等于最大遍历次数，如果是，则输出全局极值f(P_i)和整体最佳位置；如果否，则进行步骤S2.2，一直到符合遍历要求为止。

本发明通过粒子群优化法对灰色模型的背景值求解方法进行改进，不仅提高了预测准确度，而且显著提升了预测效率。

附图说明

图1是背景技术中传统背景值误差来源示意图。

图2是粒子群优化算法的流程图。

具体实施方式

以下根据图2具体说明本发明的较佳实施例。

本发明提供一种基于粒子群优化的灰色预测方法，包含以下步骤：

步骤S1、构造改进的灰色非线性模型；

步骤S2、利用粒子群优化(PSO，Particle Swarm Optimization)算法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解。

进一步，所述的步骤S1中，构造改进的灰色非线性模型的方法包含以下步骤：

步骤S1.1、灰色非线性预测模型的背景构造公式为：

z⁽¹⁾(k)＝(1-p)x⁽¹⁾(k)+px⁽¹⁾(k-1),k＝2,3,…,n (2)

设X⁽⁰⁾＝[x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)]为非负序列，X⁽¹⁾为累加序列若X⁽¹⁾可以用光滑的曲线F(t,a,u,N)拟合，其中N＞0为待定常数，并且F(t,a,u,N)满足非线性微分方程：

步骤S1.2、对一阶累加生成序列建立如下方程：

x⁰(k)+a(z¹(k))^N＝u,k＝2,3,…,n (4)

式中，

式(4)简记为GM(1,1，N)模型，其白化微分形式为：

步骤S1.3、求解灰参数估值为：

记参数B和Y分别为：

灰色非线性模型反映了拟合曲线的凹凸性和原始序列变化率之间的关系，因此更能充分利用原始序列的信息，当N＝1时灰色非线性模型转化为传统GM(1,1)模型；

步骤S1.4、从建模过程可知，灰色非线性模型仍是以紧邻均值作为背景值求解参数，但是这种取值可能会带来较大的误差，因此引入加权的背景值来求解参数，由前面公式推导可知权值p和参数a存在如下关系：

传统GM(1,1)模型背景值的构造是取一个统一的p值，即p₁＝p₂＝…＝p_n-1＝0.5，这样取得的p_i(i＝1,2,…,n-1)是导致|a|较大时模型预测失效的重要原因，也就是说p_i的取值直接影响预测模型的精度，因此在不同的时刻引入不同的p_i改进背景值的计算公式：

z⁽¹⁾(k)＝p_k-1x⁽¹⁾(k-1)+(1-p_k-1)x⁽¹⁾(k),k＝2,3,…,n (8)

灰色参数为：

其中，

步骤S1.5、对模型的初始值进行改正，其改进的白化形式为：

为求得灰参数需要先确定p₁,p₂,…,p_n-1,N_i，为使得建立的模型用于预测需要利用四阶龙格库塔法以x⁽⁰⁾(1)＝x⁽⁰⁾(1)+ξ，(ξ为待求参数)为修正的初始值进行递推预报累加序列的模拟值和预测值，再通过累减还原成原始序列的模拟值和预测值；

步骤S1.6、从(10)中第一式可知，很难用一个具体的数学表达式表示此微分方程的通解，而四阶龙格库塔法却能较好的逼近此类方程，该算法是构建在严密的数学基础上的，其实现原理较为复杂，常用于数值求解微分方程，四阶龙格库塔法经典表述如下：

y′＝f(t,y)，y(t₀)＝y₀ (11)

其中：

k₁＝f(t_n,y_n)

k₄＝f(t_n+h,y_n+hk₃)

式中，f代表需逼近的函数，待求值y_n+1由y_n加上时间间隔h与一个估算的斜率乘积决定，该斜率由以下斜率加权平均：k₁是时间开始时的斜率；k₂是时间中点的斜率，通过采用斜率k₁决定y在的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y的值；k₄是时间终点的斜率，其y值由k₃决定，当四个斜率取平均时，给中点的斜率赋较大的权值。对给定修正的初始值x⁽⁰⁾(1)＝x⁽⁰⁾(1)+ξ代替y₀作为四阶龙格库塔法的初值，利用上式(12)递推预测。

由上述可知，改进的灰色非线性模型需要使用的参变量太多，常用的方案无法满足要求。但是粒子群优化算法对求解非线性、不可微的复杂优化问题有着收敛速度快、解质量高，并且粒子群优化算法原理简单，便于通过程序实现，适用在求解多目标优化问题上。因此，利用粒子群优化算法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解。

进一步，所述的步骤S2中，如图2所示，利用粒子群优化法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解的方法包含以下步骤：

步骤S2.0、为了寻找到最佳位置和整体的最佳状态，构造适应度函数；

适应度函数选择平均相对误差均值最小，即：

步骤S2.1、初始化粒子群；

种群取值大小：M＝30、粒子维数n+1、最大遍历次数i_max＝1000、粒子Y轴坐标y_i＝[p_i1,p_i2，…,p_in-1,N_i,ξ_i]和速度大小

v_i＝[v_i1,v_i2，…,v_in-1,v_i,n,v_i,n+1]，0.1≤p_i≤1，N∈(0,2)；

步骤S2.2、根据适应度函数计算每个粒子的适应度f_i；

步骤S2.3、根据所有的粒子适应度大小选择各自的极值p_i和全局极值p_g，遍历寻找最优解，优化惯性权重，优化所有粒子的速度和坐标，同时控制其速度符合v∈(-2,2)；

设搜索的全局为D维空间，种群中共有M个粒子，粒子i在第t次遍历的X轴位置为粒子的速度为在每次迭代中粒子i跟踪个体极值点p_i＝[p_i1,p_i2,…,p_iD]^T和全局极值点p_g＝[p_g1,p_g2,…,p_gD]^T来优化自己的位置和速度，原则为：

为了满足粒子的速度在迭代时不会失控，令v∈(-2,2)，公式(13)的为粒子原来的速度；为认知方面，表示粒子原有的信息；为群体信息方面，表示各个粒子之间的彼此分享和信息传递；在式(13)中，t为目前遍历次数；c₁为粒子优化本体的最佳状态的学习因子；c₂为粒子优化所有粒子最佳状态的学习因子，c₁，c₂通常会取值为2；r₁和r₂为[0,1]范围的待定值，ω为惯性权重，它是保持先前速度快慢的因子，如果ω取值较大时，容易摆脱部分最优解的缺陷，可以用来整体搜索，反之ω取值较小时，就很适合在部分范围内进行寻找最优解，从而可以加快算法收敛速度；本发明所使用的是随迭代次数线性递减的ω值为：

ω_max＝0.95，ω_min＝0.4，i为目前遍历次数，i_max为最大遍历次数；

步骤S2.4、将各个粒子的适应度f_i同目前的最佳位置f(P_i)作对比，如果f_i＜f(P_i)，则让f_i作为当前粒子的最佳值，f(P_i)＝f_i，同时优化当前粒子的最佳位置P_i＝y_i；

步骤S2.5、把所有粒子的目前的最佳状态(该次遍历过程中的最佳状态)和全局最佳状态(在本次遍历之前的最佳全局状态)作对比，如果f(P_i)＜f(P_g)，则将f(P_i)当作目前整体的最佳状态，即f(P_i)＝f(P_g)，同时更新整体的最佳位置V_i＝V_i+c₁×rand()×(pbest_i-X_i)+c₂×rand()×(gbest_i-X_i)，X_i＝X_i+V_i(V_i为粒子的速度，pbest_i，gbest_i为粒子的个体最优值和整体最优值，rand()为介于0～1之间的随机数，X_i为粒子的当前速度，c₁和c₂为学习因子)；

步骤S2.6、判断遍历次数是否大于等于最大遍历次数，如果是，则输出全局极值f(P_i)和整体最佳位置；如果否，则进行步骤S2.2，一直到符合遍历要求为止。

本发明通过粒子群优化法对灰色模型的背景值求解方法进行改进，不仅提高了预测准确度，而且显著提升了预测效率。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (3)

1.一种基于粒子群优化的灰色预测方法，其特征在于，包含以下步骤：

步骤S1、构造改进的灰色非线性模型；

步骤S2、利用粒子群优化算法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解。

2.如权利要求1所述的基于粒子群优化的灰色预测方法，其特征在于，所述的步骤S1中，构造改进的灰色非线性模型的方法包含以下步骤：

步骤S1.1、灰色非线性预测模型的背景构造公式为：

z⁽¹⁾(k)＝(1-p)x⁽¹⁾(k)+px⁽¹⁾(k-1),k＝2,3,…,n (2)

设X⁽⁰⁾＝[x⁽⁰⁾(1),x⁽⁰⁾(2),…,x⁽⁰⁾(n)]为非负序列，X⁽¹⁾为累加序列若X⁽¹⁾可以用光滑的曲线F(t,a,u,N)拟合，其中N＞0为待定常数，并且F(t,a,u,N)满足非线性微分方程：

步骤S1.2、对一阶累加生成序列建立如下方程：

x⁰(k)+a(z¹(k))^N＝u,k＝2,3,…,n (4)

式中，

式(4)简记为GM(1,1，N)模型，其白化微分形式为：

步骤S1.3、求解灰参数估值为：

记参数B和Y分别为：

步骤S1.4、权值p和参数a存在如下关系：

在不同的时刻引入不同的p_i改进背景值的计算公式：

z⁽¹⁾(k)＝p_k-1x⁽¹⁾(k-1)+(1-p_k-1)x⁽¹⁾(k),k＝2,3,…,n (8)

灰色参数为：

其中，

步骤S1.5、对模型的初始值进行改正，其改进的白化形式为：

步骤S1.6、四阶龙格库塔法经典表述如下：

y′＝f(t,y)，y(t₀)＝y₀ (11)

其中：

k₁＝f(t_n,y_n)

k₄＝f(t_n+h,y_n+hk₃)

式中，f代表需逼近的函数，待求值y_n+1由y_n加上时间间隔h与一个估算的斜率乘积决定，该斜率由以下斜率加权平均：k₁是时间开始时的斜率；k₂是时间中点的斜率，通过采用斜率k₁决定y在的值；k₃也是中点的斜率，但是这次采用斜率k₂决定y的值；k₄是时间终点的斜率，其y值由k₃决定，当四个斜率取平均时，给中点的斜率赋较大的权值。对给定修正的初始值x⁽⁰⁾(1)＝x⁽⁰⁾(1)+ξ代替y₀作为四阶龙格库塔法的初值，利用上式(12)递推预测。

3.如权利要求1所述的基于粒子群优化的灰色预测方法，其特征在于，所述的步骤S2中，利用粒子群优化法对改进的灰色非线性模型中的参数[p_i1,p_i2,…,p_in-1,N_i,ξ_i]进行寻找最优解的方法包含以下步骤：

步骤S2.0、构造适应度函数；

适应度函数选择平均相对误差均值最小，即：

步骤S2.1、初始化粒子群；

种群取值大小M、粒子维数n+1、最大遍历次数i_max、粒子Y轴坐标y_i＝[p_i1,p_i2，…,p_in-1,N_i,ξ_i]和速度大小v_i＝[v_i1,v_i2，…,v_in-1,v_i,n,v_i,n+1]，0.1≤p_i≤1，N∈(0,2)；

步骤S2.2、根据适应度函数计算每个粒子的适应度f_i；

步骤S2.3、根据所有的粒子适应度大小选择各自的极值p_i和全局极值p_g，遍历寻找最优解，优化惯性权重，优化所有粒子的速度和坐标，同时控制其速度符合v∈(-2,2)；

设搜索的全局为D维空间，种群中共有M个粒子，粒子i在第t次遍历的X轴位置为粒子的速度为在每次迭代中粒子i跟踪个体极值点p_i＝[p_i1,p_i2,…,p_iD]^T和全局极值点p_g＝[p_g1,p_g2,…,p_gD]^T来优化自己的位置和速度，原则为：

为了满足粒子的速度在迭代时不会失控，令v∈(-2,2)，公式(13)的为粒子原来的速度；为认知方面，表示粒子原有的信息；为群体信息方面，表示各个粒子之间的彼此分享和信息传递；在式(13)中，t为目前遍历次数；c₁为粒子优化本体的最佳状态的学习因子；c₂为粒子优化所有粒子最佳状态的学习因子；r₁和r₂为[0,1]范围的待定值，ω为惯性权重，随迭代次数线性递减的ω值为：

ω_max＝0.95，ω_min＝0.4，i为目前遍历次数，i_max为最大遍历次数；

步骤S2.4、将各个粒子的适应度f_i同目前的最佳位置f(P_i)作对比，如果f_i＜f(P_i)，则让f_i作为当前粒子的最佳值，f(P_i)＝f_i，同时优化当前粒子的最佳位置P_i＝y_i；

步骤S2.5、把所有粒子的目前的最佳状态和全局最佳状态作对比，如果f(P_i)＜f(P_g)，则将f(P_i)当作目前整体的最佳状态，即f(P_i)＝f(P_g)，同时更新整体的最佳位置V_i＝V_i+c₁×rand()×(pbest_i-X_i)+c₂×rand()×(gbest_i-X_i)，X_i＝X_i+V_i，V_i为粒子的速度，pbest_i，gbest_i为粒子的个体最优值和整体最优值，rand()为介于0～1之间的随机数，X_i为粒子的当前速度，c₁和c₂为学习因子；

步骤S2.6、判断遍历次数是否大于等于最大遍历次数，如果是，则输出全局极值f(P_i)和整体最佳位置；如果否，则进行步骤S2.2，一直到符合遍历要求为止。