含时滞测量值下间歇过程的双维状态估计

祁鹏程，赵忠盖，刘飞



含时滞测量值下间歇过程的双维状态估计

祁鹏程，赵忠盖，刘飞

（江南大学轻工过程先进控制教育部重点实验室，江苏 无锡 214122）

基于粒子滤波研究了间歇过程的状态估计问题。根据间歇过程双维动态特性,针对关键参数在线检测精度低、离线分析时滞大等问题，分别建立一种双维状态转移模型和时滞测量模型，并利用贝叶斯方法及前/后向平滑，提出一种含时滞测量值下的双维状态估计算法。该算法通过融合先前批次和时滞测量值的信息提高估计精度，并且克服了离线采样周期和时滞时间不确定的问题。在数字仿真和啤酒发酵过程中的仿真应用验证了该算法的有效性。

间歇过程；双维状态空间模型；时滞；粒子滤波；贝叶斯方法；前/后向平滑

引 言

间歇过程也叫批次过程，广泛应用于化工、食品和医药等行业[1]。生产过程中，关键参数的变化，如发酵过程中残糖浓度、菌体浓度等，直接影响产品的质量[2]。但是受限于传感检测技术，通常很难获得稳定准确的关键参数在线检测值，而是通过估计的方法获取关键参数的估计值。在贝叶斯框架下，给定检测值，通过近似计算状态空间的后验概率密度函数可以解决间歇过程中关键参数的估计问题[3]。

目前，许多学者对间歇过程中的状态估计问题进行了大量研究[4-10]。考虑间歇精馏塔的非线性特性，以及在线成分分析成本高等问题，文献[5]基于扩展卡尔曼滤波（EKF）设计最优估计器从而得到系统的组成成分。然而，卡尔曼及其扩展形式只能处理高斯噪声的滤波问题，基于此，文献[7]中引入粒子滤波(PF)算法并应用到间歇过程的状态估计中，取得了良好效果。上述研究只考虑了间歇过程的非线性、非高斯特性，状态仅随时间转化，忽略了批次之间的动态特性。实际上，典型的间歇过程包含许多批次，关键参数的动态变化不仅存在于时间维度，还存在于批次维度[8]。与连续过程不同，间歇过程中所有原料分批次处理，当前批次会根据先前批次的信息来调整操作步骤，如产品的质量。文献[9]首次考虑了当前批次的初始状态与先前批次相关的情况，利用平滑的方式将先前批次的信息融合到当前状态的估计中，提高了估计精度。文献[10]在上述文献的基础上考虑了所有状态，当前状态既与当前批次的先前状态相关也与先前批次相同采样时间的状态相关，并建立双维状态空间模型对关键参数进行估计，取得了良好效果。

另一方面，离线检测在间歇过程中普遍存在，离线检测一般由人工在实验室完成且检测的样本也是由人工采集。上述过程所用的时间相比于在线采样时间间隔不可忽略，即离线检测得到的测量值有很大的时滞。但在实验室分析得到的数据比在线分析仪得到的更准确，因此融合时滞测量值的信息可以进一步提高关键参数的估计效果[11]。文献[11]针对线性系统通过卡尔曼滤波和信息融合技术，实现两种检测值的信息融合。鉴于卡尔曼滤波简单、计算快的优点，文献[12]基于EKF分别分析了全增广和样本增广在处理时滞测量值的性能。文献[13]根据时滞测量值与状态为积分关系的特性，建立新模型并提出变维无迹卡尔曼滤波解决此模型的状态估计问题。但是，卡尔曼滤波及其扩展形式在使用时噪声需满足高斯分布，有很大的局限性。

本文研究间歇过程中关键参数的估计问题,针对关键参数在线检测精度低、离线分析时滞大等问题，结合间歇过程的双维动态特性，建立新的状态空间模型，基于PF提出一种含时滞测量值下的双维状态估计算法。该算法不仅考虑了先前批次的信息，还融合了离线检测的时滞测量值。分别针对仅有在线检测值和两种检测值并存两种情况，利用前/后向平滑(FBS)对先前批次进行平滑从而融合先前批次的信息。处理时滞测量值时，运用贝叶斯方法更新粒子权重实现两种检测值的信息融合，通过调整平滑步数解决时滞时间不确定的情况，并进行了数字仿真和啤酒发酵过程的仿真应用。

