CN112562797B - 沉铁过程出口离子预测方法及系统 - Google Patents

Abstract

本发明涉及湿法炼锌沉铁过程控制领域，公开一种沉铁过程出口离子预测方法及系统，以解决现有的沉铁过程出口离子浓度预测模型精确度不高的技术问题。本发明方法包括：围绕c个中心样本选取相似度高的L个时间连续样本构成c个子训练样本集；将对应不同工况的各子训练样本集分别进行训练得到c个动态概率隐变量子模型；获取用于预测离线沉铁过程出口离子浓度特征的无标签的测试样本数据，对同一测试样本数据进行局部邻域标准化处理后输入到各动态概率隐变量子模型得到出口离子浓度预测值的概率分布；采用贝叶斯决策进行联合估计将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，得到对应所述测试样本数据最终的出口离子浓度的预测值。

Description

沉铁过程出口离子预测方法及系统

技术领域

本发明涉及湿法炼锌沉铁过程控制领域，尤其涉及一种沉铁过程出口离子预测方法及系统。

背景技术

锌是一种能应用在各个领域的重要金属。湿法冶炼由于高产出和低能耗成为最主要的锌冶炼方法，主要包括磨矿、浸出、净化、电解等工序。该方法中的硫酸锌溶液需要经过净化除杂后，方可进行电解得到锌单质。目前，针铁矿法常用来去除硫酸锌溶液中的主要杂质铁离子。沉铁过程的主要生产设备为连续搅拌釜式反应器(CSTR)，然而单个反应器不能将浸出液中过量铁离子直接降低到工艺要求范围内，因此需要在四个由高到低级联的反应器中进行除铁。前一反应器的出口离子浓度为下一个反应器的入口离子浓度，各个反应器出口铁离子浓度需按要求依次降低，保证硫酸锌溶液离开最后一个反应器后其中铁离子含量降低至工艺指标要求范围内。为了达到这个目的，需根据每个反应器入口离子浓度(即上一反应器的出口离子浓度)分别对其中添加的氧气和焙砂进行调节，然而在实际沉铁过程中，反应器自身的密闭性及检测装置的局限性使得反应器出口溶液中各离子的浓度只能通过人工定期抽样检测获得。导致对操作参数(氧气和焙砂)的调节具有很大的滞后性，降低沉铁效率。因此，建立沉铁过程模型，预测反应器出口离子浓度，对调整操作参数的盲目性、提高沉铁效率具有重大意义。

实际沉铁过程中，合理控制Fe²⁺的氧化速率即氧化反应速率是实现高效除铁的关键。Fe²⁺的氧化速率过快，会使Fe³⁺的含量过高，生成氢氧化铁胶体，影响水解反应的正常进行；氧化速率过慢，会使沉铁后液铁离子含量超标，达不到需要的除铁效果。而氧化反应速率除了受催化剂Cu²⁺含量的影响以外，主要由通入反应器的氧气流量来控制。因此，获取溶解氧浓度是对沉铁过程进行机理建模的重要的一环。可是实际生产过程中由于检测难度高等原因，对沉铁溶液进行离线抽样检测并没有检测溶解氧浓度。且氧气溶解度受到多种因素影响，机理模型难以在考虑众多影响因素的同时兼顾可计算性和精确度。由于现有的沉铁过程出口离子浓度预测模型中难以检测溶解氧浓度，造成预测的出口离子浓度精确度不高，使得根据预测出口离子调节氧气流量的控制行为，并不能实现合理控制Fe²⁺的氧化速率，达到高效除铁的作用。

因此，现有的沉铁过程出口离子浓度预测模型中难以检测溶解氧浓度，造成预测的出口离子浓度精确度不高已成为本领域技术人员亟待解决的技术问题。

发明内容

本发明目的在于公开一种沉铁过程出口离子预测方法及系统，以解决现有的沉铁过程出口离子浓度预测模型中难以检测溶解氧浓度，造成预测的出口离子浓度精确度不高的技术问题。

为达上述目的，本发明公开一种沉铁过程出口离子预测方法，包括以下步骤：

S1：构建影响沉铁过程出口离子浓度效果的特征向量，获取有标签数据样本集，采用仿射聚类传播算法得到c个有代表性的中心样本点，分别围绕c个中心样本选取相似度高的L个时间连续样本构成c个子训练样本集；

S2：将对应不同工况的各子训练样本集分别进行局部邻域标准化处理，然后以处理后的各子训练样本集数据分别对c个动态概率隐变量子模型参数进行辨识，以及根据EM算法求得各动态概率隐变量子模型的参数，所述EM算法利用现有的子训练样本集数据信息并使用矩估计对初值进行选择；

S3：基于S1所构建的特征向量获取用于预测离线沉铁过程出口离子浓度特征的无标签的测试样本数据，对同一测试样本数据进行局部邻域标准化处理后输入到各动态概率隐变量子模型得到出口离子浓度预测值的概率分布；

S4:采用贝叶斯决策进行联合估计将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，得到对应所述测试样本数据最终的出口离子浓度的预测值。

