CN111624979B - 一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法 - Google Patents

Abstract

本发明公开了一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法。本方法通过提取高维过程变量数据集中变化缓慢的特征，分析多变量闭环过程中可能出现的振荡信号与慢特征的共性，基于各慢特征的自相关函数曲线进行振荡指数的计算，对多个不同周期的振荡源进行提取和检测，并且结合振荡源重构，建立基于慢特征分析的多周期振荡检测与溯源模型。该方法不仅可以实现对多变量闭环控制系统的振荡检测，能有效提取和识别闭环系统中不同周期的多个振荡源，而且通过对振荡源的重构和溯源指标的设计，可以进一步实现对多周期振荡的溯源，判断出振荡来自于哪个控制回路，传播路径如何，完成对工业闭环系统的多周期振荡检测与溯源。

Description

一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源 方法

技术领域

本发明属于工业系统控制性能退化监测及溯源诊断，特别是针对多变量闭环控制系统多振荡共发情形的振荡检测与溯源方法。

背景技术

由于现代工业规模的增长和自动化的迅速发展，自动控制系统已经逐渐取代了基于工程师的控制策略，从而形成了自动化控制及工业过程的无人值守化管理和发展。在这种情况下，控制性能退化情况监测及溯源诊断受到了越来越多的关注，其目的是提供一种在线的自动化程序，以提供性能是否满足预期的信息。其中，最活跃的研究方向之一是振荡检测与溯源。因为振荡是工业过程中最常见的一种控制性能退化的情形，也是控制系统性能恶化的一种剧烈的表现形式，对于最终的产品质量和安全生产都有很大的负面影响。国际著名自动化整体解决方案供应商Honeywell公司曾针对流程工业中26000个控制回路进行调查统计分析，结果显示有16％左右的回路存在振荡。从某种意义上说，振荡是一种无法被噪声掩盖的周期性变化。振荡的原因可能是控制器设计不良，硬件非线性或外部干扰等。回路振荡会增加过程变量的波动，因此会导致产品质量降低，能源消耗增加以及产量和效益的降低。振荡还会缩短控制阀门的使用寿命，其增加的运行成本大致和振荡的幅值成正比。检测过程运行中的振荡和诊断振荡根源是非常重要的一项工作，因为工厂运行在产品质量界限附近比远离质量界限运行能产生更多的效益。并且，在一个有相互间作用的多回路系统中，一旦一个回路开始振荡，极有可能这个振荡会扩散到其他周围的回路当中。一个关键问题是如何在厂级的范围内检测多个振荡的存在并确定振荡的根源。

在过去的几十年时间里，振荡检测和溯源技术得到了广泛研究和发展，大量研究成果得到发表。方法大致可分为三类：(1)制定时域上的振荡检测标准：比如基于控制误差过零点之间的积分绝对误差(IAE)，基于自相关函数(ACF)的衰减比，和基于转移熵的方法等。(2)利用信号处理和分解技术：比如基于多元经验模态分解(MEMD)揭示振荡信号的谐波含量，基于非负矩阵分解确定线性独立振荡的基本形状，基于小波变和遗传算法的因式分解方法等。(3)频域分析方法：比如频谱包络方法，相干频谱法，频谱格兰杰因果关系法等。但是很少有方法涉及到对耦合控制回路中多振荡的检测，考虑到高维过程数据中的多振荡难以单一的方法对单个回路进行检测，因此应首先提取并检测所有可能的振荡。此外，多个振荡可能会通过具有直接物理连接，能量传递和预先设计的流通回路传播，导致振荡的快速传播，从而增加了溯源和诊断的难度。