1 问题描述

1.1 间歇过程的双维动态特性

一个典型的间歇过程包含许多批次。在单个批次中，与连续过程一样，状态随时间演化，具有时间动态特性。而在批次间，由于原料变化、运行调整等原因，状态在批次维度上同样具有动态特性。考虑到批次之间开始时加入的原料有缓慢的变化，且一般会根据最后产品的质量来调节下一批次的原料，因此考虑当前批次的初始状态与先前批次的初态和末态相关的情况。对于其他状态，同样会根据先前批次相同采样时间的运行状态调整操作，所以考虑其他状态既与当前批次的先前状态相关也与先前批次相同采样时间的状态相关的情况[10]。

假设间歇过程包含个批次，每个批次有个采样时刻，则双维动态模型示意图如图1所示。其中，和代表第批次在采样时刻的状态和检测值。水平和竖直的箭头代表状态的转换方向，斜箭头表示每个采样时刻得到的检测值。实线表示两个状态之间没有其他的状态，而虚线代表有其他的状态。由于第一批次没有先前批次，与连续过程一样，只有时间上的动态特性。而其他批次不仅具有时间维度的动态关系，还有批次维度的动态关系。

图1 间歇过程双维动态示意图

1.2 间歇过程的时滞特性

间歇过程中，在线检测值直接由传感器测得，即时但精度低，离线检测值通过人工离线分析得到，有时滞但精度高，且两种检测值的采样周期不同，给间歇过程中关键参数的估计带来了困难，如何融合时滞测量值是提高估计效果的关键。间歇过程含时滞测量值的示意图如图2所示，其中，表示离线检测时刻，表示第批次在时刻的在线检测值，表示第批次在时刻的离线检测值，表示离线检测值的时滞时间。垂直的箭头表示在线检测值没有时滞。虚线箭头表示离线检测值在时刻采样，检测值经过时滞在时刻得到。因为离线采样周期较大，通常最近一次离线采样的检测值在进行下一次离线采样之前就已经得到。

图2 间歇过程含时滞测量值示意图

另外，在实际生产过程中，在线采样频率固定且周期小，而离线采样由人工完成，采样频率不固定且周期大，且由于离线分析的操作时间不是确定的，因此时滞时间也可能是时变的。

定义1：第批次时刻所有在线检测值的集合为，离线检测值的集合为。

2 间歇过程状态估计的实现

首先根据间歇过程的双维动态特性和时滞特性，分别建立新的双维状态转移模型和时滞测量模型。考虑两种检测值的采样周期不同，且离线采样周期不固定，本文从两方面进行分析：仅有在线检测值的时刻和在线检测值、离线检测值同时存在的时刻。

2.1 模型描述

2.1.1 双维状态转移模型

第2及以后批次需考虑先前批次的信息。对于初始状态，即＞1，，与上一批次的初态和末态相关，其模型如下

对于其他状态，即＞1，＞1，与当前批次的先前状态和先前批次相同采样时间的状态相关，模型 如下

2.1.2 时滞测量模型 根据图2可以得到间歇过程含时滞测量值的测量模型如下

2.2 状态估计：仅有在线检测值的时刻

由模型描述可知，第1批次状态仅随时间变化，属于单维情况，而其他批次属于双维情况，所以分单维和双维两种情况分析。

2.2.1 单维估计

根据第1批次的状态空间模型，结合贝叶斯滤波，估计步骤如下[14]：

预测

更新

对于非线性非高斯过程，PF利用Monte Carlo（MC）方法估计目标状态后验分布，取得了良好效果。其原理是用样本均值代替积分运算[15]

其中，为粒子数。但是真实的状态分布情况未知，通过重要性采样分布可以获取先验粒子[16]，则

考虑粒子退化的问题，需要对粒子进行重采样。重采样后，每个粒子的权重变为1/[18]。

2.2.2 双维估计

由于其他批次的初始状态与其他状态的转移方程不同，本节又分初始状态和其他状态进行考虑。根据状态空间模型和PF算法可以得到初始状态的预测步骤如下