优选地，上述S1所构建的特征向量x_i包括：

其中，从左到右分别为：二价铁离子浓度、三价铁离子浓度、铜离子浓度、锌离子浓度、液体流速、FeOOH的质量、氧化锌的质量和通入的氧气的流量；

预测的沉铁过程出口离子浓度向量y_i具体为：

其中，下标中的in和out分别代表输入和输出，为输出的氧气的浓度，为输出的氢离子浓度。

优选地，S1具体包括：

对数据集{X,Y}＝{(x_i,y_i)}_{i＝1,2,...,n}进行仿射传播聚类得到含有c个最优聚类中心v_j＝{x_j,y_j}_{j＝1,2,....,c}的聚类集合V，其中，n为训练数据集的数量；

分别根据聚类中心v_j＝{x_j,y_j}_{j＝1,2,....,c}利用相似度选择算法在数据集{X,Y}＝{(x_i,y_i)}_{i＝1,2,...,n}中选择相似度高的L个连续样本集，用于子模型参数辨识，记做{X_i,Y_i}_{i＝1,2,...,c},其中{X_i,Y_i}＝{(x_j,y_j)}_{j＝1,2,...,L}，对第j个聚类中心v_j，计算出距离每个测试数据集数据的距离：

其中，d_j,s为聚类中心与测试数据集的距离；Δ_s为有标签样本输入数据的方差组成的对角阵的逆矩阵；x_s为样本中的输入特征数据；选择最大的d_j,s距离，然后根据相似度函数选择L个连续的子样本数据集，c个中心点将得到子样本数据集为：{X_i,Y_i}_{i＝1,2,...,c}。

可选的，本发明局部邻域标准化处理的计算公式为：

其中，局部领域N_k(x_i)表示第k个动态概率隐变量子模型数据集，即N_k(x_i)＝{X_i,Y_i}_i＝k,m(N_k(x_i))和s(N_k(x_i))分别表示样本x_i的领域内的均值和标准差，z_i表示x_i进行局部邻域标准化后的值。

优选地，S2具体包括：

以H_t∈R^K表示t时刻描述动态过程的隐变量，H_t-1为t-1时刻的隐变量，并记t时刻的输入变量x_t∈R^D，输出变量y_t∈R^M，则x_t、y_t可由隐变量H_t线性表示，相邻时刻间隐变量的关系及输入输出变量与隐变量的关系表示为：

H_tk＝A_kH_t-1+η_k

x_t＝P_kH_tk+e_k

y_t＝C_kH_tk+w_k

其中，A_k∈R^K×K为第k个动态概率隐变量子模型的不同时刻隐变量间的传递矩阵，K为隐变量的维度，P_k∈R^D×K是输入变量的负载矩阵，C_k∈R^M×K为输出变量的负载矩阵，D和M分别为输入向量和输出向量的维度；η_k是第k个动态概率隐变量子模型动态过程的噪声，e_k和w_k分别为输入输出的噪声项；

基于隐变量H_tk∈R^K先验分布服从高斯分布，得到第k个动态概率隐变量子模型的参数集为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}_{k＝1,2,...,c}，μπ_k和Σ_πk分别为初始时刻，即t＝1时隐变量的均值和方差，满足：H_1k～N(μ_πk,Σ_πk)；

根据得到的第k个工况样本集X_k，即给定输入变量序列x₁,x₂,...,x_L，用x_1:L简要表示，相应地有输出变量的序列y₁,y₂,...,y_L，用y_1:L简要表示，则隐变量序列可用H_1:Lk表示，L表示序列的时间长度；则第k个动态概率隐变量子模型的输入变量、输出变量和隐变量的在完整数据集下的对数似然函数J_k如下：

采用EM算法求解似然函数J_k最大时候的参数值Θ_k，包括求取子训练样本集的对数似然函数J_k关于隐变量分布p(H_1:Lk|x_1:L,y_1:L)的期望：

其中，第k个动态概率隐变量子模型对应的参数为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}，以下将简写为

利用现有的子训练样本集数据信息{X_i,Y_i}＝{(x₁,y₁),(x₂,y₂),......,(x_L,y_L)}，使用矩估计对EM算法的初值进行选择，数据的一阶矩和二阶矩如下：

E_H(H_tk)＝A_k*E_H(H_t-1k) E_H(x_t)＝P_k*E(H_tk)

E_H(yt)＝C_k*E_H(H_tk)

E_H(xtxt^T)＝P_k*E_H(H_tkH_tk ^T)*P_k ^T+Σ_ek E_H(ytyt^T)＝C_k*E_H(H_tkH_tk ^T)C_k ^T+Σ_fk

得到第k个动态概率隐变量子模型EM算法参数的初始值如下：

A_k＝eye(K)；Σ_ηk＝zeros(K,K)；

μπ_k＝normrnd(0,1,[K,1])；Σ_πk＝normrnd(0,0.5,[K,K])；

