一种基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法与流程

本发明属于化工过程软测量建模与应用领域，具体涉及一种基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法。

背景技术：

氨是一种非常重要的基础化工产品，其产量在各类化工产品中位居榜首，在工业上大量用于生产尿素、纯碱、铵态氮肥、硝酸，以及制成纤维、塑料、染料等有机合成工业产品。合成氨的原料包括氮气和氢气，其中氮气可从空气中大量获得，而氢气需要通过专门的制氢设备生产。在多数合成氨工艺中，一段转化炉(简称一段炉)是制备氢气的主要设备，其中的化学反应(催化剂为镍)如下：

上述化学反应是吸热反应，需要对一段炉提供热量。因此，反应温度是保持制氢反应稳定进行的重要因素。一段炉常规的加热方式是在辐射段烧咀燃烧燃料气以及回收的废气。为了维持工艺设定的反应温度，需要将一段炉内的氧气含量控制在指定的范围内。氧气含量(单位为摩尔百分比，mol％)可以通过质量分析仪测定。但质量分析仪不仅价格高昂、测量周期较长，且容易出故障。失去氧气含量的测量值，闭环控制器将无法工作，可能引起一系列的不良后果，如废品率增加、能源消耗增加等引起的环境污染和成本增加等问题，甚至导致安全隐患。

数据驱动的氧气含量软测量模型可以实现氧气含量的在线实时估计，以弥补质量分析仪的不足。其原理是在离线阶段根据氧气含量与过程中容易测量的变量(如温度、压力、流量、液位等参数，又被称之为辅助变量)间的依赖关系建立数学模型，然后在线利用该数学模型对氧气含量进行估计，因此具有无测量滞后、成本低、通用性好、易维护等优点。但由于一段炉燃烧过程十分复杂，且工况切换频繁，生产数据呈现不确定性、多模态、强非线性等特点，使得传统的软测量模型(如主成分分析模型、偏最小二乘模型，神经网络模型、支持向量机模型等)难以取得满意的估计精度。另一方面，由于质量分析仪的测量周期较长，有标签样本(即氧气含量已知的样本)数量很少，导致传统有监督的建模方法由于“过学习”或“欠学习”等原因难以获得准确的模型参数。训练不佳的氧气含量软测量模型必然无法提供满意的估计精度，而且人工整定参数耗时费力，难度非常大。

因此，研究和开发能够同时处理一段炉过程中复杂的不确定性、强非线性、多模态以及有标签样本稀少等问题的氧气含量软测量模型，有助于提高氧气含量的估计精度，从而助力合成氨企业实现安全生产、节能环保、降本增效的目标，是十分必要和迫切的。

技术实现要素：

针对现有技术的不足，本发明提供一种基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法，以贝叶斯高斯混合模型的形式建立氧气含量与辅助变量间的概率化数学模型，根据工况切换自适应分配混合模型的贡献度，有效解决不确定性、非线性、多模态等问题，并通过半监督学习同时利用有标签样本和无标签样本(即氧气含量未知、仅辅助变量已知的样本)，解决有标签样本稀少导致的模型估计精度不高的问题。具体技术方案如下：

一种基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法，其特征在于，包括以下步骤：

(1)选择与一段炉氧气含量y关联的辅助变量其中d表示辅助变量个数；

(2)收集同时包含辅助变量与氧气含量的有标签样本集与仅包含辅助变量的无标签样本集其中nl与nu分别代表有标签样本与无标签样本的数量；

(3)对(xl,yl)与xu做无量纲化处理，将辅助变量样本与氧气含量样本的样本方差转换为单位方差；

(4)给定狄利克雷过程的截断水平m，初始化模型参数的共轭先验分布参数a0、b0、c0、d0、e0、f0、β0、v0、m0、w0与后验分布参数a、b、hk、lk、ck、dk、ek、fk、βk、vk、mk、wk、ωk、ωk，其中，且

α表示狄利克雷过程的集中因子；

χk表示第k个混合模型系数的参数；

μk和λk分别表示第k个混合模型中辅助变量x的分布的均值向量和精度矩阵；

表示第k个混合模型中辅助变量x与氧气含量y之间的线性回归系数；

τk表示的精度矩阵参数；

ηk表示第k个混合模型中测量噪声的精度矩阵系数。

共轭先验分布参数和后验分布参数的含义为：

(a0,b0)和(a,b)分别表示α的先验分布参数和后验分布参数；

(hk,lk)表示χk的后验分布参数；

(m0,β0,w0,v0)和(mk,βk,wk,vk)分别表示(μk,λk)的先验分布参数和后验分布参数；

(c0,d0)和(ck,dk)分别表示τk的先验分布参数和后验分布参数；