其中，假设初始状态与末态是相互独立的，即

因为随着得到的测量值信息增多，初态对末态的影响较小，且每个批次的长度较长，可以忽略初态与末态的关系[10]。而在前一批次估计中已得到，只要计算可得，很明显是一个平滑分布。本文采用FBS平滑算法，该方法保存每个采样时刻的粒子作为平滑粒子，在当前时刻估计状态的基础上递推计算先前时刻的平滑状态，且有效地解决了平滑过程中粒子退化的问题[19]。具体如下

（2）当＞2,平滑方式如下

初始状态的更新步骤如下

对于其他状态，根据状态空间模型和PF算法可以得到其他状态的估计步骤如下：

预测

与式（11）的处理方式一样，同样需要假设

表示两个状态在时间维度上相差一个批次的长度就相互独立。通过计算平滑分布即可得到。而在计算前必须先计算，具体如下

那么可得

更新

2.3 状态估计：两种检测值并存的时刻

在处理时滞测量值前，首先根据在线检测值的信息更新当前时刻的粒子权重，然后再通过贝叶斯方法融合含时滞的离线检测值，再次更新粒子权重。具体如下

则

对于单维估计，平滑的处理方式与式（15）相同，对于双维估计与式（17）相同，只是在时，是和同时作为时刻检测值时的权重，通过递推计算即可得到分布对应的粒子和权重。

实际操作中，由于离线检测值是需要经过人工实验分析才能得到，操作时间具有不确定性，即时滞时间是不确定的。另外，离线采样周期不固定，即离线检测时刻不固定。从本节可以看出融合时滞测量值方法的关键一步是计算，在本文式是采用FBS算法进行平滑处理。因为当前时刻（时滞值到达时刻）等于离线检测时刻与时滞时间之和，在知道和下就可以确定，而和在实际过程中很容易得到，所以一旦确定和就可以根据当前调整平滑步数计算出。因此本文提出的算法可以克服时滞时间和采样周期不确定的问题。

综上，含时滞测量值下间歇过程的双维状态估计示意图如图3所示，具体步骤如下。

（1）第1批次的状态估计实现

⑤重采样：根据归一化后的权重得到个近似服从分布的粒子，并重置所有粒子权重为1/。

（2）其他批次的状态估计实现

⑦重采样：根据归一化后的权重得到个近似服从分布的粒子，并重置所有粒子权重为1/。

3 算法验证

为了比较间歇过程中考虑时滞测量值与不考虑时滞测量值的估计效果，选择每个批次的初始状态误差（ISE）和均方根误差（RMSE）[21]作为性能指标

3.1 数字仿真

根据文献[8]得到间歇过程的数字仿真模型如下。

（1）双维状态转移模型

第1批次

初始状态

除初始状态之外的其他批次

（2）时滞测量模型

图4为数字仿真的效果图。其中，黑色虚线是没有考虑时滞测量值的误差，红色点线代表融合了时滞测量值的误差。由图可知，随着批次的增加，状态的ISE和RMSE都越来越小，趋于收敛，且在考虑时滞测量值的情况下状态的ISE和RMSE比不考虑时滞测量值时的小，表明在考虑时滞测量值的情况下关键参数的估计效果较好。

图4 数字仿真无时滞值与含时滞值状态的初始状态误差和均方根误差

3.2 啤酒发酵过程

发酵是啤酒生产中耗时最长的工序，是一个典型的间歇生产过程，生产过程中两个关键参数为残糖浓度和菌体浓度。在发酵阶段，残糖浓度不可过高也不可过低，过高则发酵迟缓、啤酒发酵度低、生物稳定性差，过低则微生物耗氧能力下降，溶氧偏高。而菌体浓度高低直接影响啤酒的风味，浓度过高易致酵母衰退、自溶、酒质下降，浓度过低则发酵缓慢、发酵副产物增加。啤酒发酵过程关于残糖浓度和菌体浓度的状态转移模型如下[22]

在实际生产过程中，啤酒发酵是按批量处理的，且在发酵过程中通常会离线采样分析来得到残糖浓度和菌体浓度更准确的测量值，对生产过程进行监控、控制和优化。但离线分析有时滞，因此如何融合先前批次和时滞测量值的信息是提高残糖浓度和菌体浓度估计效果的关键。根据图1和图2建立啤酒发酵过程含时滞测量值的双维动态空间模型。