P_k＝rand(D,K)；C_k＝rand(M,K)；

Σ_fk＝E_H(y_ty_t ^T)-C_k*(Σ_πk+μ_πk*μ_πk ^T)C_k ^T

Σ_ηk＝E_H(x_tx_t ^T)-P_k*(Σ_πk+μ_πk*μ_πk ^T)*P_k ^T

利用EM算法的M步求解出优化的模型参数为：

μ_πk＝E_H(H_1k)＝m_1k

为求解优化的参数，EM算法的E步需要求得隐变量的一阶和二阶条件数学期望，E(H_tk|x_1:L,y_1:L),E(H_tkH_tk ^T|x_1:L,y_1:L)和E(H_tkH_t-1k ^T|x_1:L,y_1:L)，首先需要利用前向算法求得p(H_tk|x_1:t,y_1:t)的分布，设为N(l_tk,L_tk)，参数可由下式得到：

当t＝1时有p(H₁|x₁,y₁)＝N(l_1k,L_1k)，参数可由下式得到：

当t≥2时有p(H_tk|x_1:t,y_1:t)＝N(l_tk,L_tk)，参数可由下式得到：

然后使用后向平滑算法推导出每个时刻隐变量的后验概率分布H_tk|x_1:L,y_1:L～N(m_t,M_t)，得到：

m_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1(m_t+1k-A_kl_tk)+l_tk

M_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1[M_t+1k(A_kL_tkA_k ^T+Σ_ηk)^-1-I]A_kL_tk+L_tk

类似地，当t＝L时，p(H_tk|x_1:L,y_1:L)的平滑均值为m_Lk＝l_Lk，平滑方差为M_Lk＝L_Lk；

最后，根据以上概率分布可得到M步所需的隐变量的一阶和二阶期望如下：

E(H_tk|x_1:L,y_1:L)＝m_tk

通过EM算法反复迭代E步和M步，直到参数矩阵Θ_k收敛，即辨识得到最优的动态概率隐变量子模型参数。

优选地，S3具体包括：

首先需求出初始时刻，即t＝1时的测试样本隐变量后验分布

得到初始时刻测试样本隐变量概率函数的均值和方差分别为：

在已知输入变量序列x_1:t-1的条件下，t≥2，隐变量与输入的联合分布服从高斯分布，形式如同下式所示：

类似地，得到参数：

根据贝叶斯法则可得到，隐变量关于测试样本输入序列的后验分布的方差和均值分别为：

当t＝1时，

当t≥2时，

其中，下标m表示为测试样本的数量，上标q表示测试样本；最后，第t个测试样本序列对应的第k个动态概率隐变量子模型输出的最佳估计为：

优选地，得到各动态概率隐变量子模型的预测输出为{y_1k ^q,y_2k ^q,...,y_mk ^q}_{k＝1,2,...,c}后，S4采用BCM算法将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，具体包括：

对新的预测样本x_t ^q，每个子模型都产生一个高斯分布输出，对于第k个子模型，输出的均值为E(y_tk ^q|X_k,Y_k,x_t ^q)，其协方差为cov(y_tk ^q|X_k,x_t ^q)，其公式如下：

通过BCM算法得到的全局输出也是高斯分布的，其均值和方差如下：

得到的均值E(y_t ^q|X,Y,x_t ^q)可以判定为t时刻样本的预测值y_t ^q。

为达上述目的，本发明还公开一种沉铁过程出口离子预测系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述计算机程序时实现上述相对应方法的步骤。

本发明具有以下有益效果：

对于复杂的沉铁工业过程，由于传感器采集到的过程数据与容器的出口离子浓度存在复杂的未知函数关系，不能单单采用机理模型来得到过程函数关系。但是过程能够产生很多包含过程信息的数据，本发明采用数据驱动的软测量建模方法来预测沉铁过程出口离子浓度，充分考虑了沉铁过程存在多工况、过程的非线性、高噪音以及过程存在动态特性等，基于贝叶斯决策混合动态概率隐变量模型的沉铁过程出口离子预测方法，通过工业过程数据得到沉铁过程出口离子预测模型，提高模型的可信度。相比现有技术而言，使用本发明中的沉铁过程出口离子预测模型预测出的出口离子浓度准确性更高，能达到高效除铁的作用。

下面将参照附图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明实施例的沉铁过程出口离子预测方法流程示意图。

图2是对沉铁过程进行机理分析采用的单个反应器的CSTR系统示意图。

图3是本发明实施例预测出口离子浓度的动态概率隐变量模型结构图。

图4是本发明实施例BCM算法结构示意图。

图5是本发明实施例方法预测值与真实值的对比示意图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

实施例1

一种沉铁过程出口离子预测方法，如图1所示，包括以下步骤：

S1：构建影响沉铁过程出口离子浓度效果的特征向量，获取有标签数据样本集，采用仿射聚类传播算法得到c个有代表性的中心样本点，分别围绕c个中心样本选取相似度高的L个时间连续样本构成c个子训练样本集。

