基于粒子滤波和EM算法的PX氧化副反应尾气COx含量软测量的研究

© 2013年。出版社：Elsevier B.V.信息工程研究院负责评选和同行评议可在www.sciencedirect.com在线获取ScienceDirectIERI Procedia 6（2014）63 - 712013未来软件工程与多媒体工程国际会议基于粒子滤波和EM算法的PX氧化副反应尾气COx含量软测量刘冉冉a，b，李正明a， *a江苏大学电气与信息工程学院，江苏镇江212013b常州刘国军职业技术学院，江苏常州213004摘要基于期望最大化算法，提出了数据驱动非线性模型的参数估计方法。在此基础上，利用粒子滤波近似计算积分，推导出基于粒子滤波的EM算法。仿真结果验证了该算法在PX氧化副反应尾气COx含量软测量中的有效性© 2014年。李正明出版社：Elsevier B.V. 在CC BY-NC-ND许可下开放访问。信息工程研究所科学委员会负责同行评议关键词：期望最大化算法;粒子滤波;数据驱动;软测量;1. 介绍我们一直致力于通过对工业过程的有效控制来改善过程设备的性能，从而实现生产率的提高和能源消耗的降低。然而，在实际生产中，由于技术或经济的原因，一些关键参数* 李正明。联系电话：.+ 0-511-887-90618;传真：.+0-511-887-91245.电子邮件地址：lzming@ujs.edu.cn2212-6678 © 2014。李正明。由Elsevier B. V.在CC BY-NC-ND许可下开放获取。 信息工程研究所科学委员会负责同行评审doi：10.1016/j.ieri.2014.03.01164刘然然和李正明/ IERI Procedia 6（2014）63生物发酵罐的浓度、高炉铁水中的硅含量以及化学反应器中的反应物浓度、转化率和催化剂活性。两种传统的解决方案之一是使用在线分析仪控制。该方法虽然能实现快速检测，但设备投资大、维护费用高、滞后大，影响了慢化效果。另一种方法是离线实验分析，以保证测量数据的准确性。但该方法测量时间长，不利于实时处理。因此，软测量技术应运而生。软测量技术通过建立若干数学关系，选取几个容易测量的变量（或称辅助变量），对难以测量或暂时无法测量的重要变量（或称主导变量）进行推断和估计。该方法的实质是建立辅助变量和主导变量之间的关系模型（Brosillow 1978）。软测量已经被广泛应用（P.Kadleca，B. Gabrys和S. Strandtb 2009; S. Khatibisepehr和B. Huang 2008;W. Yan，H.邵和X。Wang 2004）。它在常规检测的基础上，利用辅助变量与主导变量之间的关系，通过软件计算得到主导变量的估计值。软测量方法的研究和应用较多，按建模方法分为机理模型、灰箱模型和数据驱动模型三类。虽然数据驱动模型的结构多种多样，但研究最多的是稳态模型和动态模型。考虑到数据随时间的变化，常用的动态模型包括状态空间模型和输入输出模型，例如自回归模型（例如，ARX，ARMAX），输出误差（OE）模型和Box-Jenkins（BJ）模型（P.K. Pearson 2004）。非线性动力学模型代表了上述模型中的一个非线性动力学系统。因此，提出的软测量技术主要采用这样的模型在这项工作中。近年来，基于SMC（顺序蒙特卡罗）（A。Doucet，N. de Freitas和N. Gordon2001），提出了一种结合粒子滤波和EM算法的参数估计方法[8]。在文学（RB。使用EM算法、粒子滤波器和平滑器来近似对数似然函数。使用逐点状态估计，每个密度函数 EM算法每次迭代都要计算状态平滑后的状态，计算量很大。在工作中，期望函数的计算使用粒子滤波。用滤波器代替平滑器，不仅可以减少计算量，而且可以得到可接受的估计性能。第二节介绍了EM算法框架下非线性状态空间模型的参数估计。在第三节中，将所提出的算法应用于PX氧化副反应中，以估算尾气中的COx含量。第四节讨论了几个问题，第五节根据仿真结果得出结论。2. 状态空间模型2.1. 模型结构给出了非线性状态空间模型xf（xt1，ut1， ）tyth（xt，）vt（一）（二）式中，λ是系统参数的矢量，x t、u t和y t分别是输入和输出的状态、测量。t和v t是协方差：分别是Q和R的白噪声。输入u1，u2......，u was known，omitted在的后续推导为简化。让X代表隐藏状态x1，x2......， xT，Y0代表y1，y2，yT，以体现相应的输出测量。刘然然和李正明/ IERI Procedia 6（2014）63650000不2.2. 期望最大化算法在统计学中，期望最大化算法在依赖于不可观测潜变量的概率模型中搜索最大似然参数估计。EM算法交替计算两个步骤：第一步是计算期望值（E），其中潜在变量的现有估计值用于计算其最大似然估计值;第二步是最大化（M），其中步骤E的最大似然估计值用于计算参数值。并且在步骤M中获得的参数的估计被用于下一步骤E，因此该过程交替地继续。对于状态空间模型（1）和（2），步骤E是计算Q函数，定义如下：Q（Q）|（k）X|Y0，阿斯克[p|（）][p |n（）] p（X |Y，k）dX（三）步骤M是确定Q函数最大时的的值，其为最大值Q（Ω）|（克罗地亚）（四）2.3. 基于EM算法对于状态空间模型（1）和（2），设p∈（x1：T，y1：T|表示状态和观测变量的似然函数。将Q函数定义为对数似然函数的期望值，并将Q函数定义为对数似然函数的期望值（ x1：T，y1：T|图）流动积分公式（Jing D e ng 2013）中所示：Q（Q）|（k）X|Y0，阿斯克[p|（）][p |n（）] p（X |Y，k）dX（五）log[p(x1:T，y1：T|)] p (x1:T|y1:T，k）dx在公式（5）中，第一部分p（x1：T，y1：T|是状态和输出的联合密度函数，其可以如下导出：TTp（x1：T，y1：T|）|）|xt 1，）p（yt|xt，）（六）t2t1将公式（6）代入公式（5），得到：Q（）|k）对数[p（x1：T，y1：T|)] p (x1:T|y1:T，k）dx最大值log[p（x1|y1：T，）]p（x1：T|y1:T，k）dx[001 pdf 1st-31files ]|，）]（ |、克罗地亚）（七）t102不p xtxt1px1：Ty1：Tdx1：T [p（y|x，n）] p（x| y ，k）dxt1t t1：T1：T1：T为了计算公式（7）中的Q函数，密度函数的值p（x1：T|y1:T，必须是得到了由于2.4. 基于粒子滤波粒子滤波是对概率密度函数p（x k）进行逼近|z k）通过寻找在状态空间中散布的一组随机样本。用样本均值代替积分运算，得到了1：T1：T1：T66刘然然和李正明/ IERI Procedia 6（2014）63N不第一章1 不tt1最小方差状态估计这些样品被称为在文献[8]和[10]的基础上，利用粒子滤波估计密度函数，计算Q函数。当给定当前参数估计时，状态的密度函数可以近似为：p（x）|y，k）（xxi）（八）t1：t t ttt1其中， 是狄拉克δ函数，N是粒子的数量，i是第i个粒子的标准化权重粒子，和1。因此，先前基于后验密度函数的积分公式（3）可以转换为基于加权样本的求和运算。在许多情况下，上述后验密度可以是多变量的、高维的、多峰的、非标准的和非分析的。从这样的密度函数中直接采样粒子是非常困难的，因此引入了重要性采样（IS）（MS Arulampalam 2002;陈哲2003）。后验概率可以通过对重要性抽样进行重抽样得到。当样本数趋于无穷大时，粒子滤波得到的估计是渐近无偏的（S.D.）。Grantham和L.H. Ungar 1990）。一般情况下，重要性密度取状态转移的概率分布，如：q（x|y ，nk）np（x|X、克罗地亚）（九）t1：t t tt 1此时，每个粒子的权重归一化可以表示如下：：我我tt1p（y |x i，k）：（十）吉吉阿吉提（十一）t1Ni第一章1不然而，该算法有一个主要的缺陷-权值的减少。其直接后果是，大多数颗粒的颗粒浓度的计算，由于其重量小，意义不大。减少这种现象的最有效方法是重新采样。其基本思想是通过对粒子和相应权值所表示的概率密度函数进行重采样来增加权值大的粒子数。在重新采样之后，每个粒子的权重将重置为“0”。tN重采样虽然解决了粒子的减少，但同时也降低了粒子的多样性。因此，重新采样仅在必要时使用，而不是在每个步骤中使用。因此，引入了有效样本量的概念来衡量权重的降低，定义为（A。X. Shao2010）：Neff 1N（NH2）2（十二）1t其中，Rai是归一化的重要性权重。阈值N 可以提前设计。重新采样tth当N eff时需要 第四，权重的减少相对严重。否则就没有必要了。结合当前参数估计，粒子滤波算法的流程图可以概括如下：第一步：从p（x）画{x i} N|k），设t=1。0i100N第二步：重要性抽样。不刘然然和李正明/ IERI Procedia 6（2014）6367不t1tt从p（x）生成{x i} N|X，J.K.）。t i1不t1第三步：重锤。分配权重{i} N使用公式（10）和（11），计算机N使用等式(12)，如果N中国转到步骤4，或转到步骤5。第四步：休息。不第一章1EFFEff th，生成{x i} N，设置 1不第一章1不N.第五步：设置t=t+1，重复第二步到第四步，直到t大于T。2.5. 基于粒子滤波的对p（x ~1：T）| y1:T，k）是所有状态的平滑问题。基于EM算法，计算所有状态的平滑。后来，文学的方法（A。Doucet 2001）被用来进一步边缘化国家。Q函数可以进一步表示为以下形式：Q（）|k）对数[p（x1|y1：T，）] p（x1|y1：T，k）dx最大值log[p（x|X，n）]p（x| y ，k）dx（十三）t102不tt1t1： t1：Tt1：t [p（y|x，n）] p（x|y，k）dxt1t tt1：T t计算p（x|y，k）和p（x，x|y，krk）也是平滑问题，t1：Tt1t1：T的计算。一个可用的解决方案是：当t=1：T时，替换p（x）|y，k），其中p（x|y，J.K.）。当t=1时：T-t1：T t1：t1、更换p（x） ，x|y，，，p（x，x |y，J. Deng和B.Huang 2012）。在实际应用中，t1t1：T不t 11：t粒子滤波算法可以大大降低计算复杂度。该算法通过迭代计算，在一定程度上提高了EM算法的性能。在公式（13）中，通过粒子滤波，|y因此，可以近似为：t1：TNp（x）|y ，nk）np（x|y，k）（xi）（十四）t1：Tt1：t t t t第一章1xt和x1的联合密度函数可以近似为：p（x，x |y ，nk）np（x，x |y，k）不t1 1：Ttt1 1：tp（x|x，k）p（x|y、克罗地亚）（十五）t1t t1：tN阿吉吉（x X（一）哪里不|t 1 t tt 1第一章1t1p（x i | x i, k) w ii（十六）不|t 1N第一章1p（xi）| x i,k) w i使用这些密度函数的近似，Q函数可以写成以下形式：168刘然然和李正明/ IERI Procedia 6（2014）63N TNt tt123456Q（）|k）|（）]log[p（x|x i ，）]11第一章1不Nt2 i1t1|ttt 1（十七） ilog[p (y |xi，）]t2 i1在得到Q函数的近似值后，可采用EM算法进行计算。在步骤期望中，利用当前估计的参数k来估计Q函数。在逐步最大化中，通过最大化Q函数来获得最新的参数1通过计算Q函数对θ的微分，可以得到Q函数最大时的参数值。在每次迭代中，当ΔQ/Δj≤0时，参数估计是最优的，其中Δj是第j个系统参数。基于PF的EM算法（简称PF-EM算法）可以概括为：第一步：初始化参数，t=0。第二步：期待。在t时刻，利用公式（17）计算Q函数的估计，给出当前系统参数Δk。步骤3：最大化。通过最大化Q函数的近似公式，可以获得最新的参数_l令k=k+1。步骤4：重复步骤2和3，直到满足收敛条件。例如，两次迭代之间系统参数的变化小于容差。3. PX氧化副反应作为聚酯纤维和塑料的重要原料，精对苯二甲酸（PTA）是由PX经AMOCO氧化反应制得的，其催化剂为Mn和Co，助催化剂为Br，溶剂为醋酸（HAc）（Prengle Jr H W and Barond N 1970）。除此之外，还有大量的副反应，产物主要是CO和CO2（张春阳1999;云明2002;吴德华2000），导致产品质量下降，催化剂和促进剂消耗剧增本文以扬子石化公司的装置为研究对象，应用该算法对COx含量进行了研究，通过主成分分析，筛选出6个变量作为指标，并对各变量的影响因素进行了分析。输入数据，包括反应温度（u1，？溶剂比（u2，Kg.HAc/Kg.PX），Co浓度（u3，wt%）、Mn浓度（u4，wt%）、Br浓度（u5，wt%）和反应温度（t）。时间（u6，s）。本文中提供的所有工业数据都已归一化。PX氧化副反应是一个连续的过程。并通过动态模型描述了尾气中COX含量与历史值的关系，最终的非线性状态空间模型如下，并证明该模型是可行的：BT u c cos（x ）t t1y t xtt1t1t（十八）其中，待识别的系统参数为[a，b，b，b，b，b，c]T，X t 代表 的 状态向量，y表示输出向量，并且u∈ [u，u，u，u，u，u]T表示输入向量。中的输入数tt1 t2 T3 T4 t5 t6 t这个型号是6。作为独立同分布的高斯噪声，Vt和Vt分别代表过程噪声和观测噪声。它们的协方差分别为Q=0.05和R=0.05。采用基于粒子滤波的EM算法对状态空间模型参数进行估计。在EM算法的步长期望中，根据公式（17）计算Q函数，其中刘然然和李正明/ IERI Procedia 6（2014）6369年实验室分析PF-EM估计误差1 1（xiaxiBuAcos（xi））2log[p（x|Xi，）[exp[tt 1t1t1]]（十九）tt 12Q1 1（y）2log[p（y|x i，）] log[exp[expt]]（二十）t t2R通过使Q函数的导数等于0，新的 新的，b1新，b2新，b3新，b4新，b5新，b6新，c]T可以在迭代期间经由先前步骤的T来确定然后利用本文提出的参数估计方法确定系统参数。粒子滤波器的数量为100，初始系统参数为[0.5，0.1，0.1，0.1，0.1，0.1，1]。利用提出的软测量方法进行了仿真分析。状态真值与PF-EM软测量估计之间的趋势比较如图1所示。平方误差的变化如图2所示。样本空间和后验密度随时间的关系如图3所示。1510.9100.8五○七0.600.5500400300200-5-10-1505101520253035404550时间0.40.30.20.100102030405060708090100时间100060时间0-2020样本空间图1.走势对比图2平方误差图3.样本空间和后验密度图1-3表明，该算法可以有效地识别非线性模型。预测变量能很好地跟踪实际状态的变化。在有限次迭代之后，系统误差趋于0.通过粒子滤波器对重要性采样函数中的粒子进行重采样，可以很好地逼近后验密度。4. 讨论提出了一种非线性状态空间模型软测量的实现方法。从而实现了基于粒子滤波和EM算法的非线性系统参数估计。在软测量应用中，现有的大多数方法都是根据预测值与实际观测值之间的更新偏差进行偏差补偿。因此，该算法具有很大的潜力。然而，为保证实际工业应用的可靠性，该方法仍存在许多问题值得进一步研究，如测试数据的延迟、预测的鲁棒性以及由于更新实验数据而导致的预测值突变等。在线滤波可以减小实验室采样中的时间延迟的影响。预测的鲁棒性可以通过机理模型来提高。通过调节状态空间模型的方差可以减小预测值的突变。因此，模型的预测质量可以通过组合这些技术来增强。5. 结论提出了一种非线性数据驱动状态空间模型的辨识方法。EM算法是2014年4010200-10COx含量后验密度70刘然然和李正明/ IERI Procedia 6（2014）63用于估计系统参数，而粒子滤波器近似期望函数。所提出的PF-EM算法由于计算量小，可以识别实际的工业过程。仿真结果表明，该算法能有效地辨识非线性系统，且估计误差在允许范围内。因此，它可以扩展到其他工业过程。确认作者感谢江苏省高校毕业生科研创新项目（CXLX12_0648）、国家自然科学基金（61273142）和江苏省高等学校重点学科建设项目（PAPD）的资助。引用[1] 布罗西洛过程推理控制Jour. AichE，1978，24（3）：485-509.[2] P. Kadleca，B. Gabrys，S. Strandtb，过程工业中的数据驱动软传感器，计算机和化学工程：33（2009）795[3] S.哈蒂比塞佩尔湾黄文，软测量中不规则数据的处理：贝叶斯方法与比较研究，工业与工程化学研究：47（2008）8713[4] W. Yan，H.邵，X.王文，基于支持向量机和贝叶斯模型选择的软测量建模，计算机与化工：28（2004）1489[5] S.D. Grantham，L.H.王文，化工厂自动化故障排除的基本原理，计算机与化学工程：14（1990）783[6] P.K.李国忠，控制系统辨识，国立台湾大学机械工程研究所硕士论文（2004）。[7] A. Doucet，N. de Freitas，N. Gordon，Sequential Monte Carlo Methods in Practice，Springer，NewYork，2001。[8] R.B. Gopaluni，粒子滤波方法在缺失观测下的非线性过程识别，加拿大化学工程杂志：86（2008）1081[9] Jing Deng，Li Xie，Lei Chen，Shima Khatibisepehr，Biao Huang，Fangwei Xub，Aris Espejo.基于贝叶斯模型在线修正策略的软测量技术的开发与工业应用。Journal of Process Control：23（2013）317[10] J. Deng，B. Huang，具有丢失输出数据的非线性参数变化系统的识别，AIChE Journal（2012），http://dx.doi.org/10.1002/aic.13735。[11] 李文，等.粒子滤波器在非线性/非高斯贝叶斯跟踪中的应用[J].信号处理，2002，50（2）：174-188.[12]陈哲。王晓刚.粒子滤波器的应用[M].北京：清华大学出版社，2003.[13] X.绍湾Huang，J. Lee，F. Xu，中国春萤叶甲A. Espejo，约束贝叶斯状态估计-一种比较研究和一种新的基于粒子滤波的方法，过程控制杂志：20（2010）143[14] 放大图片作者：J.液相氧化法制取石油化工产品[J].烃类加工，1970a，49（3）：106.[15] 放大图片作者：J.液相氧化制石油化工（2）.动力学、传质与反应器设计[J].烃类工艺，1970b，49（11）：150.[16] 张春阳通过优化PTA氧化工艺参数，降低了氧化剂消耗，刘然然和李正明/ IERI Procedia 6（2014）6371乙酸[J].聚酯工业，1999，12（3）：39-42.[17] 云冥。PTA生产中降低PX消耗的分析[J].聚酯工业，2002，15（5）：35-37.[18] 吴德华，黄超，吴文革。优化TA工艺降低醋酸消耗[J].聚酯工业，2000，13（2）：37-40.