本发明针对闭环控制系统的多振荡检测与溯源问题，从挖掘振荡变化特性的角度出发，提出了一种基于慢特征分析的闭环控制系统多振荡检测与溯源方法。实际上，慢特征分析很适合从高维过程数据中提取和恢复可能的多次振荡。首先，由于大多数控制器的稳定时间长而导致的低频滤波效应，工业控制回路中的大多数振荡具有相对较低的频率。这与慢特征分析的目标一致，该目标旨在最小化潜在特征的时间变化。因此，通过获得缓慢变化的特征的目标可以容易地进行振荡的提取。其次，慢特征分析可以提取从最慢到最快变化的特征，从而实现不同频段的分离。这提供了关于分离具有不同频率的多个振荡的优势。第三，慢特征分析的目标还可以解释为找到一组映射函数，这些映射函数的输出具有最大的单滞后自相关。由于某些噪声的单滞后自相关为零，因此可以证明慢特征分析在某种程度上减轻了噪声的影响。因此，本方法通过提取高维过程变量数据集中变化缓慢的特征，分析多变量闭环过程中可能出现的振荡信号与慢特征的共性，基于各慢特征的自相关函数曲线进行振荡指数的计算，对多个不同周期的振荡源进行提取和检测，并且结合振荡源重构，建立基于慢特征分析的多周期振荡检测与溯源模型。该方法不仅可以实现对多变量闭环控制系统的振荡检测，能有效提取和识别闭环系统中不同周期的多个振荡源，而且通过对振荡源的重构和溯源指标的设计，可以进一步实现对多周期振荡的溯源，判断出振荡来自于哪个控制回路，传播路径如何，完成对工业闭环系统的多周期振荡检测与溯源。到目前为止，尚未见到与本发明相关的研究报道。

发明内容

本发明的目的在于针对多变量工业闭环控制系统多振荡检测和溯源问题，提供了一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法。

本发明的目的是通过以下技术方案实现的：

一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法，包含以下步骤：

(1)选择闭环控制回路待检测的过程测量数据：对于一个多变量工业闭环控制系统，设其运行过程包含M个可测被控变量，在t时刻采样可以得到一个1×M的向量y(t)＝(y₁(t),y₂(t),…,y_M(t))，经过N次采样后得到正常过程下过程测量变量的数据矩阵Y_M＝(y(t),y(t+1),…,y(t+N))^T。对于M个可测被控变量的设定值，同样在t时刻进行采样得到一个1×M的向量p(t)＝(p₁(t),p₂(t),…,p_M(t))，同样经过N次采样后得到正常过程下被控变量设定值构成的数据矩阵P_M＝(p(t),p(t+1),…,p(t+N))^T。为了消除变量设定值变化的影响，将设定值数据矩阵和被控变量数据矩阵相减，得到矩阵

(2)输入数据时序拓展：考虑过程时延及动态性，对输入矩阵

进行d维时序扩展，得到输入矩阵

(3)数据标准化：对输入数据矩阵

按列减去该列的均值，并除以该列标准差进行标准化，得到标准化后的数据矩阵X。

(4)对步骤(3)中标准化后的过程数据矩阵X进行慢特征分析建模，获得变化缓慢的慢特征矩阵S。

这里采用线性慢特征分析，每一个慢特征都可以看作是输入变量的线性组合，所以从输入数据矩阵X映射到慢特征数据矩阵S(t)＝[s₁(t),s₂(t),…,s_M(t)]的表达式为：

S(t)＝W×X (1)

其中W＝[w₁,w₂,…,w_M]^T是映射矩阵，M为慢特征的个数，与输入变量个数相同。慢特征的映射矩阵的求解是一个优化问题，目标是使得经过映射后的特征尽可能地缓慢，即变化的速度最小，转换为数学表示即使特征的一阶差分Δ(s_i)最小。慢特征目标函数的优化可通过广义特征值分解进行求解：

AW＝BWΩ (2)

其中A表示输入数据矩阵X的一阶差分矩阵

的协方差矩阵，由

计算得到，B表示输入数据矩阵X的协方差矩阵<XX^T>_t。对应的，W＝[w₁,w₂,…,w_M]^T为特征向量构成的特征向量矩阵，M为输入数据矩阵的变量个数；Ω为对应特征向量矩阵的广义特征值构成的对角矩阵。求解获得映射矩阵后代入(1)式获得慢特征数据矩阵S。

(5)对慢特征矩阵中的所有慢特征进行缓慢度筛选，取J个慢特征进行振荡检测和相关性检验以去除冗余振荡源，最后提取到K个不同周期的回路振荡源，具体过程如下：

