CN107704962B - 一种基于不完整训练数据集的蒸汽流量区间预测方法 - Google Patents

Abstract

一种基于不完整训练数据集的蒸汽流量区间预测方法。针对含缺失点的蒸汽流量时间序列数据，利用相空间重构技术构造输入输出训练样本集，并建立一个初始相关向量机回归模型。然后通过贝叶斯方法推导出缺失输出的后验概率分布；并由缺失输出的后验概率的均值来填补与输出相对应的缺失输入点，不断迭代上述推导与填补过程来训练预测模型。通过单步迭代预测方法来预测未来一段趋势的预测值，其预测区间边界由预测分布的均值和方差决定。这种预测方法可以很好地应用在含有缺失点时间序列数据集中，其预测结果更好地辅助调度人员进行冶金蒸汽能源系统的平衡优化工作。

Description

一种基于不完整训练数据集的蒸汽流量区间预测方法

技术领域

本发明属于信息技术领域，涉及到相关向量机回归，贝叶斯后验概率推理等技术，是一种基于不完整数据集的冶金企业蒸汽流量区间预测方法。针对含缺失点的蒸汽流量时间序列数据，利用相空间重构技术构造输入输出训练样本集，并建立一个初始相关向量机回归模型。然后通过贝叶斯方法推导出缺失输出的后验概率分布；并由缺失输出的后验概率的均值来填补与输出相对应的缺失输入点，不断迭代上述推导与填补过程来训练预测模型。通过单步迭代预测方法来预测未来一段趋势的预测值，其预测区间边界由预测分布的均值和方差决定。这种预测方法可以很好地应用在含有缺失点时间序列数据集中，其预测结果更好地辅助调度人员进行冶金蒸汽能源系统的平衡优化工作，在冶金工业其他能源介质数据中亦可推广应用。

背景技术

钢铁工业是高能耗产业，一直以来都面临着节能降耗的问题。而蒸汽作为钢铁生产中重要的二次能源，不仅被直接应用到各个环节的生产过程中，还被用来转化为电能以保证钢厂的电能供应。因此，针对蒸汽系统进行合理调度能够最大限度的利用蒸汽资源，而蒸汽流量数据的准确预测对有效调度方案的实施有着重要作用。当预测到蒸汽供应大于蒸汽消耗时，可以及时调整消耗源(如电厂，锅炉等)使其增加对蒸汽的消耗；反之，则减少蒸汽消耗，以达到蒸汽产消平衡，提高蒸汽利用率的目的(Y.Liu,Q.Liu,J.Zhao,et al,(2012).Data-driven based model for flow prediction of steam system in steelindustry[J].Information Sciences,193(11):104-114.)。然而，钢铁工业生产过程环境复杂恶劣，通过工业传感器采集的蒸汽流量数据由于传感器故障等原因常常含有缺失值，进而造成时间序列训练数据集的不完整，采用不完整数据集很难直接构造预测模型进行蒸汽流量预测。

常见的处理含缺失值的方法是数据填补，此方案的思路是先将缺失点进行数据补全，再使用填补后的数据进行预测建模。典型的数据填补方法有单一填补，例如，最临近点填补，样条插值等，以及多重填补方法等。(J.L.Schafer,J.W.Graham.(2002)Missingdata:our view of the state of the art.Psychol Methods,7(2):147.)近年来，有改进的最近邻方法，贝叶斯网方法以及矩阵补全方法等被应用于数据填补。另外，目前针对冶金数据的预测模型包括支持向量机，神经网络等基于数据的机器学习方法。

但是，以上数据填补的方法有如下缺点：数据填补方案会导致原数据集的方差发生变化，从而影响建模的精度；而且不同的数据填补方法的填补效果不同，也会对模型的建立产生不同的影响，且针对不同的数据需要分别确定最优的填补方案，过程繁琐。因此，如果能够将含缺失点的数据集直接用来建模，在建模的过程中去填补缺失点就可以避免提前数据填补的步骤简化方法，并且可以使填补方案适用于提高模型预测精度。另一方面，现有预测模型(如，支持向量机，神经网络等)大多数只能给出预测结果的单一点估计值，并不能给出预测可靠性的度量，而工业应用现场往往需要考虑预测的不确定性的。(J.Zhao,Q.Liu,W.Wang,et al.(2012)Hybrid neural prediction and optimized adjustmentfor coke oven gas system in steel industry[J].IEEE Transactions on NeuralNetworks and Learning Systems,23(3):439.)

