引言

间歇过程因其适应多品种、多规格和高质量的市场要求，近年来被广泛应用于精细化工、生物制药、食品、聚合物反应、金属加工等领域^[1]。间歇过程的一个固有特征是具有多时段性，即在不同时段呈现出不同的过程相关特性。若忽略间歇过程的多时段过程特性，将一次间歇操作的多时段过程看作一个整体进行间歇过程建模，会导致所建立的过程模型难以准确描述间歇过程，直接影响生产过程的安全运行和产品质量，造成生产企业经济损失甚至人员伤亡。因此，进行多时段间歇过程的时段划分，对于提高间歇过程模型精度及产品质量，保障生产过程中的人身、财产安全都具有重要的意义^[2-3]。

目前应用较广泛的间歇过程时段划分方法主要有基于聚类分析^[4-5]和基于模型识别^[6-7]的方法。基于聚类分析的方法通过无监督机器学习，将具有相似特征的过程数据归为一类，从而实现时段划分^[8-9]。该方法原理简单、计算复杂度较低，但其准确性严重依赖于时段个数的设定和时段初始中心的选取，影响了时段划分结果的精度。基于模型识别的间歇过程时段划分方法不需要设定时段个数和初始时段中心，只需分析间歇过程运行时间顺序的过程特性并建立统计模型，通过评估时段划分对统计模型特性的影响确定合适的时段划分点^[10]，但是该方法在处理间歇过程数据时均需将三维过程数据展开为二维数据，从而损失了间歇过程的三维数据特性，直接影响多时段间歇过程的时段划分结果的准确性。

平行因子分解2(PARAFAC2)^[11]是主成分分析的一种高维推广，因其能够在不对三维数据作二维展开的情况下直接进行主成分分析，可以有效地保持数据内部结构的完整性，在数据挖掘、人脸识别、特征提取等领域得到广泛应用^[12]。Louwerse等^[13]首次将平行因子法应用于间歇过程的三维过程数据建模，提升了过程模型精度；Matero等^[14]利用PARAFAC2解决了三维不等长批次过程数据的建模；考虑间歇过程的多时段特性，Luo等^[15]利用聚类分析方法划分时段并基于PARAFAC2建立各个时段的过程模型，相比整体PARAFAC2模型，多时段PARAFAC2模型具有更高的准确性。以上研究虽然获得了较高精度的过程模型，但均采用基于聚类分析或模型识别方法划分间歇过程的多时段，没有充分利用PARAFAC2处理三维过程数据所具有的优势划分多时段，仅是利用PARAFAC2建立各个时段的过程模型，多时段划分结果的准确性仍然制约着间歇过程模型精度的提升。

针对多时段间歇过程的时段划分问题，本文利用间歇过程三维数据建立三维时变数据集，首先对每个时间片矩阵进行PARAFAC2建模，然后从初始时刻开始，将时间片矩阵按照时间顺序依次添加三维的时间块矩阵，并采用PARAFAC2建立三维过程数据模型；再分别计算各模型的平方预测误差(square prediction error, SPE)统计量，通过衡量新加入的时间片对时间块模型统计量的影响来判断是否具有相似的过程特性，从而获取时段划分点；然后通过利用时段评价划分指标(PPCI)得到最佳的时段划分结果。最后利用青霉素发酵过程仿真实验验证了所提方法的有效性。

1 平行因子分解2(PARAFAC2)

给定一个三维过程数据集X={X₁, X₂, …, X_I}，X_i∈R^K_i×J，i=1, 2, …, I，I为批次个数，J为变量个数，K_i为每批次的采样点个数。对矩阵X进行PARAFAC2建模可得

$ {\mathit{\boldsymbol{X}}_i} = {\mathit{\boldsymbol{F}}_i}{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{E}}_i} $

(1)

式中，F_i(K_i×R)为得分矩阵(R为主元个数)，D_i(R×R)为加权对角矩阵，A(J×R)为负载矩阵，E_i(K_i×J)为残差矩阵。三维数据矩阵X的PARAFAC2模型如图 1所示。

为了使该模型具有唯一性，需增加约束条件。设F_i^TF_i为恒定的矩阵，即

$ \mathit{\boldsymbol{F}}_a^{\rm{T}}{\mathit{\boldsymbol{F}}_a} = \mathit{\boldsymbol{F}}_b^{\rm{T}}{\mathit{\boldsymbol{F}}_b} = \mathit{\Phi },\forall a,b = 1,2, \cdots ,I $

(2)

令F_i=P_iF，其中P_i(K_i×R)为正交矩阵，F(R×R)为定值矩阵。由于

$ \mathit{\boldsymbol{F}}_i^{\rm{T}}{\mathit{\boldsymbol{F}}_i} = {\mathit{\boldsymbol{F}}^{\rm{T}}}\mathit{\boldsymbol{P}}_i^{\rm{T}}{\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{F}} = {\mathit{\boldsymbol{F}}^{\rm{T}}}\mathit{\boldsymbol{F}},i = 1,2, \cdots ,I $

(3)