以中国某锌冶炼企业中实际生产数据为例，沉铁过程中的测量变量及其采样频率如表1所示。

表1：

测量变量	采样频率
		每个反应器溶液的流速(包括返流流速,m3/h)	在线检测
每个反应器出口溶液Fe2+，Fe3+，Cu2+浓度(g/L)	每2小时离线化验一次
		每个反应器出口溶液Zn2+浓度(g/L)	每8小时离线化验一次
每个反应器出口溶液H+浓度(pH)	在线检测
		每个反应器添加的氧化锌(t/h)	在线检测
每个反应器添加的氧气流量(m3/h)	在线检测

气体溶于液体的过程就是气体分子转移到液相主体的过程。影响气体在液体中的溶解度的因素通常有压强、温度、液体流速、溶液中离子浓度、溶液中固体含量、机械搅拌等。压强和温度变化时，液相主体间隙的变化使气体溶解度改变；通入氧气速率、溶液流速、机械搅拌通过影响与液面接触的氧气分子数目，改变溶解氧浓度。溶液中离子浓度引起溶液极性变化，氧气作为非极性分子的溶解度会发生改变。此外，沉铁过程中存在反应的动态平衡，其中针铁矿的存在对氧气溶解度也有影响。

如图2所示，以单个反应器为对象，分析针铁矿法沉铁过程影响氧气溶解过程的因素。设分别是反应器入口溶液中Fe²⁺浓度，Fe³⁺浓度，H⁺浓度，为反应器出口溶液中的离子浓度。其中发生的主要反应包括：

氧化反应：4Fe²⁺+4H⁺+O₂→4Fe³⁺+2H₂O

水解反应：

中和反应：2H⁺+ZnO→Zn²⁺+H₂O

实际沉铁过程中，压强和温度的影响忽略不计，机械搅拌速度恒定。小试法中的溶解氧浓度并不都属于同一时刻，也并非每个时刻都能采到，但存在包含上一采样时刻溶解氧浓度的过程数据及同时包含当前采样时刻和上一采样时刻溶解氧浓度的过程数据。因此，为了考虑过程的动态特性，还利用其前采样时刻的输入变量来预测当前采样时刻的输出变量。因此，浸出液中溶解氧浓度的动态模型可以描述成：

y_t＝f(x_t,x_t-1,y_t-1)

在工业过程中采集连续n个有标签的数据集{X,Y}＝{(x_t,y_t)}_{t＝1,2,...,n}。

之后，用仿射聚类传播算法对数据集{X,Y}＝{(x_t,y_t)}_{t＝1,2,...,n}聚类，找到各聚类中心。

仿射传播聚类算法可以描述为：在初始时刻，每个样本点都被当成可能存在的聚类中心，然后通过相似度计算，不断进行聚类过程循环迭代，每个数据点去竞争当选聚类中心，直到找到最优的聚类中心为止。

基于上述采集得到的数据集{X}＝{(x_t)}_{t＝1,2,...,n}，第一步先算出各数据点的相似度矩阵S_n*n(对于x_i和x_j,S_i,j＝-||x_i-x_j||²)，然后不断地从数据点为x_i候选的聚类中心v_j搜集证据，同时也从候选的聚类中心v_j中为数据点x_i搜寻证据，这两个证据参数的定义分别为：

R(i,j)＝S(i,j)-max{A(i,k)+S(i,k)}(k＝1,...,n,k≠j)

A(i,j)＝min{0,R(j,J)+∑{max(0,R(k,j))}}(k＝1,...,n,k≠i,j)

其中，证据可信度R(i,j)用来描述x_i适合作为x_j的聚类中心的程度，证据可用度用来表示x_i选择x_j作为其聚类中心的适合程度。x_j数据点成为聚类中心的可能性取决于R(i,j)和A(i,j)的大小。R(i,j)和A(i,j)越大，那么，x_j成为聚类中心的可能性也就越大。在原来的R(i,j)和A(i,j)的基础上，它们的更新公式如下：

R^new(i,j)＝(1-λ)R^new(i,j)+λR^old(i,j)

A^new(i,j)＝(1-λ)A^new(i,j)+λA^old(i,j)

仿射聚类传播算法通过以上公式进行循环更新，对各个证据进行搜索、传递和整理，λ来控制算法的迭代更新速度，对于数据点x_i，如果数据点x_j能够是使R(i,j)+A(i,j)成为R(i,k)+A(i,k)，k＝1,2，…,n最大，那么认为x_j是x_i的聚类中心。通过这样迭代竞争的方式完成整个聚类过程，可以得到含有个最优聚类中心v_j＝{x_j,y_j}_{j＝1,2,...,c}的聚类集合。

利用聚类中心在数据集{X,Y}中找到样本大小为L的连续子数据集{X_i,Y_i}_i＝1,2,..,t。