(5.1)利用(4.2)中获得的慢特征矩阵S，对每个特征向量s_j，计算其缓慢度

将缓慢度大于慢特征阈值T的J个慢特征进行保留，慢特征阈值可以根据具体情况进行调整，0＜T＜1，推荐选取0.5，丢弃剩下的M-J个特征。

(5.2)对每个保留的J个慢特征进行自相关函数计算，画出每个慢特征的自相关函数曲线：

k表示时刻的间隔，为正整数，k一般取到1/3N，

表示特征向量s_j的均值。

根据每个慢特征的自相关函数曲线，若其满足振荡周期大于等于4，则为其计算振荡指数：

其中，a是自相关函数曲线上从第一个峰值点到连接前两个最小极值点的直线距离，b是从第一个最小极值点到连接自相关函数曲线的起点和第一个峰值点的直线的距离。取J个慢特征中振荡指数大于振荡阈值R的特征为可能的振荡源，振荡阈值R一般设置为0.5，可根据具体情况进行调整。将提取出来的H个可能振荡源记为y_h(h＝1,2,...H)。

(5.3)随机选取第i个可能振荡源作为参照，计算它与其他由(5)保留得到的振荡特征之间的相关性系数：

其中，i,n∈h,i≠n。

若相关性系数c_i,n大于相关性阈值C，则说明这两个特征重复冗余，删除其中振荡指数较低的一个特征，保留其中振荡指数较大的特征，上述相关性阈值C一般设置为0.5。在留下的特征中重复执行上述相关性检验和冗余振荡源删除步骤，直到剩下的特征中任意两个特征的相关性不大于阈值C。经过此步保留得到K个慢特征，记为u_k,k＝1,2,...K，作为回路中真正存在的振荡源。

(6)根据(5.3)中得到的K个振荡源，利用各个回路经过时序拓展的输入数据进行振荡源重构，具体地，是将振荡源u_k在用每个回路的数据张成的子空间中进行重构：

得到重构振荡源记为

其中m＝1,2,...,M，即对于每个振荡源u_k，都有M个重构振荡源，用于后续溯源分析。

(7)对每一个振荡源u_k，需要通过溯源指标，从M个参与建模的控制回路中找到产生该振荡信号的回路。具体通过以下步骤实现：

(7.1)对于第k个振荡源u_k，k＝1,2,...K，有M个用不同回路数据重构得到的重构信号

其中m＝1,2,...,M。将振荡源u_k进行傅里叶变换，得到其在频谱上幅值最大处的频率为f_k，最大幅值为A_k。将重构信号也进行傅里叶变换，得到其在频谱上频率为f_k处的幅值，记为a_k。计算每个重构信号

与振荡源u_k之间的溯源指标：

即每个振荡源对应M个溯源指标。

(7.2)在第k个振荡源对应的M个溯源指标d_k,m中，将d_k,m(m＝1,2,...,M)按从大到小的顺序排列，取溯源指标值最大的m对应的回路作为振荡产生源头回路，并推出振荡传播路线即为d_k,m减小的路径。

(7.3)重复对K个不同的振荡源进行(7.1)和(7.2)中的振荡溯源操作，直到每一个振荡源都完成溯源。

进一步地，所述步骤4可通过两步奇异值分解(SVD)来求解慢特征，具体过程如下：

(4.1)首先对输入数据矩阵X进行白化处理去除其变量间的相关性。具体做法是，对B矩阵进行奇异值分解。B矩阵为输入矩阵X的协方差矩阵，表达式为<XX^T>_t，对其进行奇异值分解(SVD)得：

<XX^T>_t＝UΛU^T

其中，矩阵U是协方差矩阵<XX^T>_t分解得到的特征向量组成的矩阵，而Λ是一个对角阵，特征值是其每一个对角线上的元素，白化后得数据矩阵记为Z，表达式如下：

Z＝Λ^-1/2U^TX＝QX

白化后数据矩阵Z满足<ZZ^T>_t＝Q<XX^T>_tQ^T，其中<ZZ^T>_t是它的协方差矩阵。