发明内容

本发明主要解决基于不完整数据集的冶金企业蒸汽流量区间预测问题。针对从工业现场采集的含缺失点冶金企业蒸汽流量时间序列数据，首先将用0填补后的蒸汽流量数据通过相空间重构技术构造输入输出训练样本集，建立一个初始相关向量机回归模型。并通过贝叶斯方法推导出缺失输出的后验概率分布(多维高斯分布)；随后由缺失输出的后验概率分布均值来填补与输出相对应的缺失输入点。迭代上述过程来训练并得出最终的预测模型。

本发明的技术方案：

一种基于不完整训练数据集的蒸汽流量区间预测方法，步骤如下：

1.用0来填补不完整的蒸汽流量训练数据集，再选定一个合适的嵌入维数通过相空间重构方法构造输入输出样本集。这样缺失点分别存在于输入集和输出集中，并且各个缺失点在输入输出中存在相互对应关系，记录下这种一一对应关系。

2.构造初始预测模型。首先由零填补后的输入输出样本集构造一个初始相关向量机回归模型。

3.根据贝叶斯方法，推导出缺失输出的后验概率与模型权值为联合多维高斯分布形式。并由相应的缺失输出后验概率均值来填补训练集中的缺失输入点。

4.计算相关向量机模型的边缘似然度函数值。

5.反复迭代第3和第4步骤，直至边缘似然度函数值收敛为止。输出最终预测模型。

6.预测阶段：当系统采集到新的蒸汽流量数据时，将此数据输入到预测模型中并预测出下一时刻的流量值，以此单补迭代预测未来一段时间的蒸汽流量值。同时由预测值的概率分布构造预测区间的上下限从而得出预测区间。

本发明的有益效果：本发明方法针对含缺失点的蒸汽流量数据采取直接建模的手段。这样即省去了提前数据填补的步骤，减少了挑选最优填补方法的时间；同时由于建模过程的目标是最小化预测误差，并且缺失数据点的处理是与预测建模同时进行，因此本发明的缺失填补过程能够最大程度的提高预测精度。另外，本发明不同于现有大多数预测方法只给出预测值的点估计的方式，它能给出预测精度的可靠性度量(预测区间)，从而给现场调度人员更多的分析判别依据。

附图说明

图1(a)是表示含缺失点的蒸汽流量时间序列示意图。

图1(b)是由相空间重构构造的输入输出训练样本集示意图。

图2为本发明应用流程图。

图3(a)为训练集数据在缺失10％时的贝叶斯神经网络方法效果图。

图3(b)为训练集数据在缺失10％时的Bootstrap方法效果图。

图3(c)为本发明的训练集数据在缺失10％时的区间预测效果图。

具体实施方式

为了更好地理解本发明的技术方案与具体实施方式，以下结合附图对本发明的实施方式作进一步描述。由工业现场采集的数据通常都含有高噪声，并且常有缺失值的情况，这样的不完整数据集无法直接用来构造预测模型。如附图1(a)展示了一条含缺失点的时间序列，其中白色方框表示数据缺失点，而灰色填充方框表示观测到的数据点。在现有方法中不得不先通过数据填补使数据集完整，再应用填补后的数据集参与训练预测模型。但这种方案比较繁琐，同时填补方法的选择也对最终预测结果产生很大影响。另一方面，针对工业数据含有高噪声的情况，现场操作人员除了需获得变量的实时预测值之外，对预测精度的可靠性也越来越重视。针对以上情况，本发明构造了一种基于不完整数据集的直接训练相关向量机预测模型的方法，并在预测阶段能够给出预测区间。

本发明的具体实施步骤如下：

步骤1：数据准备与初始填补