对每个聚类中心v_j＝{x_j,y_j}_{j＝1,2,....,c}利用相似度选择函数在测试数据{X,Y}＝{(x_i,y_i)}_{i＝1,2,...,n}选择相似度高的连续样本集{X_i,Y_i}_{i＝1,2,...,c}，其中，{X_i,Y_i}＝{(x_j,y_j)}_{j＝1,2,...,L}，L为选择的相似样本数，对第j个聚类中心，首先根据聚类中心v_j，算出距离每个测试数据集数据的距离，可由下式得到：

其中，d_j,s为聚类中心与测试数据集的距离计算公式；Δ_s为有标签样本输入数据的方差组成的对角阵的逆矩阵；x_s为样本中的输入特征数据，v_j为第j个聚类中心，n表示测试数据集样本个数。选择最大的d_j,s距离，然后根据相似度函数选择L个连续的子样本数据集，c个聚类中心得到的子样本数据集为，{X_i,Y_i}_{i＝1,2,...,c}，为后续建立c个动态概率隐变量子模型提供训练数据集。

在该步骤中，优选的局部邻域标准化处理的计算公式为：

其中，局部领域N_k(x_i)表示第k个动态概率隐变量子模型数据集，即N_k(x_i)＝{X_i,Y_i}_i＝k,m(N_k(x_i))和s(N_k(x_i))分别表示样本x_i的领域内的均值和标准差，z_i表示x_i进行局部邻域标准化后的值。

S2：将对应不同工况的各子训练样本集分别进行局部邻域标准化处理，然后以处理后的各子训练样本集数据分别对c个动态概率隐变量子模型参数进行辨识，以及根据EM算法求得各子动态概率隐变量子模型的参数，所述EM算法利用现有的子训练样本集数据信息并使用矩估计对初值进行选择。

优选地，参考图3所示的预测出口离子浓度的动态概率隐变量模型结构，以H_t∈R^K表示t时刻描述动态过程的隐变量，H_t-1为t-1时刻的隐变量，并记t时刻的输入变量x_t∈R^D，输出变量y_t∈R^M，则x_t、y_t可由隐变量H_t线性表示，相邻时刻间隐变量的关系及输入输出变量与隐变量的关系表示为：

H_tk＝A_kH_t-1+η_k

x_t＝P_kH_tk+e_k

y_t＝C_kH_tk+w_k

其中，A_k∈R^K×K为第k个动态概率隐变量子模型的不同时刻隐变量间的传递矩阵，K为隐变量的维度，P_k∈R^D×K是输入变量的负载矩阵，C_k∈R^M×K为输出变量的负载矩阵，D和M分别为输入向量和输出向量的维度；η_k是第k个动态概率隐变量子模型动态过程的噪声，e_k和w_k分别为输入输出的噪声项；

基于隐变量H_tk∈R^K先验分布服从高斯分布，得到第k个动态概率隐变量子模型的参数集为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}_{k＝1,2,...,c}，μπ_k和Σ_πk分别为初始时刻，即t＝1时隐变量的均值和方差，满足：H_1k～N(μ_πk,Σ_πk)；

根据得到的第k个工况样本集X_k，即给定输入变量序列x₁,x₂,...,x_L，用x_1:L简要表示，相应地有输出变量的序列y₁,y₂,...,y_L，用y_1:L简要表示，则隐变量序列可用H_1:Lk表示，L表示序列的时间长度；则第k个动态概率隐变量子模型的输入变量、输出变量和隐变量的在完整数据集下的对数似然函数J_k如下：

采用EM算法求解似然函数J_k最大时候的参数值Θ_k，包括求取子训练样本集的对数似然函数J_k关于隐变量分布p(H_1:Lk|x_1:L,y_1:L)的期望：

其中，第k个动态概率隐变量子模型对应的参数为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}，以下将简写为

利用现有的子训练样本集数据信息{X_i,Y_i}＝{(x₁,y₁),(x₂,y₂),......,(x_L,y_L)}，使用矩估计对EM算法的初值进行选择，数据的一阶矩和二阶矩如下：

E_H(H_tk)＝A_k*E_H(H_t-1k) E_H(x_t)＝P_k*E(H_tk)

E_H(y_t)＝C_k*E_H(H_tk)

E_H(x_tx_t ^T)＝P_k*E_H(H_tkH_tk ^T)*P_k ^T+Σ_ek E_H(y_ty_t ^T)＝C_k*E_H(H_tkH_tk ^T)C_k ^T+Σ_fk

得到第k个动态概率隐变量子模型EM算法参数的初始值如下：

A_k＝eye(K)；Σ_ηk＝zeros(K,K)；

μ_πk＝normrnd(0,1,[K,1])；Σ_πk＝normrnd(0,0.5,[K,K])；

P_k＝rand(D,K)；C_k＝rand(M,K)；

Σ_fk＝E_H(y_ty_t ^T)-C_k ^*(Σ_πk+μ_πk ^*μ_πk ^T)C_k ^T

Σ_ηk＝E_H(x_tx_t ^T)-P_k*(Σ_πk+μ_πk*μ_πk ^T)*P_k ^T

利用EM算法的M步求解出优化的模型参数为：

