CN106600001A - 基于混合高斯分布相关学习机的玻璃窑炉温度预报方法 - Google Patents

Abstract

本发明涉及一种基于混合高斯分布相关学习机的玻璃窑炉温度预报方法，属于自动控制、信息技术和先进制造领域，其特征在于，针对玻璃窑炉温度预报所具有的玻璃窑炉内部反应过程复杂、数据存在复杂非对称噪声及输入变量中包含时间序列变量等建模难题，提出一种基于混合高斯分布下鲁棒相关学习机的玻璃窑炉温度预报方法。该方法采用核函数回归模型作为预报模型结构，使用非零均值混合高斯分布作为预报模型残差项的概率密度分布，将时间序列变量并列列出作为模型的输入变量，并采用贝叶斯推断方法获得模型结构参数的近似后验概率分布，从而获得预报模型的结构参数。本发明可有效应用于玻璃窑炉池底温度的预报，以提高玻璃窑炉控制和操作优化的效果。

Description

基于混合高斯分布相关学习机的玻璃窑炉温度预报方法

技术领域

本发明属于自动控制、信息技术和先进制造领域，具体涉及针对玻璃窑炉池底温度指标预报所具有的玻璃窑炉内部反应过程复杂、数据存在复杂非对称噪声以及输入变量中包含连续时间序列变量等建模难题，提出一种基于混合高斯分布相关学习机的玻璃窑炉温度预报方法。

背景技术

在玻璃窑炉生产过程的控制和优化过程中，玻璃窑炉池底温度的预报起到了关键性的指导作用。但由于实际的工业生产过程涉及复杂的物理化学过程，生产过程数据受环境和测量的影响存在强不确定性和复杂的非对称噪声，以及输入变量中包含连续时间序列变量等原因，采用传统的建模方法会出现预报精度不高、泛化能力差等问题。针对以上问题，设计和使用基于混合高斯分布下鲁棒相关学习机的玻璃窑炉池底温度预报方法是解决上述指标预报难题的有效途径之一。

发明内容

本发明提出一种基于混合高斯分布相关学习机的玻璃窑炉温度预报方法，其特征在于，所述方法是在计算机上依次按以下步骤实现的：

步骤(1)：根据对生产过程的机理分析和先验知识，从历史数据数据库中选择与待预测指标相关的变量作为模型的输入向量x∈R^d，其中d表示输入变量的个数；相对应的待预测指标作为模型的输出变量y，假设一共有N个数据，则指标预报模型的模型训练数据集可以表示为：

其中，x_n和y_n分别表示第n个训练样本的输入向量和输出值；

步骤(2)：由于往往数据中不同特征之间的量纲是不同的，需要先对数据进行归一化处理，即根据如下表达式将各个特征逐一归一化为均值为0，方差为1的形式：

其中x表示某一种输入特征，μ表示该输入特征的均值，δ表示该输入特征的标准差，x^*是该特征归一化以后的值；

步骤(3)：设置混合高斯分布下鲁棒相关学习机相关变量的初始值，具体包括：

a₀，b₀：模型参数先验分布的超参数；

μ₀，β₀：混合高斯分布均值向量先验分布的超参数；

c₀，d₀：混合高斯分布精准度向量先验分布的超参数；

α₀：混合高斯分布混合向量先验分布的超参数；

K：混合高斯分布的组份数；

θ：高斯核函数的参数；

tol：停止迭代的条件参数；

maxiter：最大迭代次数；

在无先验知识下超参数a₀，b₀，β₀，c₀，d₀，α₀可设为10^-4，μ₀可设为0，参数K和θ根据交叉验证方法进行调节以得到最好的预报性能；

步骤(4)：混合高斯分布下鲁棒相关学习机的模型训练过程；

步骤(4.1)：所提出的混合高斯分布下鲁棒相关学习机方法的模型假设为：

其中参数w₁，w₂，...，W_N是基函数φ₁(x)，φ₂(x)，...，φ_N(x)的加权权重，W_N+1是模型的偏置项；基函数φ_n(x)＝k(x_n，x)，函数k(x_n，x)为高斯核函数，其表达式为：

下面用向量w表示[w₁，ω₂，...，w_N+1]组成的向量，向量φ(x)表示[φ₁(x)，φ₂(x)，...，φ_N(x)，1]组成的向量，则模型可以表示为：

f(x；w)＝w^Tφ(x)

步骤(4.2)：初始化模型参数w和误差隶属度参数r_nk；其初始化方法为：模型参数w可使用随机初始化的方法；在初始化误差隶属度参数r_nk时，对第n个训练样本，对r_n1，r_n2，…，r_nK，随机初始化其中的一个为1，其余为0；

