[发明专利]一种基于快速核独立成分分析的精馏塔故障检测方法

权利要求书

1.一种基于快速核独立成分分析的精馏塔故障检测方法，其特征在于，包括以下步骤：

首先，离线建模阶段包括如下所示步骤(1)至步骤(4)；

步骤(1)：利用精馏塔设备所安装的测量仪表，在精馏塔正常运行状态时采集N个样本数据x₁，x₂，…，x_N，其中第i个采样时刻的样本数据x_i∈R^m×1由m个采样数据组成，具体包括：塔釜液位、塔釜压力，塔釜底部产品流量，进料流量，进料温度，顶部回流流量，冷凝器液位，和各层塔板的温度；其中i∈{1，2，…，N}；

步骤(2)：对N个样本数据x₁，x₂，…，x_N实施标准化处理，得到N个m×1维的数据向量

步骤(3)：对N个数据向量进行筛选，得到n个聚类中心向量z₁，z₂，…，z_n，其中n小于N，具体的实施过程如下所示；

步骤(3.1)：初始化k＝0，j＝0，和c＝m；

步骤(3.2)：根据c均值聚类算法对N个数据向量进行聚类，得到c个聚类中心向量，分别记作z_k+1，z_k+2，…，z_k+c后，再设置k＝c与j＝j+1；

步骤(3.3)：判断是否满足条件j＜ε；若是，则设置c＝c+5后，返回步骤(3.2)；若否，则得到n个聚类中心向量z₁，z₂，…，z_n；其中ε表示c均值聚类算法的执行次数；

步骤(4)：利用z₁，z₂，…，z_n建立基于核独立成分分析的故障检测模型，保留核矩阵K，载荷矩阵A，矩阵φ，和控制上限Q_lim；

其次，离线建模阶段完成后，即可按照如下所示步骤不间断的对精馏塔实施在线故障检测；

步骤(5)：在最新采样时刻t，利用精馏塔设备所安装的测量仪表测量得到由m个采样数据组成的数据向量x_t∈R^m×1，并对其实施与步骤(2)中相同的标准化处理，得到新数据向量

步骤(6)：根据如下所示公式计算核向量k_t∈R^1×n中的第b个元素k_t(b)：

其中，b∈{1，2，…，n}，R^1×n表示1×n维的实数向量，δ为核参数，上标号T表示矩阵或向量的转置；

步骤(7)：根据如下所示公式对核向量k_t实施中心化处理得到

上式中，向量II_t∈R^1×n中所有元素都为1，矩阵II_n∈R^n×n中全部元素都是1；

步骤(8)：根据公式计算核独立成分向量s_t，再根据公式D_t＝s_tφs_t^T计算故障检测指标Q_t；

步骤(9)：判断是否满足条件：Q_t≤Q_lim；若是，则当前采样时刻精馏塔运行状态正常，返回步骤(5)继续实施对下一最新采样时刻的故障检测；若否，则执行步骤(10)从而决策精馏塔是否进入故障运行状态；

步骤(10)：返回步骤(5)继续实施对下一最新采样时刻样本数据的故障检测，若连续3个采样时刻的故障检测指标都不满足步骤(9)中的判断条件，则精馏塔进入故障工况；否则，返回步骤(5)继续实施对下一最新采样时刻的故障检测。

2.根据权利要求1所述的一种基于快速核独立成分分析的精馏塔故障检测方法，其特征在于，所述步骤(4)的具体实施过程如下所示：

步骤(4.1)：根据如下所示公式计算核矩阵K∈R^n×n中第a行第b列的元素K(a，b)：

其中，δ为核参数，||z_a-z_b||²＝(z_a-z_b)^T(z_a-z_b)，a∈{1，2，…，n}，b∈{1，2，…，n}，R^n×n表示n×n维的实数矩阵；

步骤(4.2)：根据如下所示公式对核矩阵K进行中心化处理得到矩阵

其中，矩阵II_n∈R^n×n中全部元素都是1；

步骤(4.3)：求解矩阵所有特征值所对应的特征向量，并按照大小对特征值进行降序排列得到λ₁≥λ₂≥…≥λ_n，与特征值λ₁，λ₂，…，λ_n对应的特征向量分别为p₁，p₂，…，p_n，再根据公式对特征向量p₁，p₂，…，p_n实施处理，其中b∈{1，2，…，n}；

步骤(4.4)：根据公式γ_d＝λ_d/(λ₁+λ₂+…+λ_n)计算特征值比例γ₁，γ₂，…，γ_n后，再确定出γ₁，γ₂，…，γ_n中大于阈值ζ的个数，该个数记为D，并将特征向量p₁，p₂，…，p_D组建成变换矩阵P＝[p₁，p₂，…，p_D]；

步骤(4.5)：根据如下所示公式计算核矩阵J∈R^N×n中第i行第b列元素J(i，b)：

其中，i∈{1，2，…，N}，b∈{1，2，…，n}，

步骤(4.6)：根据如下所示公式对J实施中心化处理得到后，再根据公式计算得分矩阵Z∈R^N×D，并初始化Z₀＝Z与g＝1：

上式中，矩阵II_N∈R^N×n中所有元素都是1，R^N×n表示N×n维的实数矩阵；

步骤(4.7)：将向量w_g设置为D×D维单位矩阵的第g列向量；

步骤(4.8)：根据如下所示公式更新向量w_g：

w_g＝E{Z₀^Tf(Z₀w_g)}-E{h((Z₀w_g))}w_g ⑦

其中，E{ }表示计算大括号内向量的均值，函数f(u)＝tanh(u)与函数g(u)＝sech(u)²，u为函数自变量；

步骤(4.9)：对向量w_g实施归一化处理w_g＝w_g/||w_g||后，判断向量w_g是否收敛；若否，则返回步骤(4.8)；若是，则执行步骤(4.10)；其中||w_g||表示计算向量w_g的长度；

步骤(4.10)：根据公式更新Z₀后，判断是否满足条件g＜D；若是，则设置g＝g+1后返回步骤(4.7)；若否，则将w₁，w₂，…，w_D与β₁，β₂，…，β_D分别组建成矩阵W＝[w₁，w₂，…，w_D]与矩阵B＝[β₁，β₂，…，β_D]；其中

步骤(4.11)：根据公式A＝PW(B^TW)^-1计算载荷矩阵A后，再根据公式计算核独立成分矩阵S∈R^N×D，并计算矩阵φ＝(S^TS)^-1；

步骤(4.12)：根据公式Q＝diag{SφS^T}计算故障检测指标向量Q∈R^N×1，再使用核密度估计法计算得到在置信度等于99％条件下的控制上限Q_lim，其中diag{ }表示将大括号内的矩阵对角线元素转变成列向量的操作运算；

步骤(4.13)：保留核矩阵K，载荷矩阵A，矩阵φ，和控制上限Q_lim。