μ_πk＝E_H(H_1k)＝m_1k

为求解优化的参数，EM算法的E步需要求得隐变量的一阶和二阶条件数学期望，E(H_tk|x_1:L,y_1:L),E(H_tkH_tk ^T|x_1:L,y_1:L)和E(H_tkH_t-1k ^T|x_1:L,y_1:L)，首先需要利用前向算法求得p(H_tk|x_1:t,y_1:t)的分布，设为N(l_tk,L_tk)，参数可由下式得到：

当t＝1时有p(H₁|x₁,y₁)＝N(l_1k,L_1k)，参数可由下式得到：

当t≥2时有p(H_tk|x_1:t,y_1:t)＝N(l_tk,L_tk)，参数可由下式得到：

然后使用后向平滑算法推导出每个时刻隐变量的后验概率分布H_tk|x_1:L,y_1:L～N(m_t,M_t)，得到：

m_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1(m_t+1k-A_kl_tk)+l_tk

M_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1[M_t+1k(A_kL_tkA_k ^T+Σ_ηk)^-1-I]A_kL_tk+L_tk

类似地，当t＝L时，p(H_tk|x_1:L,y_1:L)的平滑均值为m_Lk＝l_Lk，平滑方差为M_Lk＝L_Lk；

最后，根据以上概率分布可得到M步所需的隐变量的一阶和二阶期望如下：

E(H_tk|x_1:L,y_1:L)＝m_tk

在该步骤中，由于传统的EM算法求解的不是全局最优解，故本实施例对其进行了改进，通过利用现有的子训练样本集数据信息并使用矩估计对初值进行选择，有效避免了陷入局部最优。

通过EM算法反复迭代E步和M步，直到参数矩阵Θ_k收敛，即辨识得到最优的动态概率隐变量子模型参数。

S3：基于S1所构建的特征向量获取用于预测离线沉铁过程出口离子浓度特征的无标签的测试样本数据，对同一测试样本数据进行局部邻域标准化处理后输入到各动态概率隐变量子模型得到出口离子浓度预测值的概率分布。

优选地，该步骤具体包括：

首先需求出初始时刻，即t＝1时的测试样本隐变量后验分布

得到初始时刻测试样本隐变量概率函数的均值和方差分别为：

在已知输入变量序列x_1:t-1的条件下，t≥2，隐变量与输入的联合分布服从高斯分布，形式如同下式所示：

类似地，得到参数：

根据贝叶斯法则可得到，隐变量关于测试样本输入序列的后验分布的方差和均值分别为：

当t＝1时，

当t≥2时，

其中，下标m表示为测试样本的数量，上标q表示测试样本；最后，第t个测试样本序列对应的第k个动态概率隐变量子模型输出的最佳估计为：

S4:采用贝叶斯决策进行联合估计将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，得到对应所述测试样本数据最终的出口离子浓度的预测值。

得到各动态概率隐变量子模型的预测输出为{y_1k ^q,y_2k ^q,...,y_mk ^q}_{k＝1,2,...,c}后，该步骤采用BCM算法，其对应的结构图参照图4，将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，具体包括：

对新的预测样本x_t ^q，每个子模型都产生一个高斯分布输出，对于第k个子模型，输出的均值为E(y_tk ^q|X_k,Y_k,x_t ^q)，其协方差为cov(y_tk ^q|X_k,x_t ^q)，其公式如下：

通过BCM算法得到的全局输出也是高斯分布的，其均值和方差如下：

得到的均值E(y_t ^q|X,Y,x_t ^q)可以判定为t时刻样本的预测值y_t ^q。

基于本实施例贝叶斯决策混合动态概率隐变量模型的沉铁过程出口离子预测方法构造模型预测值与真实值对比图参照图5；由图可知，本实施例预测值与实测值基本吻合，准确度可靠。

实施例2

一种沉铁过程出口离子预测系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述计算机程序时实现上述方法所对应的各步骤。

综上，本发明上述各实施例所分别公开的沉铁过程出口离子预测方法及系统，至少具有以下有益效果：

对于复杂的沉铁工业过程，由于传感器采集到的过程数据与容器的出口离子浓度存在复杂的未知函数关系，不能单单采用机理模型来得到过程函数关系。但是过程能够产生很多包含过程信息的数据，本发明采用数据驱动的软测量建模方法来预测沉铁过程出口离子浓度，充分考虑了沉铁过程存在多工况、过程的非线性、高噪音以及过程存在动态特性等，基于贝叶斯决策混合动态概率隐变量模型的沉铁过程出口离子预测方法，通过工业过程数据得到沉铁过程出口离子预测模型，提高模型的可信度。相比现有技术而言，使用本发明中的沉铁过程出口离子预测模型预测出的出口离子浓度准确性更高，能达到高效除铁的作用。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种沉铁过程出口离子预测方法，其特征在于，包括以下步骤：

S1：构建影响沉铁过程出口离子浓度效果的特征向量，获取有标签数据样本集，采用仿射聚类传播算法得到c个有代表性的中心样本点，分别围绕c个中心样本选取相似度高的L个时间连续样本构成c个子训练样本集；其中,所构建的特征向量x_i包括：