步骤(4.3)：根据如下公式分别更新模型内部参数μ_w，∑_w，a_n，b_n，α_k，β_k，m_k，c_k，d_k，r_nk：

其中：

在上述各表达式中，符号<x>表示求变量x的期望；计算上述各式所需的各个期望可通过如下各式进行计算：

<μ_k>＝m_k

其中[S]_nn表示矩阵S对角线上第n个元素；所用到的函数ψ(a)的表达式为：

其中Γ表示伽马函数。

步骤(4.4)：计算下界函数：

其中各个概率分布的定义为：

p(zn：|π)＝Multinomial(zn：|π)

p(π)＝Dir(π|α₀)

q(π)＝Dir(π|α)

q(τ)＝Gam(τ_k|c_k，d_k)

q(θ)＝Gam(θ_n|a_n，b_n)

其中符号Gam，Dir，Multinomial分别表示正态分布，伽马分布，狄里克雷分布和多项式分布；

步骤(4.5)：判断训练过程是否收敛；如果未收敛，则继续迭代步骤(4.3)和步骤(4.4)；如果收敛，则停止迭代，所得模型参数w＝μ_w；判断收敛的准则为：下界函数的变化率小于停止迭代的条件参数tol或迭代次数超过最大迭代次数maxiter。

附图说明

图1：本发明的算法流程图。

图2：本发明的实施流程图。

具体实施方式

为了更好地理解本发明的技术方案，图1中给出了本发明的算法流程图。

以某大型玻璃生产企业为例来说明本发明的实施过程，该实施过程的流程图见图2。首先从该公司的玻璃窑炉生产线的MES系统和数据库中收集相关的数据。然后进行数据预处理，例如输入特征选择，时滞选择，填补缺失数据等。然后按照图1所示的算法流程进行模型的训练以得到预报模型的参数。当需要进行指标预报的时候就将当前的模型输入传递到生产指标预报模块中，通过计算得到最终的模型预测值。由于实际生产过程是实时变化的，需要定期对模型进行更新，即用最新的数据重新对模型进行训练。

为说明本发明所提方法的有效性，从该玻璃生产企业的生产系统中读取一段时间的生产数据，进行数据预处理以后得到400条完整的按时间顺序排列的数据。使用这些训练数据对所提出的鲁棒相关学习机方法进行训练并建立窑炉池底的温度预报模型。

将所提方法与经典相关学习机进行对比，所选择的预测性能指标包括平均绝对误差(Mean Absoulte Error，MAE)，均方根误差(Root Mean Square Error，RMSE)以及R²，其数学定义为：

由定义可以看出，MAE和RMSE越小，说明预报性能越好。而对R²而言，约接近1，模型的预报性能越好。

在此次试验中，把所有400个数据分为3份，首先使用前200个数据训练模型，后201-300的数据作为测试数据调整训练算法的参数。通过交叉验证调整模型参数，所得的模型参数为：经典相关学习机的高斯核函数参数为0.04；所提方法的高斯核函数参数为0.00032，混合高斯分布组分数为3。为了准确的表明模型的泛化性能，使用101-300的数据训练模型，在最后100个未用过的数据作为测试数据，每次测试分别从训练数据和和测试数据中随机抽取50％的数据进行测试，实验重复50次，最终结果为50次测试结果的平均值。表1中列出了实验结果。从结果中可以看出，相比经典的相关学习机算法，本发明所提出的鲁棒相关学习机算法有更好的预报性能。

表1所提方法与经典相关学习机算法在池底温度预报问题上的性能比较

	经典相关学习机	所提出的方法
			MAE	1.2574	0.6313
RMSE	1.4621	0.7563
			R2	0.2040	0.7863

Claims (4)

1.一种基于混合高斯分布相关学习机的玻璃窑炉温度预报方法，其特征在于，包括：

采集并处理与玻璃窑炉池底温度预报相关的生产历史数据；

设置混合高斯分布下鲁棒相关学习机相关变量的初始值；

根据所采集的生产历史数据，对混合高斯分布下鲁棒相关学习机进行训练，得到玻璃窑炉池底温度智能预报模型。

2.根据权利要求1所述的方法，其特征在于，所述采集并处理与玻璃窑炉池底温度预报相关的生产历史数据，具体包括：预报模型输入变量包括天然气用量、窑炉拉引量、单耗、窑炉顶温度和窑炉流液洞温度，设为输入向量x∈R^d，其中，d为输入变量的个数d＝5，预报模型输出变量包括窑炉池底温度，设为y；假设一共有N个数据，则玻璃窑炉池底温度预报模型的训练数据集可以表示为：