(e0,f0)和(ek,fk)分别表示ηk的先验分布参数和后验分布参数；

ωk和ωk表示的后验分布参数；

(5)构建有标签样本(xl,yl)和无标签样本xu及其对应的隐变量的似然函数，其中，zi＝(zi1,…,zim)^t和zj＝(zj1,…,zjm)^t分别表示第i个有标签样本(xi,yi)和第j个无标签样本xj对应的二进制隐变量，且满足

(6)将步骤(3)处理后的训练样本集，步骤(4)中的初始模型参数，以及步骤(5)中构建的似然函数输入半监督贝叶斯高斯混合模型中，通过变分推理学习各个模型参数的最优后验分布q(α)与这里q(·)表示对应变量的最优后验分布。

(7)采集仅包含辅助变量的未知样本，按步骤(3)消除辅助变量的量纲，利用步骤(6)中获得的模型参数的最优后验分布，对氧气含量进行估计。

进一步地，所述的步骤(5)中构建的有标签样本(xl,yl)和无标签样本xu及其对应的隐变量zl、zu的似然函数为：

其中χ＝(χ1,…,χm)，μ＝(μ1,…,μm)，λ＝(λ1,…,λm)，η＝(η1,…,ηm)，表示均值为μk、协方差矩阵为的高斯概率密度函数,

进一步地，所述步骤(6)的模型参数的最优后验分布的参数a，b，hk，lk，ck，dk，ek，fk，βk，vk，mk，wk，ωk和ωk的迭代公式具有如下形式：

a＝a0+m-1

ck＝c0+(d+1)/2

其中，ψ(·)表示digamma函数，i表示相应维数的单位矩阵，1为全1列向量，tr(·)矩阵的迹，表示第k个混合模型的估计误差，这里表示根据的分布计算的期望；κik和κjk的计算方式为

其中

进一步地，所述的步骤(7)具体为：

根据步骤(6)中计算出来的α的后验分布以及狄利克雷过程的性质，各模型混合系数π＝(π1,…,πm)的后验分布可计算为

q(π)＝dir(π|φ1,…,φm)

其中dir(π|φ1,…,φm)代表参数为(φ1,…,φm)的狄利克雷分布，且

然后，根据步骤(6)中计算出的模型参数的后验分布，可得去量纲后的辅助变量xt的边缘分布为

其中表示参数为的学生氏t分布。进而可得xt对应的隐变量zt＝(zt1,…,ztm)的后验分布为

其中zt1,…,ztm均为0-1变量，且满足

然后可求得氧气含量的概率分布，从而得到氧气含量的估计值。

进一步地，所述的氧气含量yt的概率分布为：

其中

因此，可得氧气含量的估计值为

与现有技术相比，本发明的有益效果如下：

1、以混合模型的形式建立氧气含量和辅助变量的数学模型，能够有效解决由工况切换及复杂的燃烧过程引起的多模态、强非线性问题；

2、通过半监督学习可以同时利用有标签样本和无标签样本，解决有标签样本不足引起的模型参数学习不佳的问题，从而提高氧气含量的估计精度；

3、可在一轮训练中同时解决参数学习和模型选择问题，不必遍历所有候选混合模型数量，从而提高训练效率。

附图说明

图1为本发明的基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法的流程图；

图2为某合成氨工厂一段炉装置的工艺原理图；

图3为本发明对氧气含量的估计结果示意图，其中，横坐标代表氧气含量，单位为摩尔百分比(mol％)，纵坐标代表测试样本序号，实线代表氧气含量真实值，虚线代表氧气含量估计值；

图4为高斯混合模型对氧气含量的估计结果示意图，其中，横坐标代表氧气含量，单位为摩尔百分比(mol％)，纵坐标代表测试样本序号，实线代表氧气含量真实值，虚线代表氧气含量估计值；

图5为偏最小二乘模型对氧气含量的估计结果示意图，其中，横坐标代表氧气含量，单位为摩尔百分比(mol％)，纵坐标代表测试样本序号，实线代表氧气含量真实值，虚线代表氧气含量估计值。

具体实施方式

下面结合具体的实施例对本发明的基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法作进一步的阐述。需要指出的是，所描述的实施例仅旨在加强对本发明的理解，并不对本发明起任何限定作用。

一种基于半监督贝叶斯高斯混合模型的合成氨过程一段炉氧气含量在线估计方法，如图1所示，具体包括如下步骤：

(1)在一段炉生产设备中选择与氧气含量y关联的辅助变量其中d表示辅助变量个数；