(4.2)求解线性慢特征分析问题：具体的，可以看成是去寻找一个合适的矩阵P＝WQ^-1使得S＝P×Z，并且特征矩阵S应满足<SS^T>_t＝P<ZZ^T>_tP^T＝I，其中I是单位矩阵。

第二步奇异值分解是求得矩阵P的关键步骤，奇异值分解对象是白化后数据矩阵Z的差分矩阵

的协方差矩阵

分解后可求得：

P矩阵是特征向量矩阵，Ω是特征值构成的对角阵，其中特征值从小到大按序排列，然后可以通过下式计算得到慢特征分析的映射矩阵W为：

W＝PΛ^-1/2U^T

对应的慢特征矩阵S可以通过对输入矩阵X进行投影得到：

S＝WX＝PΛ^-1/2U^TX

本发明的有益效果为：本发明可以有效地将多变量工业闭环控制系统产生的高维数据集中所有可能的振荡信号提取出来，无需一一检测控制回路是否振荡，并且能准确追溯其振荡产生源头回路，便于现场工程师直接通过对振荡源所在回路进行运维就能够迅速排除振荡故障，实现闭环系统的自动振荡检测和溯源。

附图说明

图1是本发明的流程图；

图2是本方法具体应用的火力发电过程工艺流程图，图中，灰渣泵1，引风机2，除尘器3，省煤器4，空预器5，送风机6，汽包7，排粉机8，磨煤机9，高加10，给水泵11，除氧器12，低加13，凝泵14，凝汽器15，循泵16，高压缸17，中压缸18，低压缸19，发电机20，主变21，冷却塔22；

图3是本方法应用在火力发电过程一次风机回路卡涩引起振荡情况下的12个慢特征提取结果图；

图4是本方法应用在火力发电过程一次风机回路卡涩引起振荡情况下的12个慢特征自相关函数计算结果图；

图5是本方法应用在火力发电过程一次风机回路卡涩引起振荡情况下的振荡慢特征的振荡指数计算结果图；

图6是本方法应用在火力发电过程一次风机回路卡涩引起振荡情况下的12个重构特征和振荡源的对比图，图中实线表示重构振荡源，虚线表示振荡源；

图7是本方法应用在火力发电过程一次风机回路卡涩引起振荡情况下的溯源指标计算结果图；

具体实施方式

下面结合附图及具体实例，对本发明作进一步详细说明。