从工业现场实时关系数据库提取一段时间的冶金蒸汽流量时间序列数据。检测并记录时间序列中缺失点的位置，并用零值填补所有缺失点。

步骤2：基于相空间重构的输入输出样本集构建

相空间重构是时间序列预测中常用的训练样本集构造方法。假设有一条含缺失点的离散时间序列{u(1),u(2),...,u(k),...,u(N_u)}，其中u(k)表示第k个采样时刻的时间序列值，N_u是序列的总长度；由式(1)构造第i个相矢量s_i(p)：

s_i(p)＝[u(i),u(i+τ),u(i+2τ),...,u(i+(p-1)τ)]^T (1)

其中i＝1,2,...,N，而N＝N_u-(p-1)τ/Δt，τ表示延迟时间，Δt表示数据的采样间隔，p是相矢量的嵌入维数，也是矢量s_i(p)的长度；第i个相矢量s_i(p)同时表示输入样本集中的第i个输入样本，即x_i＝s_i(p)，有x_i＝[x_i1,x_i2,...,x_ic,...,x_ip]^T，其中x_ic是向量x_i中的第c个分量，c＝1,2,…,p。向量x_i对应的目标值t_i为式(2)所示：

t_i＝u(i+p×τ) (2)

此外，N个输入样本的矩阵形式为X＝[x₁,x₂,...,x_i,...,x_N]^T；N个目标值所组成的向量为t＝[t₁,t₂,...,t_i,...,t_N]^T。附图1(b)表示了由相空间重构技术构造的输入输出样本集，延迟时间τ取为1，此样本集对应于图1(a)中所示的含缺失点的时间序列。可以看出，输入和输出中的缺失点相互耦合，如图1(b)中的缺失点“miss k”不管在输入还是输出中它们的值都是相等的。需要事先将这种对应关系记录下来。

步骤3：构建初始相关向量机回归模型

构建初始相关向量机模型即构建输入输出变量之间的函数关系。相关向量机回归模型描述如下。假设第i个目标值t_i与第i个输入向量x_i有如下的映射关系。

其中，j＝0,1,2,…,N，x_j是第j个输入向量，且规定K(x_i,x₀)＝1。另外，第i个基函数向量φ_i表示为φ_i＝[1,K(x_i,x₁),K(x_i,x₂),...,K(x_i,x_j),...,K(x_i,x_N)]^T，w＝[w₀,w₁,...,w_j,...,w_N]^T表示权值向量，K(x_i,x_j)＝exp{-||x_i-x_j||²/2l²}是高斯核函数，其中l≠0是核函数参数；ε_i表示第i个均值为0，方差为σ²(β≡σ^-2)的高斯分布噪声项。将权值向量w视为连续随机变量，假设w中的每一个元素分别服从是均值为0而方差不同的正态分布

其中，参数α＝[α₀,α₁,...,α_j,...,α_N]^T，α_j是第j个权值w_j的先验方差的倒数，且有对角阵A＝diag(α)。并且同时将参数α和β视为随机变量，并假设其服从伽马分布。

步骤4：基于不完整数据集的相关向量机回归模型构建

本步骤由步骤2所得输入输出样本集来构建基于不完整数据集的相关向量机回归模型。假设训练输出集由缺失点和可观测点组成，有

其中t_[m]和t_[obs]分别表示输出集中缺失点和可观测点的向量，输出缺失点和观测点的个数分别为N_[m]和N_[yobs]，有N＝N_[m]+N_[yobs]。这两个向量分别对应的输入样本集矩阵是

和

而且有

此外，所有缺失的和观测到的输入点组成的集合分别表示为向量x_[m]和x_[obs](需要注意的是，矩阵

和

与向量x_[m]和x_[obs]不同，它们中的元素既包含缺失点也有可观测点)。

本发明将噪声方差向量写为β＝[β_[m],β_[m],...,β_[m],β_[obs],...,β_[obs],β_[obs]]^T，且有β_[obs]≡β≡σ^-2。另外，β还可以写成