其中，从左到右分别为：二价铁离子浓度、三价铁离子浓度、铜离子浓度、锌离子浓度、液体流速、FeOOH的质量、氧化锌的质量和通入的氧气的流量；

预测的沉铁过程出口离子浓度向量y_i具体为：

其中，下标中的in和out分别代表输入和输出，为输出的氧气的浓度，为输出的氢离子浓度；

S2：将对应不同工况的各子训练样本集分别进行局部邻域标准化处理，然后以处理后的各子训练样本集数据分别对c个动态概率隐变量子模型参数进行辨识，以及根据EM算法求得各子动态概率隐变量子模型的参数，所述EM算法利用现有的子训练样本集数据信息并使用矩估计对初值进行选择；

S3：基于S1所构建的特征向量获取用于预测离线沉铁过程出口离子浓度特征的无标签的测试样本数据，对同一测试样本数据进行局部邻域标准化处理后输入到各动态概率隐变量子模型得到出口离子浓度预测值的概率分布；

S4:采用贝叶斯决策进行联合估计将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，得到对应所述测试样本数据最终的出口离子浓度的预测值。

2.根据权利要求1所述的沉铁过程出口离子预测方法，其特征在于，S1具体包括：

对数据集{X,Y}＝{(x_i,y_i)}_{i＝1,2,...,n}进行仿射传播聚类得到含有c个最优聚类中心v_j＝{x_j,y_j}_{j＝1,2,....,c}的聚类集合V，其中，n为训练数据集的数量；

分别根据聚类中心v_j＝{x_j,y_j}_{j＝1,2,....,c}利用相似度选择算法在数据集{X,Y}＝{(x_i,y_i)}_{i＝1,2,...,n}中选择相似度高的L个连续样本集，用于子模型参数辨识，记做{X_i,Y_i}_{i＝1,2,...,c},其中{X_i,Y_i}＝{(x_j,y_j)}_{j＝1,2,...,L}，对第j个聚类中心v_j，计算出距离每个测试数据集数据的距离：

其中，d_j,s为聚类中心与测试数据集的距离；Δ_s为有标签样本输入数据的方差组成的对角阵的逆矩阵；x_s为样本中的输入特征数据；选择最大的d_j,s距离，然后根据相似度函数选择L个连续的子样本数据集，c个中心点将得到子样本数据集为：{X_i,Y_i}_{i＝1,2,...,c}。

3.根据权利要求2所述的沉铁过程出口离子预测方法，其特征在于，局部邻域标准化处理的计算公式为：

其中，局部领域N_k(x_i)表示第k个动态概率隐变量子模型数据集，即N_k(x_i)＝{X_i,Y_i}_i＝k,m(N_k(x_i))和s(N_k(x_i))分别表示样本x_i的领域内的均值和标准差，z_i表示x_i进行局部邻域标准化后的值。

4.根据权利要求3所述的沉铁过程出口离子预测方法，其特征在于，S2具体包括：

以H_t∈R^K表示t时刻描述动态过程的隐变量，H_t-1为t-1时刻的隐变量，并记t时刻的输入变量x_t∈R^D，输出变量y_t∈R^M，则x_t、y_t可由隐变量H_t线性表示，相邻时刻间隐变量的关系及输入输出变量与隐变量的关系表示为：

H_tk＝A_kH_t-1+η_k

x_t＝P_kH_tk+e_k

y_t＝C_kH_tk+w_k

其中，A_k∈R^K×K为第k个动态概率隐变量子模型的不同时刻隐变量间的传递矩阵，K为隐变量的维度，P_k∈R^D×K是输入变量的负载矩阵，C_k∈R^M×K为输出变量的负载矩阵，D和M分别为输入向量和输出向量的维度；η_k是第k个动态概率隐变量子模型动态过程的噪声，e_k和w_k分别为输入输出的噪声项；

基于隐变量H_tk∈R^K先验分布服从高斯分布，得到第k个动态概率隐变量子模型的参数集为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}_{k＝1,2,...,c}，μ_πk和Σ_πk分别为初始时刻，即t＝1时隐变量的均值和方差，满足：H_1k:N(μ_πk,Σ_πk)；

根据得到的第k个工况样本集X_k，即给定输入变量序列x₁,x₂,...,x_L，用x_1:L简要表示，相应地有输出变量的序列y₁,y₂,...,y_L，用y_1:L简要表示，则隐变量序列可用H_1:Lk表示，L表示序列的时间长度；则第k个动态概率隐变量子模型的输入变量、输出变量和隐变量的在完整数据集下的对数似然函数J_k如下：

采用EM算法求解似然函数J_k最大时候的参数值Θ_k，包括求取子训练样本集的对数似然函数J_k关于隐变量分布p(H_1:Lk|x_1:L,y_1:L)的期望：