${\{x_n,y_n\}}_{n=1}^N$

其中，x_n和y_n分别表示第n个训练样本的输入向量和输出值；之后，根据如下表达式将各个输入变量逐一归一化为均值为0，方差为1的形式：

$x^{*}=\frac{x-\mathrm{\&mu;}}{\mathrm{\&delta;}}$

其中，x表示一种输入变量，μ表示该输入特征的均值，δ表示该输入特征的标准差,x^*是该特征归一化以后的值。

3.根据权利要求1所述的方法，其特征在于，所述设置混合高斯分布下鲁棒相关学习机相关变量的初始值，具体包括：

a₀，b₀：模型参数先验分布的超参数；

μ₀，β₀：混合高斯分布均值向量先验分布的超参数；

c₀，d₀：混合高斯分布精准度向量先验分布的超参数；

α₀：混合高斯分布混合向量先验分布的超参数；

K：混合高斯分布的组份数；

θ：高斯核函数的参数；

tol：停止迭代的条件参数；

maxiter：最大迭代次数；

在无先验知识下超参数a₀，b₀，β₀，c₀，d₀，α₀可设为10^-4，μ₀可设为0，参数K和θ根据交叉验证方法进行调节以得到优化后的预报模型。

4.根据权利要求1所述的方法，其特征在于，所述对混合高斯分布下鲁棒相关学习机进行训练，是按如下步骤进行的：

步骤1：所提出的混合高斯分布下鲁棒相关学习机方法的模型假设为：

$f(x;w_1,w_2,\mathrm{...},w_{N+1})=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_n{\mathrm{\&phi;}}_n\left(x\right)+w_{N+1}$

其中，参数w₁，w₂，...，w_N是基函数φ₁(x)，φ₂(x)，...，φ_N(x)的加权权重，w_N+1是模型的偏置项；基函数φ_n(x)＝k(x_n，x)，函数k(x_n，x)为高斯核函数，其表达式为：

${\mathrm{\&phi;}}_n\left(x\right)=\exp \&lsqb;-\frac{{(x_n-x)}^T(x_n-x)}{\mathrm{\&theta;}}\&rsqb;$

下面用向量w表示[w₁，w₂，...，w_N+1]组成的向量，向量φ(x)表示[φ₁(x)，φ₂(x)，...，φ_N(x)，1]组成的向量，则模型可以表示为：

f(x；w)＝w^Tφ(x)；

步骤2：初始化模型参数ω和误差隶属度参数r_nk；其初始化方法为：模型参数w可使用随机初始化的方法；在初始化误差隶属度参数r_nk时，对第n个训练样本，对r_n1，r_n2，…,r_nK,随机初始化其中的一个为1，其余为0；

步骤3：根据如下公式分别更新模型内部参数μ_w，∑_w，a_n，b_n，α_k，β_k，m_k，c_k，d_k，r_nk：

${\mathrm{\&mu;}}_w={\mathrm{\&Sigma;}}_w\underset{k=1}{\overset{K}{\&Sigma;}}\{<{\mathrm{\&tau;}}_k>\underset{n=1}{\overset{N}{\&Sigma;}}\&lsqb;<z_{nk}>\mathrm{\&phi;}\left(x_n\right)<y_n-{\mathrm{\&mu;}}_k>\&rsqb;\}$

${\mathrm{\&Sigma;}}_w={\&lsqb;<\mathrm{\&Theta;}>+\underset{n=1}{\overset{N}{\&Sigma;}}\underset{k=1}{\overset{K}{\&Sigma;}}(<z_{nk}><{\mathrm{\&tau;}}_k>\mathrm{\&phi;}(x_n\left)\mathrm{\&phi;}{\left(x_n\right)}^T\right)\&rsqb;}^{-1}$

$a_n=\frac{1}{2}+a_0$

$b_n=b_0+\frac{1}{2}<{w_n}^2>$

$\begin{array}{cc}\mathrm{\&alpha;}={({\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,...,{\mathrm{\&alpha;}}_K)}^T& {\mathrm{\&alpha;}}_k={\mathrm{\&alpha;}}_{0k}+\underset{n=1}{\overset{N}{\&Sigma;}}<z_{nk}>\end{array}$

${\mathrm{\&beta;}}_k=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}<z_{nk}>+{\mathrm{\&beta;}}_0$

$m_k=\frac{1}{{\mathrm{\&beta;}}_k}\&lsqb;\underset{n=1}{\overset{N}{\&Sigma;}}(<z_{nk}><y_n-w^T\mathrm{\&phi;}(x_n)>)+{\mathrm{\&beta;}}_0{\mathrm{\&mu;}}_0\&rsqb;$