向量β_[m]和β_[obs]的长度分别为N_[m]和N_[yobs]。写成对角矩阵的形式是B_[m]＝diag(β_[m])，B_[obs]＝diag(β_[obs])和B＝diag(β)。由于t_[m]是缺失的，并没有受到测量噪声的干扰，因此本发明将β_[m]中的元素β_[m]限定为大于10⁸的值以此来保证对应的缺失输出点的噪声方差很小。似然函数表达为如下形式。

其中，核函数矩阵Φ＝[φ₁,φ₂,...,φ_i,...,φ_N]^T。将β_[obs]的先验概率分布设置为伽马分布。由于

和w的先验分布是两个高斯分布的卷积的形式，本发明给出以下联合概率分布。

p(w,t_[m],t_[obs]|α,β_yobs)＝N(w,t_[m],t_[obs]|μ_[w,m,obs],Σ_[w,m,obs]) (6)

其中，μ_[w,m,obs]＝[0 0 0]^T，

而且有

由式(6)可以看出t_[m]，t_[obs]和w是联合高斯分布，推导出w和t_[m]的联合后验概率分布也是一个多维高斯分布的形式，如下所示。

p(w,t_[m]|t_[obs],α,β_[obs])＝N(w,t_[m]|μ_[w,m|obs],Σ_[w,m|obs]) (7)

其中，

针对缺失输入点x_[m]的优化问题，考虑到相空间重构方法构造的样本集的特点，每一个输入缺失点都有在输出缺失集中相对应的等值缺失点(输入和输出缺失对应点的值相等)，因此，本发明将x_[m]中每一个元素的后验概率都用由(7-9)式计算出的相应的缺失输出点的后验概率来代替，并将x_[m]的后验均值填补到输入缺失点的相应位置，并进一步更新核矩阵Φ。

其中，

x_[m]集合中共有

个输入缺失点，x_[m],v表示x_[m]中的第v个元素，

表示x_[m],v的后验概率均值；而t_{[m],corr_v}表示第v个缺失输入点在缺失输出集t_[m]中的对应变量；t_{[m],corr_v}的后验概率p(t_{[m],corr_v}|t_[obs],α,β_[obs])由式(8)和(9)得出。

将w和t_[m]积分消去，本发明给出对数边缘似然度函数。

其中，协方差矩阵

矩阵大小为N_[yobs]×N_[yobs]，I是单位矩阵，log(·)表示自然对数运算。通过求偏导最大化

来找到α和β_[obs]的极值点α^*和

从而得出它们的迭代式如下所示。

其中

且

以下简要整理建模步骤：

①初始化：向量α中的元素初始化为值1/N²，β_[obs]初始化为10倍的t_[obs]的方差，相关向量集初始化为所有输入样本的集合(此时相关向量集中有N个样本)；

②使用(8)和(9)式更新w和t_[m]的联合后验概率的均值和协方差矩阵μ_[w,m|obs]和Σ_[w,m|obs]。

③使用(12)和(13)式更新α和β_[obs]。

④如果α_j＞M(M是一个大于10¹⁰的值)，将第j个样本从相关向量集中除去。

⑤通过(10)式来得到x_[m]中每一个元素的后验概率均值，并填补到输入缺失点的相应位置，并进一步更新核矩阵Φ。

⑥由(11)式计算对数边缘似然度函数的值。

⑦如果α中有元素相对上次迭代的变化或边缘似然度函数值的变化小于给定的阈值，终止迭代过程，得出最终预测模型，相关向量集中剩余的向量为最终的相关向量；否则，转到第②步。

步骤5：单步迭代预测得出未来一段时间的预测区间

当给定一个新的输入x_*(此输入向量中没有缺失点)，β_[obs]的值使用可

来代替。得出对应的输出t_*的预测分布近似为一个正态分布

其中，

φ_*＝[1,K(x_*,x₁),K(x_*,x₂),...,K(x_*,x_N)]^T。μ_[w|obs]是权值向量w的后验概率的均值。此时，预测上下限分别表示为

其中，U_limit和L_limit分别表示置信区间为(1-χ％)的预测区间上下限，同时

是标准高斯分布中对应概率为