其中，第k个动态概率隐变量子模型对应的参数为Θ_k＝{μ_πk,Σ_πk,A_k,Σ_ηk,P_k,C_k,Σ_ek,Σ_wk}，以下将简写为

利用现有的子训练样本集数据信息{X_i,Y_i}＝{(x₁,y₁),(x₂,y₂),......,(x_L,y_L)}，使用矩估计对EM算法的初值进行选择，数据的一阶矩和二阶矩如下：

E_H(H_tk)＝A_k*E_H(H_t-1k) E_H(x_t)＝P_k*E(H_tk)

E_H(y_t)＝C_k*E_H(H_tk) E_H(H_tkH_tk ^T)＝A_k*E_H(H_t-1kH_t-1k ^T)A_k ^T+Σ_ηk

E_H(x_tx_t ^T)＝P_k*E_H(H_tkH_tk ^T)*P_k ^T+Σ_ek E_H(y_ty_t ^T)＝C_k*E_H(H_tkH_tk ^T)C_k ^T+Σ_fk

得到第k个动态概率隐变量子模型EM算法参数的初始值如下：

A_k＝eye(K)；Σ_ηk＝zeros(K,K)；

μ_πk＝normrnd(0,1,[K,1])；Σ_πk＝normrnd(0,0.5,[K,K])；

P_k＝rand(D,K)；C_k＝rand(M,K)；

Σ_fk＝E_H(y_ty_t ^T)-C_k*(Σ_πk+μ_πk*μ_πk ^T)C_k ^T

Σ_ηk＝E_H(x_tx_t ^T)-P_k*(Σ_πk+μ_πk*μ_πk ^T)*P_k ^T

利用EM算法的M步求解出优化的模型参数为：

μ_πk＝E_H(H_1k)＝m_1k

为求解优化的参数，EM算法的E步需要求得隐变量的一阶和二阶条件数学期望，E(H_tk|x_1:L,y_1:L),E(H_tkH_tk ^T|x_1:L,y_1:L)和E(H_tkH_t-1k ^T|x_1:L,y_1:L)，首先需要利用前向算法求得p(H_tk|x_1:t,y_1:t)的分布，设为N(l_tk,L_tk)，参数可由下式得到：

当t＝1时有p(H₁|x₁,y₁)＝N(l_1k,L_1k)，参数可由下式得到：

当t≥2时有p(H_tk|x_1:t,y_1:t)＝N(l_tk,L_tk)，参数可由下式得到：

然后使用后向平滑算法推导出每个时刻隐变量的后验概率分布H_tk|x_1:L,y_1:L:N(m_t,M_t)，得到：

m_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1(m_t+1k-A_kl_tk)+l_tk

M_tk＝L_tkA_k ^T(A_kL_tkA_k ^T+Σ_ηk)^-1[M_t+1k(A_kL_tkA_k ^T+Σ_ηk)^-1-I]A_kL_tk+L_tk

类似地，当t＝L时，p(H_tk|x_1:L,y_1:L)的平滑均值为m_Lk＝l_Lk，平滑方差为M_Lk＝L_Lk；

最后，根据以上概率分布可得到M步所需的隐变量的一阶和二阶期望如下：

E(H_tk|x_1:L,y_1:L)＝m_tk

通过EM算法反复迭代E步和M步，直到参数矩阵Θ_k收敛，即辨识得到最优的动态概率隐变量子模型参数。

5.根据权利要求4所述的沉铁过程出口离子预测方法，其特征在于，S3具体包括：

首先需求出初始时刻，即t＝1时的测试样本隐变量后验分布

得到初始时刻测试样本隐变量概率函数的均值和方差分别为：

在已知输入变量序列x_1:t-1的条件下，t≥2，隐变量与输入的联合分布服从高斯分布，形式如同下式所示：

类似地，得到参数：

根据贝叶斯法则可得到，隐变量关于测试样本输入序列的后验分布的方差和均值分别为：

当t＝1时，

当t≥2时，

其中，下标m表示为测试样本的数量，上标q表示测试样本；最后，第t个测试样本序列对应的第k个动态概率隐变量子模型输出的最佳估计为：

6.根据权利要求5所述的沉铁过程出口离子预测方法，其特征在于，得到各动态概率隐变量子模型的预测输出为{y_1k ^q,y_2k ^q,...,y_mk ^q}_{k＝1,2,...,c}后，S4采用BCM算法将c个动态概率隐变量子模型输出的预测值的概率分布结果进行融合，具体包括：

对新的预测样本x_t ^q，每个子模型都产生一个高斯分布输出，对于第k个子模型，输出的均值为E(y_tk ^q|X_k,Y_k,x_t ^q)，其协方差为cov(y_tk ^q|X_k,x_t ^q)，其公式如下：

通过BCM算法得到的全局输出也是高斯分布的，其均值和方差如下：

得到的均值E(y_t ^q|X,Y,x_t ^q)可以判定为t时刻样本的预测值y_t ^q。

7.一种沉铁过程出口离子预测系统，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述权利要求1至6任一所述方法的步骤。