（1）双维状态转移模型

对于啤酒发酵的第1批次，没有先前批次的信息，当作连续过程处理

对于初始状态，在上一批次的发酵结束后，会根据啤酒成品的质量来调整本批次原料的糖度和酵母浓度

虽然当前批次的发酵会根据原料变化、过程飘移等问题调整发酵的相应操作，但总体上还是重复上一批次的整个过程。因此当前批次的状态与当前批次的先前状态和先前批次相同采样时间的状态相关

（2）时滞测量模型

生产过程通常配有在线分析仪，如液相色谱分析仪、YSI2700等，从而得到残糖浓度和菌体浓度的测量值对生产过程进行监控，但误差较大。因此为了提高啤酒成品的质量，大概每隔0.5 h对生产过程离线检测一次，通过实验室分析得到更准确的测量值。模型如下

图5为啤酒发酵过程估计效果图，分别为残糖浓度和菌体浓度的ISE和RMSE。黑色虚线是没有考虑时滞测量值的误差，红色点线代表融合了时滞测量值的误差。从图中可以看出，随着批次的增加，状态的ISE和RMSE都越来越小，趋于收敛。由于每个批次只离线检测两次，即只能得到两个时滞测量值的信息，导致估计效果没有明显提高，但在考虑时滞测量值的情况下状态估计的总体效果要好于不考虑时滞测量值时的情况。

图5 啤酒发酵过程无时滞值与含时滞值状态的初始状态误差和均方根误差

4 结 论

本文针对间歇过程，建立一种双维状态空间模型，研究了含时滞测量值下关键参数的估计问题。分别分析了仅有在线检测值和两种检测值并存两种情况。考虑到间歇过程非线性、非高斯特性，引入PF算法进行状态估计。为充分利用离线检测值，以及间歇过程沿批次维度上历史批次信息，通过贝叶斯方法对粒子权重进行更新，并采用FBS算法进行平滑，实现了多源信息的融合。数字仿真和啤酒发酵的应用仿真结果说明：在间歇过程中融合时滞测量值可以有效地提高关键参数的估计效果。

References

[1] 侯迪波. 间歇生产过程中的知识发现方法研究 [D]. 杭州: 浙江大学, 2004.
HOU D B. Research on knowledge discovery method in batch process [D]. Hangzhou: Zhejiang University, 2004.

[2] 李慧娟. 青霉素发酵过程建模与优化 [D]. 沈阳: 东北大学, 2010.
LI H J. Modeling and optimization of penicillin fermentation process [D]. Shenyang: Northeastern University, 2010.

[3] 姜蕾蕾. 啤酒发酵过程状态监测及控制策略的研究 [D]. 哈尔滨: 哈尔滨理工大学, 2014.
JIANG L L. Study on state monitoring and control strategy of beer fermentation process [D]. Harbin: Harbin University of Science and Technology, 2014.

[4] KOZUB D J, MACGREGOR J F. State estimation for semibatch polymerization reactors [J]. Chemical Engineering Science, 1992, 47 (5): 1047-1062.

[5] OISIOVICI R M, CRUZ S L. State estimation of batch distillation columns using an extended Kalman filter [J]. Chemical Engineering Science, 2000, 55 (20): 4667-4680.

[6] WANG J L, ZHAO L Q, YU T. On-line estimation in fed-batch fermentation process using state space model and unscented Kalman filter [J]. Chinese Journal of Chemical Engineering, 2010, 18 (2): 258-264.

[7] CHEN T, JULIAN M, ELAINE M. Particle filters for state and parameter estimation in batch processes [J].Journal of Process Control, 2005, 15 (6): 665-673.

[8] YAO Y, GAO F. Batch process monitoring in score space of two-dimensional dynamic principal component analysis (PCA) [J]. Industrial & Engineering Chemistry Research, 2007, 46 (24): 8033-8043.

[9] ZHAO Z G, HUANG B, LIU F. Bayesian method for state estimation of batch process with missing data [J]. Computers and Chemical Engineering, 2013, 53 (11): 14-24.