的分位数。

未来一段时间的预测均值由单步迭代预测的方式给出。在如此的机制下进行预测，每次预测输出的均值是一个点，且下一轮的预测会将上一轮的预测均值结果作为输入，进行循环滚动预测。

附图3(a)、(b)、(c)分别表示训练集数据在缺失10％时的基于Bayesian的神经网络方法，Bootstrap方法和本发明的区间预测效果比较图。对比方法包括常用的基于Bayesian的回声状态网络(ESN,Echo State Network)区间预测，以及基于Bootstrap的ESN区间预测方法。预测结果比较可参见表1，其中RMSE为均方根误差，用来表示预测均值的误差情况，CWC是一种基于预测区间覆盖概率和平均区间宽度的综合指标。

其中Y_b为第b个真实观测值，F_b为第b个模型输出值，N_test为被估计的样本点个数。

CWC＝NMPIW(1+γ×exp[-η(PICP-μ)]) (18)

其中，

而η和μ是两个大于0的超参数。其中，当第b个目标值在预测范围内时c_b等于1，反之，c_b等于0；而U_b和L_b分别是第b个点的预测上限和下限。PICP表示被预测区间覆盖的真实目标值的个数，NMPIW是归一化的预测区间平均宽度，其中，G是目标值的取值范围。CWC指标值越小区间预测效果越好。并且当PICP值越大且NMPIW值越小时，区间预测的效果越好。

表1三种方法区间预测结果比较(NNI表示最近邻点填补方法；SI表示样条插值填补方法)

Claims (1)

1.一种基于不完整训练数据集的蒸汽流量区间预测方法，其特征在于如下步骤：

(1)数据准备与初始填补

从工业现场实时关系数据库提取一段时间的冶金企业蒸汽流量时间序列数据，检测并记录时间序列中缺失点的位置，并用零值填补所有缺失点；

(2)基于相空间重构的输入输出样本集构建

假设有一条含缺失点的离散时间序列{u(1),u(2),...,u(k),...,u(N_u)}，其中u(k)表示第k个采样时刻的时间序列值，N_u是序列的总长度；由式(1)构造第i个相矢量s_i(p)：

s_i(p)＝[u(i),u(i+τ),u(i+2τ),...,u(i+(p-1)τ)]^T (1)

其中，i＝1,2,...,N，而N＝N_u-(p-1)τ/Δt，τ表示延迟时间，Δt表示数据的采样间隔，p是相矢量的嵌入维数，也是矢量s_i(p)的长度；第i个相矢量s_i(p)同时表示输入样本集中的第i个输入样本，即x_i＝s_i(p)，有x_i＝[x_i1,x_i2,...,x_ic,...,x_ip]^T，其中x_ic是向量x_i中的第c个分量，c＝1,2,…,p；向量x_i对应的目标值t_i为式(2)所示：

t_i＝u(i+p×τ) (2)

此外，N个输入样本的矩阵形式为X＝[x₁,x₂,...,x_i,...,x_N]^T；N个目标值所组成的向量为t＝[t₁,t₂,...,t_i,...,t_N]^T；对于由含缺失点的时间序列构造的输入输出样本集，事先记录下输入输出中缺失点的对应关系；

(3)构建初始相关向量机回归模型

假设第i个目标值t_i与第i个输入向量x_i有如下的映射关系：

其中，j＝0,1,2,…,N，x_j是第j个输入向量，且规定K(x_i,x₀)＝1；另外，第i个基函数向量φ_i表示为φ_i＝[1,K(x_i,x₁),K(x_i,x₂),...,K(x_i,x_j),...,K(x_i,x_N)]^T，w＝[w₀,w₁,...,w_j,...,w_N]^T表示权值向量，K(x_i,x_j)＝exp{-||x_i-x_j||²/2l²}是高斯核函数，其中，l≠0是核函数参数；ε_i表示第i个均值为0，方差为σ²的高斯分布噪声项，β≡σ^-2；对w中的每一个元素分别附加均值为0而方差不同的正态先验分布：