火力发电过程是一个大型复杂工业过程，如图2所示，主要由三大主要设备——锅炉、汽轮机、发电机及相应辅助设备组成，它们通过管道或线路相连构成生产主系统，即燃烧系统、汽水系统和电气系统。具体的工作过程是先将原煤由皮带输送到锅炉车间的煤斗，进入磨煤机系统进行研磨和干燥，制成具有一定细度和水分要求的煤粉然后与经过预热器预热的空气一起喷入炉内燃烧，烟气经除尘器清除灰分后由引风机抽出，经高大的烟囱排入大气。将煤粉送入锅炉后通过锅炉系统的加热将水转化为蒸汽，水在锅炉中加热后蒸发成蒸汽，经过加热器进一步加热，成为具有规定压力和温度的过热蒸汽，然后经过管道送入汽轮机，从而将煤的化学能转化为热能。在蒸汽涡轮机系统中，蒸汽不断膨胀，高速流动，冲击汽轮机的转子，以额定转速(旋转，将热能转换成机械能，燃烧的热能被蒸汽涡轮发电机转换成机械能，带动与汽轮机同轴的发电机发电，发电机将机械能转化为电能。

如图1所示，本发明方法包括以下步骤：

(1)数据选择:选择火力发电过程一次风机动叶开度卡涩案例的过程测量数据，包括共14个锅炉侧的控制变量，如表1所示，即M＝14。每次采样可以经过得到一个1×14的向量y(t)＝(y₁(t),y₂(t),…,y₁₄(t))，对于14个可测被控变量的设定值，同样每个时刻进行采样得到一个1×14的向量p(t)＝(p₁(t),p₂(t),…,p₁₄(t))，为了消除变量设定值变化的影响，将每个时刻的设定值数据和被控变量数据相减，即以控制误差(设定值-被控变量)作为过程输入变量，经过N次采样后得到输入数据矩阵

表1火力发电过程一次风机回路卡涩引起振荡案例的14个变量所在的控制回路

1	A侧一过减温水控制	7	B侧再热器减温水控制	13	B侧烟气挡板控制
						2	B侧一过减温水控制	8	B侧再热器减温水控制	14	燃油回油调节阀控制
3	A侧二过减温水控制	9	氧量校正
						4	B侧二过减温水控制	10	送风机动叶控制
5	B侧二过减温水控制	11	一次风机动叶控制
						6	A侧再热器减温水控制	12	A侧烟气挡板控制

(2)输入数据时序拓展：考虑过程时延及动态性，对输入矩阵

进行d维时序扩展，d＝40，得到输入矩阵

(3)数据标准化：对输入数据矩阵

按列减去该列的均值，并除以该列标准差进行标准化，得到标准化后的数据矩阵X。

(4)对(3)中标准化后的过程数据矩阵X进行慢特征分析建模，获得变化缓慢的慢特征矩阵S。

这里采用线性慢特征分析，每一个慢特征都可以看作是输入变量的线性组合，所以从输入数据矩阵X映射到慢特征数据矩阵S(t)＝[s₁(t),s₂(t),…,s_M(t)]的表达式为：

S(t)＝W×X (1)

其中W＝[w₁,w₂,…,w_M]^T是映射矩阵，M为慢特征的个数，与输入变量个数相同。慢特征的映射矩阵的求解是一个优化问题，目标是使得经过映射后的特征尽可能地缓慢，即变化的速度最小，转换为数学表示即使特征的一阶差分Δ(s_i)最小。慢特征目标函数的优化可通过广义特征值分解进行求解：

AW＝BWΩ (2)

其中A表示输入数据矩阵X的一阶差分矩阵

的协方差矩阵，由

计算得，B表示输入数据矩阵X的协方差矩阵<XX^T>_t。对应的，W＝[w₁,w₂,…,w_M]^T为特征向量构成的特征向量矩阵，M为输入数据矩阵的变量个数；Ω为对应特征向量矩阵的广义特征值构成的对角矩阵。本实施方式中通过两步奇异值分解(SVD)来求解慢特征，具体过程如下：

(4.1)首先对输入数据矩阵X进行白化处理去除其变量间的相关性。具体做法是，对B矩阵进行奇异值分解。B矩阵为输入矩阵X的协方差矩阵，表达式为<XX^T>_t，对其进行奇异值分解(SVD)得：

<XX^T>_t＝UΛU^T (3)

其中，矩阵U是协方差矩阵<XX^T>_t分解得到的特征向量组成的矩阵，而Λ是一个对角阵，特征值是其每一个对角线上的元素，白化后得数据矩阵记为Z，表达式如下：

Z＝Λ^-1/2U^TX＝QX (4)

白化后数据矩阵Z满足<ZZ^T>_t＝Q<XX^T>_tQ^T，其中<ZZ^T>_t是它的协方差矩阵。

(4.2)求解线性慢特征分析问题：具体的，可以看成是去寻找一个合适的矩阵P＝WQ^-1使得S＝P×Z，并且特征矩阵S应满足<SS^T>_t＝P<ZZ^T>_tP^T＝I，其中I是单位矩阵。

第二步奇异值分解是求得矩阵P的关键步骤，奇异值分解对象是白化后数据矩阵Z的差分矩阵

的协方差矩阵

分解后可求得：

P矩阵是特征向量矩阵，Ω是特征值构成的对角阵，其中特征值从小到大按序排列，然后可以通过下式计算得到慢特征分析的映射矩阵W为：

W＝PΛ^-1/2U^T (6)

对应的慢特征矩阵S可以通过对输入矩阵X进行投影得到：

S＝WX＝PΛ^-1/2U^TX (7)

(5)对慢特征矩阵中的所有慢特征进行缓慢度筛选，取J个慢特征进行振荡检测和相关性检验以去除冗余振荡源，最后提取到K个不同周期的回路振荡源，具体过程如下：

(5.1)利用(4.2)中获得的慢特征矩阵S，对每个特征向量s_j，计算其缓慢度

将缓慢度大于慢特征阈值T的J个慢特征进行保留，慢特征阈值可以根据具体情况进行调整，本实施例中，T取0.5，丢弃剩下的M-J个特征。

(5.2)对每个保留的J个慢特征进行自相关函数计算，画出每个慢特征的自相关函数曲线：

根据每个慢特征的自相关函数曲线，若其满足振荡周期大于等于4，则为其计算振荡指数：

其中，a是自相关函数曲线上从第一个峰值点到连接前两个最小极值点的直线距离，b是从第一个最小极值点到连接自相关函数曲线的起点和第一个峰值点的直线的距离。取J个慢特征中振荡指数大于振荡阈值R的特征为可能的振荡源，振荡阈值R一般设置为0.5，可根据具体情况进行调整。将提取出来的H个可能振荡源记为y_h(h＝1,2,...H)。

(5.3)随机选取第i个可能振荡源作为参照，计算它与其他振荡特征之间的相关性系数：

其中，i,n∈h,i≠n。

若相关性系数c_i,n大于相关性阈值C，则说明这两个特征重复冗余，删除其中振荡指数较低的一个特征，保留其中振荡指数较大的特征，上述相关性阈值C一般设置为0.5。在留下的特征中重复执行上述相关性检验和冗余振荡源删除步骤，直到剩下的特征中任意两个特征的相关性不大于阈值C。经过此步保留得到K个慢特征，记为u_k,k＝1,2,...K，作为回路中真正存在的振荡源。

(6)根据(5.3)中得到的K个振荡源，利用各个回路经过时序拓展的输入数据进行振荡源重构，具体地，是将振荡源u_k在用每个回路的数据张成的子空间中进行重构：

得到重构振荡源记为

其中m＝1,2,...,M，即对于每个振荡源u_k，都有M个重构振荡源，用于后续溯源分析。

(7)对每一个振荡源u_k，需要通过溯源指标，从M个参与建模的控制回路中找到产生该振荡信号的回路。具体通过以下步骤实现：

(7.1)对于第k个振荡源u_k，k＝1,2,...K，有M个用不同回路数据重构得到的重构信号

其中m＝1,2,...,M。将振荡源u_k进行傅里叶变换，得到其在频谱上幅值最大处的频率为f_k，最大幅值为A_k。将重构信号也进行傅里叶变换，得到其在频谱上频率为f_k处的幅值，记为a_k。计算每个重构信号

与振荡源u_k之间的溯源指标：

即每个振荡源对应M个溯源指标。

(7.2)在第k个振荡源对应的M个溯源指标d_k,m中，将d_k,m(m＝1,2,...,M)按从大到小的顺序排列，取溯源指标值最大的m对应的回路作为振荡产生源头回路，并推出振荡传播路线即为d_k,m减小的路径。

(7.3)重复对K个不同的振荡源进行(7.1)和(7.2)中的振荡溯源操作，直到每一个振荡源都完成溯源。

图3是本方法用于火电过程由一次风机回路卡涩引起厂级振荡案例的慢特征提取结果，包括12个慢特征，图4、5是慢特征的自相关函数曲线以及对应的振荡指数，其中自相关函数曲线上振荡周期不足4个周期的不满足计算振荡指数的要求，故在图5中前7个慢特征的振荡指数为0。表2是5个振荡特征的相关性检验结果，根据表可以看出，5个振荡特征全部相关，因此取振荡指数最大的特征为振荡源，即SF12，说明该数据集中存在一个振荡源，需要通过溯源方法定位到初始产生振荡的回路。图6是用各个回路的数据集对该振荡源进行重构得到振荡信号与振荡源的对比图。图7是溯源指标的结果图，根据图上指标最大值所在回路，可以判断出回路源头所在，即回路11(一次风机动叶控制回路)，可以看出对比实验案例的描述，振荡检测与溯源结果的判断是正确的。传统的振荡检测一般都对单个回路进行检测，而且无法很好地追溯起振荡起源，而通过本方法可以对厂级回路中的多个振荡进行提取和检测，通过溯源指标一一追溯起振荡源所在回路，有利于现场工程师快速调整故障回路。

表2在火力发电过程一次风机回路卡涩引起振荡案例振荡慢特征的相关性检验结果

	SF8	SF9	SF10	SF11	SF12
						SF8	1	0.62	/	/	/
SF9	/	1	0.67	/	/
						SF10	/	/	1	0.62	/
SF11	/	/	/	1	0.7
						SF12	/	/	/	/	1

Claims (2)

1.一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法，其特征在于，包含以下步骤：

(1)选择闭环控制回路待检测的过程测量数据：对于一个多变量工业闭环控制系统，设其运行过程包含M个可测被控变量，在t时刻采样可以得到一个1×M的向量y(t)＝(y₁(t),y₂(t),…,y_M(t))，经过N次采样后得到正常过程下过程测量变量的数据矩阵Y_M＝(y(t),y(t+1),…,y(t+N-1))^T；对于M个可测被控变量的设定值，同样在t时刻进行采样得到一个1×M的向量p(t)＝(p₁(t),p₂(t),…,p_M(t))，同样经过N次采样后得到正常过程下被控变量设定值构成的数据矩阵P_M＝(p(t),p(t+1),…,p(t+N-1))^T；将正常过程下过程测量变量的数据矩阵和正常过程下被控变量设定值构成的数据矩阵相减，得到矩阵

(2)输入数据时序拓展：考虑过程时延及动态性，对输入矩阵

进行d维时序扩展，得到输入矩阵

(3)数据标准化：对输入数据矩阵

按列减去该列的均值，并除以该列标准差进行标准化，得到标准化后的数据矩阵X；

(4)对步骤(3)中标准化后的过程数据矩阵X进行慢特征分析建模，获得变化缓慢的慢特征矩阵S：采用线性慢特征分析，将每一个慢特征看作是输入变量的线性组合，所以从输入数据矩阵X映射到慢特征数据矩阵S(t)＝[s₁(t),s₂(t),…,s_M(t)]的表达式为：

S(t)＝W×X (1)

其中W＝[w₁,w₂,…,w_M]^T是映射矩阵，M为慢特征的个数，与输入变量个数相同；慢特征的映射矩阵的求解是一个优化问题，目标是使得经过映射后的特征尽可能地缓慢，即变化的速度最小，转换为数学表示即使特征的一阶差分Δ(s_i)最小；慢特征目标函数的优化可通过广义特征值分解进行求解：

AW＝BWΩ (2)

其中A表示输入数据矩阵X的一阶差分矩阵

的协方差矩阵，由

计算得到，B表示输入数据矩阵X的协方差矩阵<XX^T>_t；对应的，映射矩阵W＝[w₁,w₂,…,w_M]^T为特征向量构成的特征向量矩阵，M为输入数据矩阵的变量个数；Ω为对应特征向量矩阵的广义特征值构成的对角矩阵；

求解获得映射矩阵后代入(1)式获得慢特征数据矩阵S；

(5)对慢特征矩阵中的所有慢特征进行缓慢度筛选，取J个慢特征进行振荡检测和相关性检验以去除冗余振荡源，最后提取到K个不同周期的回路振荡源，具体过程如下：

(5.1)利用获得的慢特征矩阵S，对每个特征向量s_j，计算其缓慢度

将缓慢度大于慢特征阈值T的J个慢特征进行保留，其中，0＜T＜1，丢弃剩下的M-J个特征；

(5.2)对保留的J个慢特征进行自相关函数计算，画出每个慢特征的自相关函数曲线：

k表示时刻的间隔，为正整数，

表示特征向量s_j的均值；

根据每个慢特征的自相关函数曲线，若其满足振荡周期大于等于4，则为其计算振荡指数：

其中，a是自相关函数曲线上从第一个峰值点到连接前两个最小极值点的直线的距离，b是从第一个最小极值点到连接自相关函数曲线的起点和第一个峰值点的直线的距离；从J个慢特征中提取振荡指数大于振荡阈值R的特征为可能的振荡源，振荡阈值R设置为0.5；将提取出来的H个可能振荡源记为y_h(h＝1,2,...H)；

(5.3)随机选取第i个可能振荡源作为参照，计算它与其他振荡源之间的相关性系数：

其中，i,n∈h,i≠n；

若相关性系数c_i,n大于相关性阈值C，则说明这两个特征重复冗余，删除其中振荡指数较低的一个特征，保留其中振荡指数较大的特征，上述相关性阈值C设置为0.5；在留下的特征中重复执行上述相关性检验和冗余振荡源删除步骤，直到剩下的特征中任意两个特征的相关性不大于阈值C；经过此步保留得到K个慢特征，记为u_k,k＝1,2,...K，作为回路中真正存在的振荡源；

(6)根据(5.3)中得到的K个振荡源，利用各个回路经过时序拓展的输入数据进行振荡源重构，具体为：将振荡源u_k用每个回路的数据进行重构：

得到重构振荡源记为

其中m＝1,2,...,M，即对于每个振荡源u_k，都有M个重构振荡源，用于后续溯源分析；

(7)对每一个振荡源u_k，需要通过溯源指标，从M个参与建模的控制回路中找到该振荡源的回路；具体通过以下步骤实现：

(7.1)对于第k个振荡源u_k，k＝1,2,…K，有M个用不同回路数据重构得到的重构信号

其中m＝1,2,…,M；将振荡源u_k进行傅里叶变换，得到其在频谱上幅值最大处的频率为f_k，最大幅值为A_k；将重构信号也进行傅里叶变换，得到其在频谱上频率为f_k处的幅值，记为a_k；计算每个重构信号

与振荡源u_k之间的溯源指标：

即每个振荡源对应M个溯源指标；

(7.2)在第k个振荡源对应的M个溯源指标d_k,m中，将d_k,m(m＝1,2,…,M)按从大到小的顺序排列，取溯源指标值最大的m对应的回路作为振荡产生源头回路，并推出振荡传播路线即为d_k,m减小的路径；

(7.3)重复对K个不同的振荡源进行(7.1)和(7.2)中的振荡溯源操作，直到每一个振荡源都完成溯源。

2.根据权利要求1所述的基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法，其特征在于，所述步骤4可通过两步奇异值分解(SVD)来求解慢特征，具体包括如下子步骤：

(4.1)首先对输入数据矩阵X进行白化处理去除其变量间的相关性；具体做法是，对B矩阵进行奇异值分解；B矩阵为输入矩阵X的协方差矩阵，表达式为<XX^T>_t，对其进行奇异值分解(SVD)得：

<XX^T>_t＝UΛU^T

其中，矩阵U是协方差矩阵<XX^T>_t分解得到的特征向量组成的矩阵，而Λ是一

个对角阵，特征值是其每一个对角线上的元素，白化后得数据矩阵记为Z，表达式如下：

Z＝Λ^-1/2U^TX＝QX

白化后数据矩阵Z满足<ZZ^T>_t＝Q<XX^T>_tQ^T，其中<ZZ^T>_t是它的协方差矩阵；

(4.2)求解线性慢特征分析问题：通过寻找一个矩阵P＝WQ^-1使得S＝P×Z，并且特征矩阵S应满足<SS^T>_t＝P<ZZ^T>_tP^T＝I，其中I是单位矩阵；

第二步奇异值分解是求得矩阵P的关键步骤，奇异值分解对象是白化后数据矩阵Z的差分矩阵

的协方差矩阵

分解后可求得：

P矩阵是特征向量矩阵，Ω是对应特征向量矩阵的广义特征值构成的对角矩阵，其中特征值从小到大按序排列，然后可以通过下式计算得到慢特征分析的映射矩阵W为：

W＝PΛ^-1/2U^T

对应的慢特征矩阵S可以通过对输入矩阵X进行投影得到：

S＝WX＝PΛ^-1/2U^TX。

Images