1 引言
间歇过程是一种重要的工业生产方式,被广泛应用于生产医药、食品、染料、香料及生化制品等小批量、高附加值产品[1].近年来,间歇过程的建模、监测、控制和优化已逐渐成为过程工业界关注的焦点[2].
间歇过程相对于连续过程具有显著的非线性、动态特性和时变特性,给间歇过程的建模带来了困难.间歇过程的机理复杂,建立精确的机理模型往往需要耗费大量的人力物力[3].对于间歇过程建模,很多传统方法只选取间歇过程的单个批次,将其当作一个连续过程来处理[4-6],无法很好地描述整个间歇过程动态.因此,Lu等人考虑间歇过程的多批次问题,提出基于批次-时间双维驱动模型[7-8].在此基础上,Zhao考虑间歇过程中各个批次之间的联系,研究了基于粒子滤波的间歇过程参数估计[9],但这些方法中假设模型结构已知,只是进行参数估计.同时,近年来,数据驱动建模的方法在间歇过程的建模中得到了大量的应用,例如为了有效利用偏最小二乘回归(PLS)、主成分分析(PCA) 多变量投影方法在处理高维、高度耦合数据时的独特优势,Nomikos和MacGrego提出多向主元分析(multiway PCA,MPCA) 和多向偏最小二乘回归(multiway PLS,MPLS)[10-12],将PCA和PLS推广应用到间歇过程.
但是,这些数据建模方法仍属于线性、静态建模方法,无法准确描述大多数间歇过程中强非线性、动态性和多阶段等特征.线性变参数(liner parameter varying,LPV) 建模方法是一种将多阶段、非线性的过程建模转化为线性多模型的方法,但传统LPV方法没有考虑间歇过程的多批次问题[13].
本文针对多批次、多阶段的间歇过程,综合上述研究,扩展LPV方法,提出多批次融合建模方法.首先根据间歇过程的多阶段动态特性,选取过程中的关键变量为调度变量[14],确定相应的工作点;然后依据调度变量与工作点将间歇过程的多批次数据直接划分为多个阶段;针对划分后的各阶段建立子模型,各阶段子模型均采用能较好逼近任意线性动态系统的ARX模型[15],并采用期望最大化算法(expectation maximization,EM) 对子模型参数进行估计,最后通过符合高斯分布的权重函数合成全局模型.采用该方法对青霉素发酵过程进行建模,结果表明基于多批次融合LPV模型具有更高的精度和更广泛的适用性.
2 多批次融合LPV建模 2.1 LPV模型LPV系统,是一种用多个线性系统逼近复杂非线性系统的方法.其模型系数矩阵是关于一个调度变量的函数,调度变量Zk(k=1,2,…,N) 是能反映过程中不同阶段特征的过程变量,其中N表示样本个数.多模型方法中通过选取合适的调度变量Z对过程进行划分、确定工况点后,可在各个工况点附近建立其相应的子模型.
选用一组带外加输入的自回归模型(AutoRegressive model with eXogenous input,ARX) 作为局部线性模型来逼近过程的动态,ARX模型具有较好的逼近线性动态系统的性能,表示如下:
(1) |
(2) |
式(1) 中yk∈Y⊆R代表输出,xk∈Rn+1代表系统在第k次采样时的回归量,ek∈R是均值为0、方差为σ2的高斯分布噪声,θZk∈Rn+1表示子模型的参数向量.式(2) 中na和nb分别表示输出和输入的阶数,uk∈U⊆Rt表示输入,na+t×nb=n.考虑到前一时刻采样数据与当前时刻数据相关性较大以及简略记,取一阶ARX模型,即na=nb=1,单输入单输出过程中t=1.
假设工况点个数记为M,即子模型的个数;每个工作点处的调度变量记为τm(m=1,2,…,M). τm本质上是非线性过程局部线性化的操作条件,在实际生产过程中,操作阶段明确,因此此处τm假设已知.
为使各子模型平滑融合,采用高斯分布函数为权重函数[16],k时刻第m个子模型的权重wkm表示如下:
(3) |
式中,Zk表示k时刻调度变量的值,om表示第m个模型的有效宽度.归一化后的权重函数可以表示为
(4) |
因此,对于输出的预测可以表示如式(5) 所示:
(5) |
ŷk表示k时刻全局模型的预测输出,
(6) |
间歇过程和连续过程有着明显的区别,具有极强非线性、生产反应周期短等特点.但传统LPV方法用于间歇过程中,仅选择单个批次过程,将间歇过程当作连续过程建模.此时一旦初始条件变化,则该模型难以适用于新的批次过程.例如对青霉素发酵过程多阶段融合建模中[13],选用青霉素发酵过程的一个批次,针对发酵过程的多阶段特点,采用基于EM算法的多模型方法建模.此方法中没有考虑间歇过程的多批次性,难以描述整个间歇过程的动态,同时当改变生产初始条件时,建立的模型适用性降低.下面将此问题展开详细描述.
本文针对青霉素发酵单批次过程采用LPV模型建模,仿真的实验数据来自于Pensim仿真平台.该平台的软件内核基于Birol模型,能够有效地模拟青霉素发酵过程[17]. Pensim仿真平台提供的部分变量和设定值的初始条件及默认值如表 1和表 2所示.在使用Pensim仿真平台时,设定相应的发酵时间和采样频率,在一定的条件下可得到青霉素发酵过程变量和产物浓度数据.
测量变量 | 范围 | 默认值 |
底物浓度/(g/L) | 5~50 | 15 |
青霉素浓度/(g/L) | 0 | 0 |
生物质浓度/(g/L) | 0~0.2 | 0.1 |
溶解氧浓度/(mmol/L) | 1~1.2 | 1.16 |
培养液体积/L | 100~200 | 100 |
CO2浓度/(mmol/L) | 0.5~1 | 0.5 |
pH值 | 4~6 | 5 |
反应器温度/K | 298~300 | 298 |
控制变量 | 设定值 | 范围 |
通风率/(L/h) | 8.6 | 3~10 |
搅拌功率/W | 30 | 20~50 |
底物流加速率/(L/h) | 0.042 | 0.035~0.045 |
底物流加温度/K | 296 | 208~296 |
pH值 | 5 | 5~6 |
反应器温度/K | 298 | 298~300 |
在仿真试验中,选取青霉素浓度为输出,由于青霉素关键过程变量--冷水流加速率--能较好地反应发酵过程的阶段特性,故将其选作模型的调度变量和模型的输入变量,其它变量均按照平台默认值设置.设定相应的发酵时间和采样频率,在平台默认的条件下得到一个批次中600组青霉素发酵过程变量和产物浓度数据,图 1表示青霉素发酵过程的输入-输出数据(输入也为调度变量).鉴于人为划分青霉素发酵过程的不同阶段依赖人工经验并且很难验证[18-19],本文对调度变量的采样值进行聚类,并将聚类中心作为主要工况点[13].青霉素发酵过程基本要经历菌体生长、青霉素合成、菌体自溶3个阶段,故采用k均值法将冷水流速率划分为3个阶段,聚类中心分别对应发酵过程的3个工况点:2.59 L/min,27.69 L/min,45.34 L/min.
对上述LPV模型进行求解,其仿真结果如图 2~图 4所示. 图 2表示的是采用单个批次建模的自我验证图.可以看出自我验证中输出的拟合度较高.为了进一步验证模型的有效性,通过设置底物浓度初始值分别为10 g/L、43 g/L,获得其它两个批次的数据进行交叉验证.验证结果如图 3所示. 图 4为图 3(b)的部分放大图.从图 3、图 4可以看出利用单个批次数据建立的模型交叉验证结果不够理想,图 3(a)中出现部分预测输出偏离实际输出较大,甚至在图 3(b)中出现了部分数据预测无效的情况,这是因为单个批次难以包含整个间歇过程的所有状态,在改变生产初始条件不同的生产过程下,该模型难以适用于所有的过程.因此本文提出多批次融合LPV建模方法对多批次、多阶段间歇过程进行建模.
2.3 多批次LPV模型假设用二维矩阵形式表示间歇过程多个批次数据,其中I表示操作批次,N表示单个批次内采样数据个数,如图 5所示.多批次融合LPV建模方法的主要思想是:首先根据间歇过程的多阶段动态特性,选取调度变量并确定相应的工作点;然后依据调度变量与工作点将间歇过程的多批次数据划分为多个阶段,运用多模型方法,针对每个阶段建立一个子模型,最后将局部模型融合成全局模型.
子模型如下:
(7) |
(8) |
式(7) 中yi,k∈Y⊆R代表第i个批次第k次采样的输出,xi,k⊆Rn+1代表第i个批次第k次采样时的回归量,ei,k∈R是均值为0、方差为σ2的高斯分布噪声,θZi,k∈Rn+1表示子模型的参数向量.式(8) 中na和nb分别表示输出和输入的阶数,ui,k∈U⊆Rt表示第i个批次第k次采样时的输入,na+t×nb=n.取一阶ARX模型,即na=nb=1,单输入单输出过程中t=1.
同LPV模型中一样,假设工况点个数记为M(即子模型的个数为M),每个工作点处的调度变量记为τm(m=1,2,…,M),假设τm已知且每个批次工况点相同.
为使各子模型平滑融合,采用高斯分布函数为权重函数[16],k时刻第m个子模型的权重wkm表示如下:
(9) |
式中,Zi,k表示第i个批次第k时刻调度变量的值,om表示第m个模型的有效宽度.归一化后的权重函数可以表示为
(10) |
因此,对于输出的预测可以表示如式(11) 所示:
(11) |
ŷk表示k时刻全局模型的预测输出,
(12) |
本文采用EM算法估计子模型参数. EM算法由Dempster等人[20]于1977年提出,旨在求解数据有缺失或隐含情况下参数的极大似然估计. EM算法由E (Expectation) 步和M (Maximization) 步组成. E步利用对隐藏变量的现有估计值和可观测数据,计算完整数据的对数似然函数的期望,记为Q函数;M步中,通过最大化Q函数校正参数向量. E步和M步不断迭代直至似然函数收敛.
引入隐藏变量Hk(k=1,2,…,N) 代表k时刻的子模型身份(该变量是离散型变量,取值为1,2,…,M中的一个).将可观测到的数据yi,k,ui,k,Zi,k(i=1,2,…,I,k=1,2,…,N) 表示为可观测数据集Cobs,其中yi,k,ui,k,Zi,k分别为青霉素发酵过程的输出、输入以及调度变量.隐藏变量Hk(k=1,2,…,N) 表示为不可观测数据集Cmis,过程中完整的数据集可以表示为C={Cobs,Cmis}={yi,k,ui,k,Zi,k,Hk}(i=1,2,…,I,k=1,2,…,N).
则EM算法可以归纳为以下步骤[19]:
第1步:初始化参数Θold.
第2步:E步,使用当前参数Θold估计近似的Q函数:
(13) |
第3步:M步,将Q函数最大化,获得新的迭代参数
(14) |
第4步:迭代过程,将获得新的参数代入第2步中,不断重复第2步和第3步,直至满足收敛条件为止.
下面详细介绍本文中EM算法的计算步骤.
E步:利用隐藏变量的现有估计值和可观测数据求完整数据集C期望的表达式如式(15) 所示.
(15) |
根据yi,k,ui,k,Zi,k,Hk之间的相关性,式(15) 可进一步化简为
(16) |
输入和调度变量与参数Θ之间不相关,则式(16) 中P(ui,N:1,Zi,N:1) 可视为常数,且对下述的计算不产生影响,可记作C1=P(ui,N:1,Zi,N:1). P(yi,k|xi,k,θm,σ2) 表示ARX模型分布概率,符合高斯分布,可以表示为
(17) |
上式中P(Hk=m|Zi,k,om) 表示权重函数,即P(Hk=m|Zi,k,om)=αkm,可由式(10) 得到. P(Hk=m|Θold,Cobs) 表示子模型身份后验分布,可由贝叶斯公式计算得到,记为
(18) |
将式(17)、式(18) 代入式(16) 中得到Q函数最终表达式如下:
(19) |
M步:通过Q函数极大化,求取新的参数.分别针对子模型参数θm和过程噪声方差σ2求偏导得到:
(20) |
(21) |
由于采用指数函数形式作为权函数来合成模型,无法得到局部模型有效宽度om的解析形式,可以将其转化为有约束情况下的极值求解问题[22].该最优值问题为
(22) |
式中的P(Hk=m|Θold,Cobs) 和P(Hk=m|Zi,k,om) 同前文描述一致.
根据式(20),(21) 和(22) 得到当前的参数、噪声方差和子模型有效宽度的估计,将其代入Q函数进行下一次迭代,直至收敛.
3 仿真结果为了解决间歇过程中传统LPV建模的不足,采用多批次LPV方法进行建模.通过设置底物浓度初值分别为8 g/L、30 g/L、45 g/L,获得3个批次的输入输出数据.同LPV模型中所述,在仿真试验中选取青霉素浓度作为输出,冷水流加速率作为模型的输入变量和调度变量. 3个批次输入输出数据如图 6所示.
为了方便下文验证模型的有效性,将青霉素3个批次数据融合成一组数据,如图 7所示.
对调度变量(输入) 的采样值进行聚类,并将聚类中心作为主要工况点[13].采用K均值法将其划分为3个阶段,其聚类中心分别对应发酵过程的3个工况点:8.29 L/min,46.3 L/min,76.33 L/min.应用多批次融合LPV建模方法,得到如图 8所示的青霉素浓度拟合效果.
为了进一步说明模型的可靠性,同LPV模型中一样,设定底物浓度初值分别为10 g/L、43 g/L,通过仿真平台产生其它两组数据对模型进行交叉验证,验证结果如图 9所示.可以从图 9中看出,预测输出与真实输出拟合程度较高,模型可靠性更高.从传统LPV建模交叉验证(图 3) 和多批次融合LPV建模的交叉验证(图 9) 中可以看出,对于非线性、多阶段性的间歇过程,本文提出的多批次LPV模型比传统LPV模型精度更高、适用性更加广泛,对于间歇过程建模具有一定的指示作用.
4 结论本文针对间歇过程采用多批次融合LPV方法进行建模.针对间歇过程多批次数据,依据调度变量将所有批次数据划分为不同阶段,基于EM算法建立基于不同工况点的子模型,并将归一化的指数函数作为权重函数,用来表示每个局部模型有效的概率,最后将所有的局部模型合成.针对青霉素发酵过程分别进行单个批次和多个批次融合建模的仿真,仿真结果表明多个批次融合建模的精度更高、适用性更广泛,验证了多批次融合LPV建模方法的有效性.
[1] | 杨志才. 化工生产中的间歇过程原理、工艺及设备[M]. 北京: 化学工业出版社, 2001. Yang Z C. The principle, process and equipment of batch process in chemical production[M]. Beijing: Chemical Industry Press, 2001. |
[2] | 钱路丰.贝叶斯核学习建模及在间歇过程中的应用研究[D].杭州:浙江大学, 2012. Qian L F. Research on Bayesian-Kernal learning modeling for batch processes[D]. Hangzhou:Zhejiang University, 2012. |
[3] | 贾立, 程大帅, 曹鲁明, 等. 基于数据的间歇过程时变神经模糊模型研究[J]. 计算机与应用化学, 2011, 28 (7): 915–918. Jia L, Cheng D S, Cao L M, et al. Research on data-based time-varying neuro-fuzzy model for batch processes[J]. Computer and Applied Chemistry, 2011, 28 (7): 915–918. |
[4] | Tao C, Morris J, Martin E. Particle filters for state and parameter estimation in batch processes[J]. Journal of Process Control, 2005, 15 (6): 665–673. DOI:10.1016/j.jprocont.2005.01.001 |
[5] | Leppävuori J T, Domach M M, Biegler L T, et al. Parameter estimation in batch bioreactor simulation using metabolic models:Sequential solution with direct sensitivities[J]. Industrial & Engineering Chemistry Research, 2011, 50 (21): 12080–12091. |
[6] | Wang G, Feng E, Xiu Z. Modeling and parameter identification of microbial bioconversion in fed-batch cultures[J]. Journal of Process Control, 2008, 18 (5): 458–464. DOI:10.1016/j.jprocont.2007.08.005 |
[7] | Lu N, Yuan Y, Gao F, et al. Two dimensional dynamic PCA for batch process monitoring[J]. Aiche Journal, 2005, 51 (12): 3300–3304. DOI:10.1002/aic.v51:12 |
[8] | Yuan Y, Gao F. Subspace identification for two-dimensional dynamic batch process statistical monitoring[J]. Chemical Engineering Science, 2008, 63 (13): 3411–3418. DOI:10.1016/j.ces.2008.04.007 |
[9] | Zhao Z, Huang B, Liu F. Parameter estimation in batch process using EM algorithm with particle filter[J]. Computers & Chemical Engineering, 2013, 57 (41): 159–172. |
[10] | Nomikos P, MacGregor J F. Monitoring batch processes using multiway principal component analysis[J]. American Institute of Chemical Engineers Journal, 1994, 40 (8): 1361–1375. DOI:10.1002/(ISSN)1547-5905 |
[11] | Nomikos P, Macgregor J F. Multi-way partial least squares in monitoring batch processes[J]. Chemometrics & Intelligent Laboratory Systems, 1995, 30 (1): 97–108. |
[12] | Nomikos P, MacGregor J F. Multivariate SPC charts for monitoring batch processes[J]. Technometrics, 1995, 37 (1): 41–59. DOI:10.1080/00401706.1995.10485888 |
[13] | 熊伟丽, 姚乐, 徐保国. 基于EM算法的青霉素发酵过程多阶段融合建模[J]. 化工学报, 2014, 65 (12): 4935–4941. Xiong W L, Yao L, Xu B G. Multi-stage fusion modeling for penicillin fermentation process based on EM algorithm[J]. CIESC Journal, 2014, 65 (12): 4935–4941. |
[14] | Zhu Y C, Xu Z. A method of LPV model identification for control[C]//Preprint of IFAC World Congress. Kington, UK:IFAC, 2008:5018-5023. |
[15] | Ljung L. System identification:Theory for the user[M]. Boston, NJ, USA: Birkhäuser, 1987: 163-173. |
[16] | Huang J, Ji G, Zhu Y, et al. Identification of multi-model LPV models with two scheduling variables[J]. Journal of Process Control, 2012, 22 (7): 1198–1208. DOI:10.1016/j.jprocont.2012.05.006 |
[17] | Liu Y. Modelling of the penicillin fermentation process via LS-SVM based on PENSIM simulator[J]. Chemical Reaction Engineering & Technology, 2006, 22 (3): 252–258. |
[18] | Wei S, Yan M, Palazoglu A, et al. A method for multiphase batch process monitoring based on auto phase identification[J]. Journal of Process Control, 2011, 21 (4): 627–638. DOI:10.1016/j.jprocont.2010.12.003 |
[19] | Yuan Y, Gao F. Phase and transition based batch process modeling and online monitoring[J]. Journal of Process Control, 2009, 19 (5): 816–826. DOI:10.1016/j.jprocont.2008.11.001 |
[20] | Dempster A P, Rubin D B. Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society, 1977, 39 (1): 1–38. |
[21] | Mclachlan G, Krishnan T. The EM algorithm and extensions (Wiley series in probability and statistics)[J]. Journal of Classification, 2007, 15 (1): 154–156. |
[22] | Rasmussen C E. Gaussian processes in machine learning[M]. Berlin, Germany: Springer, 2004: 63-71. |