因此PARAFAC2模型可表示为

$ {\mathit{\boldsymbol{X}}_i} = {\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{E}}_i} = {\mathit{\boldsymbol{T}}_i}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{E}}_i} $

(4)

式(4)中T_i=P_iFD_i为组合得分矩阵。

采用直接拟合法^[16]可以得到P_i，F，D_i和A^T，其计算式为

$ \sigma \left( {\mathit{\boldsymbol{F}},{\mathit{\boldsymbol{D}}_1},{\mathit{\boldsymbol{D}}_2}, \cdots ,{\mathit{\boldsymbol{D}}_I},{\mathit{\boldsymbol{P}}_1},{\mathit{\boldsymbol{P}}_2}, \cdots ,{\mathit{\boldsymbol{P}}_I},\mathit{\boldsymbol{A}}} \right) = \sum\limits_{i = 1}^I {{{\left\| {{\mathit{\boldsymbol{X}}_i} - {\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{A}}^{\rm{T}}}} \right\|}^2}} $

(5)

式中，P_i^TP_i=I，D_i为对角矩阵。

2 基于PARAFAC2的多时段间歇过程时段划分 2.1 三维数据标准化过程

实际间歇过程中采集的数据均为间歇过程中可被测量的状态参数，为消除变量间的量纲关系，需要对间歇过程三维数据X={X₁, X₂, …, X_I}进行标准化。首先将间歇过程三维数据X沿采样点时刻方向展开为K个二维数据矩阵$ \mathit{\boldsymbol{\widetilde X}}$(I×J); 然后再对每个$ \mathit{\boldsymbol{\widetilde X}}$(I×J)进行中心化和归一化处理，得到标准化后的二维数据矩阵$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over X} }}}$(I×J)；最后为保证间歇过程三维数据的完整性，还需将K个二维数据矩阵$ {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over X} }}}$(I×J)按照采样点时刻方向依次叠加，恢复原始的三维结构，X={X₁, X₂, …, X_I}为标准化的间歇过程三维数据。标准化过程的具体流程如图 2所示。

2.2 时段划分算法

对标准化预处理后的间歇过程数据X进行时段划分，获取基于PARAFAC2的时间片和时间块矩阵模型，并确定时段划分点。

2.2.1 时间片建模

设一个间歇过程具有K个采样点和J个测量变量，则每个测量批次可获得一个K×J的数据矩阵，重复I批次的测量步骤后，得到的数据可以表述为一个三维数据矩阵X(I×J×K)；提取一个批次内的各个采样点上的变量得到二维时间片矩阵X_k(I×J)。

对每一个时间片矩阵X_k(I×J)标准化后得到的结果X_k进行PARAFAC2分解，建立时间片PARAFAC2模型，其分解过程可表示为

$ {{\mathit{\boldsymbol{\bar X}}}_k} = {\mathit{\boldsymbol{F}}_k}{\mathit{\boldsymbol{D}}_k}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {{\mathit{\boldsymbol{\bar E}}}_k} = {{\mathit{\boldsymbol{\bar T}}}_k}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {{\mathit{\boldsymbol{\bar E}}}_k} = \sum\limits_{r = 1}^R {{{\mathit{\boldsymbol{\bar T}}}_{k,r}}\mathit{\boldsymbol{A}}_r^{\rm{T}}} + {{\mathit{\boldsymbol{\bar E}}}_k} $

(6)

式中，r表示不同的分解方向。设每个时间片的主元个数为R_k，选择出现次数最多的R_k值作为所有时间片PARAFAC2模型的主元个数R。

2.2.2 时间块建模

从间歇过程初始点开始，依次将每一个时间片按照k的方向进行叠加(k代表当前过程时间)，得到三维的时间块矩阵X_k(I×J×K)，其构建过程如图 3所示。

类似地，对每一个时间块矩阵X_k进行PARAFAC2分解，最终建立时间块PARAFAC2模型。

2.2.3 基于模型控制限变化的时段划分

计算每一个时间片矩阵X_k、时间块矩阵X_k对应于各个批次的SPE统计量，其计算式为

$ {\mathit{\boldsymbol{S}}_{k,i}} = {\mathit{\boldsymbol{e}}_{k,i}}\mathit{\boldsymbol{e}}_{k,i}^{\rm{T}} $

(7)

式中，i是批次序号，k是运行时间，e_{k, i}是残差矩阵E_k对应于k时刻第i批次的列向量。通过核密度估计方法^[17]分别得出X_k、X_k在k时刻对应的SPE控制限CLS_k、CLB_k。

引入动态调节因子α来反映时间片模型与时间块模型的精度折中程度。时段划分的步骤为：将时间片控制限CLS_k与从初始时刻到k时刻的时间块控制限CLB_k进行比较，如果存在连续3个数据满足CLB_k>αCLS_k，说明新加入的时间片对该时间块模型的重大影响导致当前时间块模型无法精确表征该时间块内的所有时间片数据，由此判断当前k时刻之前的时间块隶属于同一个时段；然后根据所获得的k时刻指示信息将该时段数据全部移除，再将余下的间歇过程运行时间作为新的数据起始点，重复进行模型控制限的判断，直到所有的时段划分点被确定。