本实施例根据对某ici-amv工艺(产量为1000t/d)合成氨生产过程一段炉装置(如图2所示)的工艺机理分析，选择对氧气含量影响最大的13个变量作为辅助变量，分别为：到03b001的燃料气流量(x1,位号：fr03001.pv)、到03b001的燃料废气流量(x2,位号：fr03002.pv)、03e005出口处燃料废气的压力(x3,位号：pc03002.pv)、03b001出口处燃料气的压力(x4,位号：pc03007.pv)、03e005出口处燃料废气的温度(x5,位号：ti03001.pv)、03b002e06出口处燃料气的温度(x6,位号：ti03009.pv)、03b001入口处过程气的温度(x7,位号：tr03012.pv)、03b001左上方燃料气的温度(x8,位号：ti03013.pv)、03b001右上方燃料气的温度(x9,位号：ti03014.pv)、03b001正上方混合气的温度(x10,位号：tr03015.pv)、03b001左侧出口转换气的温度(x11,位号：tr03016.pv)、03b001右侧出口转换气的温度(x12,位号：tr03017.pv)、03b001出口转换气的温度(x13,位号：tr03020.pv)。因此辅助变量x＝[x1,…,x13],即

(2)收集同时包含辅助变量与氧气含量的有标签样本集与仅包含辅助变量的无标签样本集其中nl与nu分别代表有标签样本与无标签样本的数量；