其中，参数α＝[α₀,α₁,...,α_j,...,α_N]^T，α_j是第j个权值w_j的先验方差的倒数，且有对角阵A＝diag(α)；并且同时将参数α和β视为随机变量，并假设它们服从伽马分布；

(4)基于不完整数据集的相关向量机回归模型构建

构建最终预测模型就是要通过输入输出样本集来迭代计算如下变量：

权值向量w；

缺失输出向量t_[m]，即输出集中所有缺失点组成的向量；

缺失输入向量x_[m]，即输入集中所有缺失点组成的向量；

噪声方差参数β的点估计值；

具体步骤如下：

①初始化：向量α中的元素初始化为值1/N²，β_[obs]初始化为10倍的t_[obs]的方差，β_[obs]是输出集中可观测点的噪声方差的倒数，有β_[obs]≡σ^-2；其中，t_[obs]表示输出集中所有可观测值组成的向量；相关向量集初始化为所有输入样本的集合，此时相关向量集中有N个样本；

②分别使用式(5)和(6)更新w和t_[m]的联合后验概率的均值μ_[w,m|obs]和协方差矩阵Σ_[w,m|obs]；

其中，β＝[β_[m],β_[m],...,β_[m],β_[obs],...,β_[obs],β_[obs]]^T，β还可写成

写成对角矩阵的形式是B_[m]＝diag(β_[m])，B_[obs]＝diag(β_[obs])和B＝diag(β)；矩阵B_[m]和B_[obs]大小分别为N_[yobs]×N_[yobs]和N_[m]×N_[m]，而N_[yobs]和N_[m]分别表示输出集中可观测点和缺失点的个数，有N＝N_[yobs]+N_[m]，将β_[m]设置为不小于10⁸的值；另外，Φ表示核函数矩阵，有Φ＝[φ₁,φ₂,...,φ_i,...,φ_N]^T，并且将Φ表示成分块矩阵的形式为

其中

和

是由缺失输出点和输出可观测点分别在输出样本集中的位置索引所提取的核矩阵Φ中以行为单位的行向量所组成的各自的新矩阵，大小分别为N_[m]×(N+1)和N_[yobs]×(N+1)；

③采用(7)和(8)式更新α_j和β_[obs]；

其中，

且

Tr(·)表示矩阵的求迹运算；

④如果α_j＞M，M设置为大于10¹⁰的值，将第j个样本从相关向量集中除去；

⑤通过(9)式来计算x_[m]中每一个元素的后验概率均值，并填补到输入集中输入缺失点的相应位置，并进一步更新核矩阵Φ；

其中，

x_[m]集合中共有

个输入缺失点，x_[m],v表示x_[m]中的第v个元素，

表示x_[m],v的后验概率均值；而t_{[m],corr_v}表示第v个缺失输入点在缺失输出集t_[m]中的对应变量值；t_{[m],corr_v}的后验概率p(t_{[m],corr_v}|t_[obs],α,β_[obs])由式(5)和(6)得出；

⑥由(10)式计算对数边缘似然度函数的值；

其中，协方差矩阵

矩阵大小为N_[yobs]×N_[yobs]，I是单位矩阵，log(·)表示自然对数运算；

⑦如果α中有元素相对上次迭代的变化或边缘似然度函数值的变化小于给定的阈值，终止迭代过程，得出最终预测模型；否则，转到第②步；

(5)单步迭代预测得出未来一段时间的预测区间

给定一个新的输入样本x_*，此输入向量中没有缺失点，得出对应的输出t_*的预测分布近似为一个正态分布

其中，

φ_*＝[1,K(x_*,x₁),K(x_*,x₂),...,K(x_*,x_N)]^T，μ_[w|obs]是权值向量w的后验概率的均值，

由式(8)计算得到；此时，预测上下限分别表示为

其中，U_limit和L_limit分别表示置信区间为(1-χ％)的预测区间上下限，同时

是标准高斯分布中对应概率为

的分位数；未来一段时间的预测均值由单步迭代的方式给出；每一次预测输出的预测均值是一个点，且下一轮的预测会将上一轮的预测结果作为输入，进行循环滚动预测。

Images