2.2.4 分段效果评价

通过调节α的值可以获取不同时段个数的间歇过程时段划分结果。采用时段划分评价指标(PPCI)^[9]来确定最佳时段个数。首先，对α设定一定的范围和步长，利用本文方法依次对间歇过程进行时段划分，得到一组对应不同α值的时段划分结果；其次，分别计算每个α对应结果的PPCI指标，PPCI指标最小值对应的时段划分结果即为最佳时段分段结果。

基于PARAFAC2的多时段间歇过程时段划分算法流程如图 4所示。

3 实验验证及结果分析

为了验证本文方法的有效性，基于青霉素发酵过程进行了间歇过程的多时段划分实验，并与K均值、FCM、WKM、SCFCM 4种间歇过程时段划分方法进行结果比较。

3.1 实验数据

青霉素是一种临床应用广泛的抗生素，其生产过程是典型的多时段间歇过程。利用Pensim仿真平台^[18]为青霉素发酵过程测量数据时段划分提供标准的过程仿真数据，选取如表 1所示的6个时段敏感变量，生成20个正常批次数据，采样时间400 h，采样间隔1 h，过程数据为X(20×6×400)。

3.2 实验过程

设置SPE统计量控制限的置信度为95%，采用基于PARAFAC2模型进行时段划分，通过改变动态调节因子α的取值来获取不同的时段划分结果。当α的取值范围为[2.0, 7.7]，步长为0.1时，对应的时段划分结果如表 2所示。

计算不同α对应的PPCI值如图 5所示，得出α取值为3.8时，时段划分结果的PPCI值最小。α=3.8的划分结果如图 6所示，可以看出时段划分点分别为p₁=28，p₂=59，p₃=216，p₄=369，划分出1~28、29~59、60~216、217~369和370~400共5个时段。

3.3 方法比较

本文选取Davies-Bouldin(DB)指标和Compactness(CP)指标^[19]作为比较各种划分方法性能的评价指标。DB指标计算式为

$ {x_{{\rm{DB}}}} = \frac{1}{L}\sum\limits_{i = 1}^L {\mathop {\max }\limits_{j \ne i} \left( {\frac{{{b_i} + {b_j}}}{{\left\| {{v_i} - {v_j}} \right\|}}} \right)} $

(8)

式中，L是时段个数，ν_i和ν_j分别是时段i和时段j的中心，b_i和b_j分别为时段i、j的平均误差，计算式为

$ {b_i} = \frac{1}{{{N_i}}}\sum\limits_{{x_j} \in {\mathit{\Omega }_i}} {{{\left\| {{x_j} - {v_i}} \right\|}^2}} $

(9)

式中，Ω_i、N_i分别为时段i包含的样本集合和样本个数。DB指标的值越小，代表时段内部各个样本越紧凑，各时段间的距离越大，时段划分结果的准确性就越高。

CP指标计算为

$ {x_{{\rm{CP}}}} = \frac{1}{L}\sum\limits_{i = 1}^L {\frac{1}{{{N_i}}}\sum\limits_{{x_j} \in {\mathit{\Omega }_i}} {\left\| {{x_j} - {v_i}} \right\|} } $

(10)

CP指标的值越小，说明类内聚类距离越近，时段划分结果越准确。

根据3.2节得到的最佳时段个数，分别利用K均值、FCM、WKM和SCFCM方法将3.1节中青霉素发酵过程测量数据(X(20×6×400))划分为5个时段，再代入式(8)~(10)计算时段划分点对应的DB和CP指标值，结果如表 3所示。

可以看出，无论从DB指标还是CP指标来看，本文方法的时段划分结果均优于其他4种方法。原因是基于K均值的方法从数据相似性出发实现时段划分，基于FCM的方法利用类间模糊距离约束实现分段，两种方法均存在时段划分结果依赖初始设定值、时序跳变等情况，时段划分结果的准确性较低；基于WKM和SCFCM的两种算法分别在K均值和FCM算法基础上增加了时序约束，但本质上依然受制于聚类方法依赖初始参数和对异常点敏感的缺点，影响了时段划分结果准确性的提升。本文方法在保证间歇过程数据完整性和结构特性的同时无需设定初始值，并充分考虑了批次间的时序性，因此时段划分结果具有更高的准确性。

4 结束语

本文提出一种基于PARAFAC2的多时段间歇过程时段划分方法。利用PARAFAC2可以对三维数据矩阵进行不展开建模的优势，在不破坏间歇过程三维数据完整性和结构特征的情况下，按照时序依次添加时间块矩阵，保证了时段划分结果的时序性和准确性，通过识别过程模型时间块矩阵控制限和时段评价划分指标(PPCI)的变化，获得了最优的时段划分结果。实验结果表明，与基于K均值、FCM、WKM、SCFCM的间歇过程时段划分方法相比，本文方法有效地提高了间歇过程时段划分的准确性。