[10] ZHAO Z G, HUANG B, LIU F. State estimation in batch process based on two-dimensional state-space model [J]. Industrial & Engineering Chemistry Research, 2014, 53 (50): 19573-19582.

[11] GUO Y F, ZHAO Y, HUANG B. Development of soft sensor by incorporating the delayed infrequent and irregular measurements [J]. Journal of Process Control, 2014, 24 (11):1733-1739.

[12] GOPALAKRISHNAN A, KAISARE N S, NARASIMHAN S. Incorporating delayed and infrequent measurements in extended Kalman filter based nonlinear state estimation [J]. Journal of Process Control, 2011, 21 (1): 119-129.

[13] GUO Y F, HUANG B. State estimation incorporating infrequent, delayed and integral measurements [J]. Automatica, 2015, 58: 32-38.

[14] FERDINAND E U. State-space estimation with a Bayesian filter in a coupled PDE system for transient gas flows [J]. Original Research Article, 2015, 39 (2): 682-692.

[15] ANDRIEU C, DOUCET A, SINGH S,. Particle methods for change detection, system identification, and control [J]. Proceedings of the IEEE, 2004, 92 (3): 423-443.

[16] DOUCET A, GODSILL S, ANDRIEU C. On sequential Monte Carlo sampling methods for Bayesian filtering [J]. Statistics and Computing, 2000, 10 (3): 197-208.

[17] STORVIK G. Particle filters for state-space models with the presence of unknown static parameters [J]. IEEE Trans. Signal Processing, 2002, 50 (2): 281-289.

[18] RIGATOS G G. particle filtering for state estimation in nonlinear industrial systems [J]. Instrumentation and Measurement, 2009, 58 (11): 3885-3900.

[19] NURMINEN H, RISTIMAKI A, ALI-LOYTTY S,. Particle filter and smoother for indoor localization[C]//4th International Conference on Indoor Positioning and Indoor Navigation. Montbéliard, France: IEEE, 2013: 1-10.

[20] ORGUNER U, GUSTAFSSON F. Storage efficient particle filters for the out of sequence measurement problem[C]//11th International Conference on Information Fusion. Piscataway, NJ, USA: IEEE, 2008: 1-8.

[21] KOTTAKKI K K, BHARTIYA S, BHUSHAN M. State estimation of nonlinear dynamical systems using nonlinear update based unscented Gaussian sum filter [J]. Journal of Process Control, 2014, 24 (9): 1425-1443.

[22] ZHAO Z G, HUANG B, LIU F. Parameter estimation in batch process using EM algorithm with particle filter [J]. Computers and Chemical Engineering, 2013, 57 (15): 156-172.

Two dimensional state estimation in batch process with delayed measurements

QI Pengcheng, ZHAO Zhonggai, LIU Fei

(Key Laboratory of Advanced Control for Light Industry Processes, Ministry of Education, Jiangnan University,
Wuxi 214122, Jiangsu, China)

This work investigates the state estimation in a batch process based on the particle filter method. Considering that the two dimensional dynamics and the key parameters are obtained online with low accuracy or offline with large time delay, a two dimensional state transition model and delayed measurement model are developed. In addition, a two dimensional state estimation algorithm is proposed by using Bayesian method and the forward-backward smoothing algorithm. The proposed algorithm improves the estimation accuracy by fusing information of previous batches and delayed measurements, and overcomes the influence of the uncertainty of offline sampling period and time delay. The applications in a numerical example and a beer fermentation show the effectiveness of the proposed method.

batch process; two dimensional state space model; time delay; particle filter; Bayesian method; forward-backward smoothing

supported by the National Natural Science Foundation of China (61134007, 61573169).

date: 2015-12-14.

Prof. ZHAO Zhonggai, gaizihao@ jiangnan.edu.cn

TQ 920.6

A

0438—1157（2016）09—3784—09

10.11949/j.issn.0438-1157.20151894

国家自然科学基金项目（61134007，61573169）。

2015-12-14收到初稿，2016-04-28收到修改稿。

联系人：赵忠盖。第一作者：祁鹏程（1991—），男，硕士研究生。