$c_k=c_0+\frac{1}{2}+\frac{1}{2}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}<z_{nk}>$

$d_k=d_0+\frac{1}{2}\{\underset{n=1}{\overset{N}{\&Sigma;}}\&lsqb;<z_{nk}><{(y_n-w^T\mathrm{\&phi;}(x_n\left)\right)}^2>\&rsqb;-{\mathrm{\&beta;}}_k{m_k}^2+{\mathrm{\&beta;}}_0{{\mathrm{\&mu;}}_0}^2\}$

$r_{nk}=\frac{{\mathrm{\&rho;}}_{nk}}{{\mathrm{\&Sigma;}}_{k=1}^K{\mathrm{\&rho;}}_{nk}}$

其中：

${\mathrm{ln\&rho;}}_{nk}=\frac{1}{2}<{\mathrm{ln\&tau;}}_k>-\frac{1}{2}\ln 2\mathrm{\&pi;}-\frac{1}{2}<{\mathrm{\&tau;}}_k><{(y_n-w^T\mathrm{\&phi;}(x_n)-{\mathrm{\&mu;}}_k)}^2>+<{\mathrm{ln\&pi;}}_k>$

在上述各表达式中，符号<x>表示求变量x的期望；计算上述各式所需的各个期望可通过如下各式进行计算：

$<(y_n-w^T\mathrm{\&phi;}{\left(x_n\right)}^2>={y_n}^2-2y_n{\mathrm{\&mu;}}_w^T\mathrm{\&phi;}(x_n)+\mathrm{\&phi;}{\left(x_n\right)}^T\&lsqb;{\mathrm{\&Sigma;}}_w+{\mathrm{\&mu;}}_w{\mathrm{\&mu;}}_w^T\&rsqb;\mathrm{\&phi;}\left(x_n\right)$

$<{w_n}^2>={\&lsqb;{\mathrm{\&Sigma;}}_w+{\mathrm{\&mu;}}_w{\mathrm{\&mu;}}_w^T\&rsqb;}_{nn}$

<μ_k>＝m_k

$\begin{array}{c}<{(y_n-w^T\mathrm{\&phi;}(x_n)-{\mathrm{\&mu;}}_k)}^2>={y_n}^2-2y_n{\mathrm{\&mu;}}_w^T\mathrm{\&phi;}\left(x_n\right)+\mathrm{\&phi;}{\left(x_n\right)}^T\&lsqb;{\mathrm{\&Sigma;}}_w+{\mathrm{\&mu;}}_w{\mathrm{\&mu;}}_w^T\&rsqb;\mathrm{\&phi;}\left(x_n\right)\\ -2m_k\left(y_n-{\mathrm{\&mu;}}_w^T\mathrm{\&phi;}\left(x_n\right)\right)+{\left({\mathrm{\&beta;}}_k{\mathrm{\&tau;}}_k\right)}^{-1}+{m_k}^2\end{array}$

$<{\mathrm{\&theta;}}_n>=\frac{a_n}{b_n}$

其中，[S]_nn表示矩阵S对角线上第n个元素；所用到的函数ψ(a)的表达式为：

$\mathrm{\&psi;}\left(a\right)=\frac{d}{da}\ln \mathrm{\&Gamma;}\left(a\right)$

其中，表示伽马函数；

步骤4：计算下界函数：

其中，各个概率分布的定义为：

p(z_n：|π)＝Multinomial(z_n：|π)

p(π)＝Dir(π|α₀)

$p\left(\mathrm{\&tau;}\right)=\underset{k=1}{\overset{K}{\&Pi;}}Gamma\left({\mathrm{\&tau;}}_k\right|c_0,d_0)$

$p\left(\mathrm{\&theta;}\right)=\underset{n=1}{\overset{N+1}{\&Pi;}}Gamma\left(\mathrm{\&theta;}\right|a_0,b_0)$

$\ln q\left(Z\right)=\underset{n=1}{\overset{N}{\&Sigma;}}\underset{k=1}{\overset{K}{\&Sigma;}}{r}_{nk}\ln r_{nk}$

q(π)＝Dir(π|α)

q(τ)＝Gam(τ_k|c_k，d_k)

q(θ)＝Gam(θ_n|a_n，b_n)

其中，符号Gam，Dir，Multinomial分别表示正态分布，伽马分布，狄里克雷分布和多项式分布；

步骤5：判断训练过程是否收敛；如果未收敛，则继续迭代步骤3和步骤4；如果收敛，则停止迭代，所得模型参数ω＝μ_ω；判断收敛的准则为：下界函数的变化率小于停止迭代的条件参数tol或迭代次数超过最大迭代次数maxiter。