本发明从计算机分散控制系统数据库中收集同时包含辅助变量与氧气含量的有标签样本集2000组(记为与仅包含辅助变量的无标签样本集5000组(记为作为训练数据集，其中nl＝2000与nu＝5000分别代表有标签样本与无标签样本的数量。

(3)对(xl,yl)与xu做无量纲化处理，将辅助变量样本与氧气含量样本的样本方差转换为单位方差；

其中去量纲的方法为：

式中，分别代表第l个辅助变量和氧气含量的样本标准差，xn(l)表示第n个样本的第l个辅助变量的采样值。

(4)给定狄利克雷过程的截断水平m，初始化模型参数的共轭先验分布参数与后验分布参数，模型参数的含义为：

α表示狄利克雷过程的集中因子；

χk表示第k个混合模型系数的参数；

μk和λk分别表示第k个混合模型中辅助变量x的分布的均值向量和精度矩阵；

表示第k个混合模型中辅助变量x与氧气含量y之间的线性回归系数；

τk表示的精度矩阵参数；

ηk表第k个混合模型中测量噪声的精度矩阵系数。

在本发明中，各模型参数的共轭先验分布和后验分布确定为：

α的先验分布p(α)和后验分布q(α)均为伽马分布，即p(α)＝gam(α|a0,b0)，q(α)＝gam(α|a,b)，其中gam(α|a0,b0)和gam(α|a,b)分别表示参数为(a0,b0)和(a,b)的伽马分布；

χk的先验分布p(χk)和后验分布q(χk)均为贝塔分布，即p(χk)＝beta(χk|1,α)，q(χk)＝beta(χk|hk,lk)，其中beta(χk|1,α)和beta(χk|hk,lk)分别表示参数为(1,α)和(hk,lk)的贝塔分布；

μk,λk的先验分布p(μk,λk)和后验分布q(μk,λk)均为高斯-威沙特分布，即其中和分别表示参数为(m0,β0,w0,v0)和(mk,βk,wk,vk)的高斯-威沙特分布；

的先验分布和后验分布均为高斯分布，即其中表示均值向量为0、协方差矩阵为的高斯分布，表示均值向量为ωk、协方差矩阵为ωk的高斯分布，i表示相应维数的单位矩阵，；

τk的先验分布p(τk)和后验分布q(τk)均为伽马分布，即p(τk)＝gam(τk|c0,d0)，q(τk)＝gam(τk|ck,dk)，其中gam(τk|c0,d0)和gam(τk|ck,dk)分别表示参数为(c0,d0)和(ck,dk)的伽马分布；

ηk的先验分布p(ηk)和后验分布q(ηk)均为伽马分布，即p(ηk)＝gam(ηk|e0,f0)，q(ηk)＝gam(ηk|ek,fk)，其中gam(ηk|e0,f0)和gam(ηk|ek,fk)分别表示参数为(e0,f0)和(ek,fk)的伽马分布。

因此，在该步骤中，需要初始化先验分布参数，包括以及后验分布参数，包括在本实例中，先验分布的参数设定为a0＝1，b0＝1，c0＝1，d0＝1，e0＝1，f0＝1，β0＝1，v0＝1，m0＝0，w0＝i；后验分布的参数a，b，hk，lk，ck，dk，ek，fk，βk，vk，mk，wk，ωk，ωk的初始值为随机值。

(5)构建有标签样本(xl,yl)和无标签样本xu及其对应的隐变量的似然函数，其中，zi＝(zi1,…,zim)^t和zj＝(zj1,…,zjm)^t分别表示第i个有标签样本(xi,yi)和第j个无标签样本xj对应的二进制隐变量，且满足具有如下形式：

(6)将步骤(3)处理后的训练样本集，步骤(4)中的初始模型参数，以及步骤(5)中构建的似然函数输入半监督贝叶斯高斯混合模型中，通过变分推理学习各个模型参数的最优后验分布q(α)以及具体过程包括变分期望部分和变分最大化部分。

在变分期望部分，需要计算隐变量zl和zu的后验分布q(zl)和q(zu)。根据变分推理的原理可得

其中表示根据的分布计算的期望，χ＝(χ1,…,χm)，μ＝(μ1,…,μm)，λ＝(λ1,…,λm)，η＝(η1,…,ηm)，表示均值为μk、协方差矩阵为的高斯概率密度函数,且

其中ψ(·)表示digamma函数。因此，

其中为简便期间，公式(7)中常数项被省略了；后续计算各参数后验分布时仍省略常数项。

同理，可得zu的后验分布q(zu)如下：

其中

从而可得

其中

在变分最大化部分，需要计算模型参数的后验分布q(θ)。仍然采用变分推理的原理。具体来说，q(α)的求解方法为

因此，α的后验分布q(α)＝gam(α|a,b)的参数更新公式为

lnq(χk)可按照下式计算

因此，χk的后验分布q(χk)＝beta(χk|hk,lk)的参数更新公式为

lnq(μk,λk)可按照下式计算

其中

上式即μk,λk的后验分布的参数更新公式，tr(·)矩阵的迹；

可按照下式计算

其中，1为全1列向量，表示第k个混合模型的估计误差，因此，的后验分布的参数更新公式为

lnq(τk)可按照下式计算

因此，τk的后验分布q(τk)＝gam(τk|ck,dk)的参数更新公式为

lnq(ηk)可按照下式计算

因此，ηk的后验分布q(ηk)＝gam(ηk|ek,fk)的参数更新公式为

通过迭代执行变分期望部分和变分最大化部分，模型参数的后验分布将会收敛。本例中收敛的准则为变分下界的相对增量低于设定阈值(10^-7)。

(7)在线阶段，采集仅包含辅助变量的未知样本xt，按步骤(3)消除辅助变量的量纲，利用步骤(6)中获得的模型参数的最优后验分布，对氧气含量进行估计。具体来说，根据步骤(6)中计算出来的α的后验分布以及狄利克雷过程的性质，各模型混合系数π＝(π1,…,πm)的后验分布可计算为

q(π)＝dir(π|φ1,…,φm)(25)

其中dir(π|φ1,…,φm)代表参数为(φ1,…,φm)的狄利克雷分布，且

然后，根据步骤(6)中计算出的模型参数的后验分布，可得去量纲后的辅助变量xt的边缘分布为

其中表示参数为的学生氏t分布。进而可得xt对应的隐变量zt＝(zt1,…,ztm)的后验分布为

其中zt1,…,ztm均为0-1变量，且满足

最后，可得氧气含量yt的概率分布为

其中

因此，根据式(29)，可得氧气含量的估计值为

为了验证本发明的有效性，从该合成氨工厂一段炉装置计算机分散控制系统中收集额外的有标签样本4000组作为校验样本集，按照步骤(7)，对氧气含量进行估计，平均估计结果如图3所示。同时，图4和图5分别给出了传统的高斯混合模型与偏最小二乘模型对氧气含量的平均估计结果。在高斯混合模型中，混合组分数量通过贝叶斯信息准则设置为12；在偏最小二乘模型中，主成分数量通过交叉验证法设置为10。可以看出，由于不能处理非线性对象，偏最小二乘模型提供的氧气含量的估计值显著偏离真值；而传统的高斯混合模型的估计结果虽然较偏最小二乘模型有所改善，但仍不能令人满意，尤其在第三、四个操作区域(第2500-4000个样本处)。相比之下，本发明提供的方法在所有操作区域提供的氧气含量估计值基本符合其真值。

采用均方根误差(rmse)量化本发明与传统的高斯混合模型、偏最小二乘模型的估计精度，定义如下

其中yt和分别代表第t个校验样本的真实氧气含量与估计值。本发明提供的方法与高斯混合模型、偏最小二乘模型的估计rmse分别为0.6933、1.1515、1.7143。可见，本发明较高斯混合模型、偏最小二乘模型对氧气含量的估计精度有显著提高，估计误差分别降低了约40％和60％。

上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，都落入本发明的保护范围。