基于熱擴散核密度確定密度峰值法的歷史工況識別

引 言

近年來，隨著物聯(lián)網、大數(shù)據和人工智能等技術的興起，數(shù)據驅動的方法在工業(yè)智能化的進程中扮演著重要角色。在實際生產過程中，原料性質、生產方案或操作條件等因素的變動將導致生產過程的多模態(tài)化[1]，如發(fā)酵過程[2]、冶金過程和鍋爐燃燒過程等，對其過程進行數(shù)字化時往往存在著非線性、多模態(tài)和變量間的強相關性等問題[3-4]。因此，深入研究多模態(tài)過程的特點對實際生產有著重要作用。通過獲取歷史工況特征，不僅可以為當前裝置選擇合適的工況模型及參數(shù)進行優(yōu)化，也能為生產決策提供重要的數(shù)據參考，如污水處理裝置的智能優(yōu)化、管道泄漏的自動化檢測和生產運行狀況的有效評估[5-6]等。

在對多模態(tài)過程的研究中，由于不同工況間存在著較大的差異，研究者通常假設每種工況下的過程數(shù)據近似服從一種高斯分布，運用主成分分析（PCA）、偏最小二乘（PLS）、獨立成分分析（ICA）和支持向量數(shù)據描述（SVDD）模型等方法提取工況數(shù)據的特征，然后建立模型應用于過程故障檢測、過程控制和過程優(yōu)化等[7-10]。由于每種工況下的數(shù)據具有相似性，有學者將數(shù)據聚類的方法用于多模態(tài)過程的特征提取[11]。常用的聚類方法包括模糊C均值法[12]、K-均值法[13]、高斯混合模型（GMM）[14-15]和隱馬爾可夫模型（HMM）[16]等，這些方法在獲取數(shù)據特征時具有一定的有效性，但仍存在一些無法避免的缺陷。如K-均值法需要事先確定聚類數(shù)量，對數(shù)據中的噪聲點敏感；模糊C均值法存在聚類數(shù)量和參數(shù)選取的問題；HMM模型需要事先知道各種模態(tài)的概率且固定不變；GMM模型在使用期望最大法求解時，存在計算量較大、對模型參數(shù)的初值敏感和容易陷入局部極值等問題，這些缺點都將導致無法準確地識別工況[17-18]。有學者對GMM模型進行深入研究，提出了給定模型參數(shù)初值[19]和基于信息準則確定聚類數(shù)量[20]的方法，其中F-J的方法較為著名[21-22]，它通過在迭代計算中不斷剔除冗余的高斯分量得出聚類結果，但是該方法需要一個較大的聚類數(shù)量導致計算量大且收斂困難，其結果的準確性也不能保證。

快速搜索發(fā)現(xiàn)密度峰[23]（CFSFDP）是基于局部密度的一種聚類技術，它根據聚類中心點密度較大且與其他中心點距離較遠的特點，引入高斯核密度估計函數(shù)（KDE）計算數(shù)據點的密度，再通過歐氏距離計算數(shù)據點間的距離，從而完成數(shù)據聚類。但是該方法的聚類效果取決于截距參數(shù)，為避免這一點，有學者對其進行改進并提出了無須事先確定截距參數(shù)的熱擴散核密度確定密度峰的技術[24]（CFSFDP-HD）。本文提出將CFSFDP-HD技術與GMM模型結合的方法，首先通過CFSFDP-HD方法對多模態(tài)過程數(shù)據進行聚類，然后將聚類結果作為GMM模型的初值，從而對多模態(tài)過程的工況進行較準確的估計。

1 工況識別方法

1.1 高斯混合模型

過程數(shù)據 X n×d 是d維的n個樣本數(shù)據，且X={x1,x2,?,xn}$X=\left(x_1,x_2,?,x_n\right)$，其概率密度函數(shù)可表示為：

p(x|θ)=∑i=1kτig(x|θi)$p\left(x|\theta \right)=\sum _{i=1}^k{\tau }_{i}g\left(x|{\theta }_i\right)$(1)

其中，k為高斯模型的數(shù)量，τi 和θi={μi,Σi}${\theta }_i=\left({\mu }_i,{\Sigma }_i\right)$分別為第i個高斯模型的權重和參數(shù)（平均值和協(xié)方差）。

第i個高斯模型對應的高斯密度函數(shù)為：

g(x|θi)=1(2π)d/2|Σi|1/2exp[?12(x?μi)TΣ?1i(x?μi)]                         (i=1,2,?,k)$\begin{array}{c}g\left(x|{\theta }_i\right)=\frac{1}{{\left(2\pi \right)}^{d/2}{\left({\Sigma }_i\right)}^{1/2}}exp\left(-\frac{1}{2}{\left(x-{\mu }_i\right)}^{{\rm T}}{\Sigma }_i^{-1}\left(x-{\mu }_i\right)\right)\\ \left(i=\mathrm{1,2},?,k\right)\end{array}$(2)

模型的參數(shù)θi 常用EM法[25]求解，通過不斷地更新后驗概率和模型參數(shù)，直到模型參數(shù)幾乎不變。針對數(shù)據X={x1,x2,?,xn}$X=\left(x_1,x_2,?,x_n\right)$和模型初始參數(shù)θ(0)={{τ(0)1,θ(0)1},{τ(0)2,θ(0)2},…,{τ(0)k,θ(0)k}}${\theta }^{\left(0\right)}=\left(\left({\tau }_1^{\left(0\right)},{\theta }_1^{\left(0\right)}\right),\left({\tau }_2^{\left(0\right)},{\theta }_2^{\left(0\right)}\right),\dots ,\left({\tau }_k^{\left(0\right)},{\theta }_k^{\left(0\right)}\right)\right)$，其迭代計算步驟如下。

E步驟：

P(s)(Ck|xj)=τ(s)kg(xj|μ(s)k,Σ(s)k)∑i=1kτ(s)ig(xj|μ(s)i,Σ(s)i)         (j=1,2,?,n)$P^{\left(s\right)}\left(C_k|x_j\right)=\frac{{\tau }_k^{\left(s\right)}g\left(x_j|{\mu }_k^{\left(s\right)},{\Sigma }_k^{\left(s\right)}\right)}{\sum _{i=1}^k{\tau }_i^{\left(s\right)}g\left(x_j|{\mu }_i^{\left(s\right)},{\Sigma }_i^{\left(s\right)}\right)}\left(j=\mathrm{1,2},?,n\right)$(3)

P(s)(Ck|xj)$P^{\left(s\right)}\left(C_k|x_j\right)$表示第j個樣本屬于第k個高斯模型的后驗概率，s表示第s次迭代。

M步驟：

μ(s+1)k=∑j=1nP(s)(Ck|xj)xj∑j=1nP(s)(Ck|xj)${\mu }_k^{\left(s+1\right)}=\frac{\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)x_j}{\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)}$(4)Σ(s+1)k=∑j=1nP(s)(Ck|xj)(xj?μ(s+1)k)(xj?μ(s+1)k)T∑j=1nP(s)(Ck|xj)${\Sigma }_k^{\left(s+1\right)}=\frac{\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)\left(x_j-{\mu }_k^{\left(s+1\right)}\right){\left(x_j-{\mu }_k^{\left(s+1\right)}\right)}^{{\rm T}}}{\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)}$(5)τ(s+1)k=∑j=1nP(s)(Ck|xj)n${\tau }_k^{\left(s+1\right)}=\frac{\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)}{n}$(6)

其中，μ(s+1)k、Σ(s+1)k、τ(s+1)k分別為第k個高斯模型在第（s+1）次迭代的平均值、協(xié)方差矩陣和先驗概率。

基于最短信息長度準則的F-J方法只需對式（6）進行如下修改，即可得到較為理想的聚類結果。

τ(s+1)k=max{0,(∑j=1nP(s)(Ck|xj))?v2}∑i=1kmax{0,(∑j=1nP(s)(Ck|xj))?v2}${\tau }_k^{\left(s+1\right)}=\frac{max\left(0,\left(\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)\right)-\frac{v}{2}\right)}{\sum _{i=1}^{k}max\left(0,\left(\sum _{j=1}^n{P}^{\left(s\right)}\left(C_k|x_j\right)\right)-\frac{v}{2}\right)}$(7)

其中，v=12d2+32d$v=\frac{1}{2}{d}^2+\frac{3}{2}d$，d為變量的個數(shù)，通過迭代將任意兩個相同的高斯模型進行合并，最終獲得多個工況模型及其參數(shù)。

1.2 熱擴散核密度確定密度峰技術

基于熱擴散的高斯核函數(shù)為：

P(di,dj,t)=1n∑j=1n12πt√e?(di?dj)22t$P\left(d_i,d_j,t\right)=\frac{1}{n}\sum _{j=1}^n\frac{1}{\sqrt[]{2\pi t}}{e}^{-\frac{{\left(d_i-d_j\right)}^2}{2t}}$ （j =1,2,…,n）(8)

P(di,dj,t)$P\left(d_i,d_j,t\right)$為樣本點i到j的轉移概率，t為核函數(shù)的帶寬，di - dj 為樣本j到i的距離。

估算任意樣本點i的概率密度函數(shù)為：

ρi=f?(d;t)≈∑k=0n?1αke?k2π2t/2cos(kπd)${\rho }_i=\stackrel{?}{f}\left(d;t\right)\approx \sum _{k=0}^{n-1}{\alpha }_k{e}^{-k^2{\pi }^{2}t/2}cos\left(k\pi d\right)$(9)

式（9）為KDE的完全自適應形式，考慮了最佳帶寬選擇和邊界校正。其中n為一個較大的整數(shù)，本文取n為樣本數(shù)量，αk 為：

αk=???????1,k=01n∑i=1ncos(kπdi),k=1,2,?,n?1${\alpha }_k=\left(\begin{array}{c}\begin{array}{cc}1,& k=0\end{array}\\ \begin{array}{cc}\frac{1}{n}\sum _{i=1}^{n}cos\left(k\pi d_i\right),& k=\mathrm{1,2},?,n-1\end{array}\end{array}\right)$(10)

最佳帶寬的選擇使用了改進的Sheather–Jones(ISJ)方法[26]，其計算步驟如下：

t=ξγ[l](t)$t=\xi {\gamma }^{\left(l\right)}\left(t\right)$(11)ξ=62√?37$\xi =\frac{6\sqrt[]{2}-3}{7}$(12)γ[l](t)=γ1(γ2(?γl(t)?))${\gamma }^{\left(l\right)}\left(t\right)={\gamma }_1\left({\gamma }_2\left(?{\gamma }_l\left(t\right)?\right)\right)$(13)γ[l](t)=1+0.5(l+0.5)31×3×?×(2l?1)nπ/2√∥∥f(l+1)∥∥2${\gamma }^{\left(l\right)}\left(t\right)=\frac{1+0.5^{\left(l+0.5\right)}}{3}\frac{1\times 3\times ?\times \left(2l-1\right)}{n\sqrt[]{\pi /2}{\left(f^{\left(l+1\right)}\right)}^2}$(14)∥∥f(l+1)∥∥2=∑k=1l?1(kπ)(2l+2)α2kexp(?(kπ)2t)${\left(f^{\left(l+1\right)}\right)}^2={\sum _{k=1}^{l-1}\left(k\pi \right)}^{\left(2l+2\right)}{\alpha }_k^{2}exp\left(-{\left(k\pi \right)}^{2}t\right)$(15)

其中，當l ≥ 5時，l的取值對式（11）的計算結果影響很小，故本文中取l = 5。

帶寬t的詳細求解步驟如下：

（1）設置一個較小的容差ε = 10-9，令yq=ε，q = 0；

（2）計算yq+1=ξγ[l](yq)$y_{q+1}=\xi {\gamma }^{\left(l\right)}\left(y_q\right)$；

（3）如果∣∣yq+1?yq∣∣<ε$\left(y_{q+1}-y_q\right)<\epsilon$，t = yq+1停止，否則yq = yq+1，q = q +1，返回步驟（2）。另外，令t = sqrt(t)/3.3，可對邊界點進行修正。

計算每一樣本點i到最近的高密度點j的距離：

δ={min(dij), if ? j, ρj>ρimax(dij), otherwise$\delta =\left(\begin{array}{c}min\left(d_{ij}\right),if?j,{\rho }_j>{\rho }_i\\ max\left(d_{ij}\right),otherwise\end{array}\right)$(16)

1.3 提出方法的計算步驟

本文提出的方法對近似服從高斯分布的未知多模態(tài)穩(wěn)態(tài)工況進行識別時，首先利用CFSFDP-HD技術對多模態(tài)過程數(shù)據進行聚類，確定聚類中心點及其個數(shù)（即工況個數(shù)），然后將每一類數(shù)據的平均值和協(xié)方差作為GMM模型的初值，迭代求出不同工況的特征參數(shù)。其計算過程如下：

（1）將數(shù)據標準化處理，求取參數(shù)αk；

（2）由參數(shù)αk 和式（11）~式（15）得出最佳帶寬t；

（3）由式（9）和式（16）的結果畫出決策圖，并由此完成聚類；

（4）將每一類的特征參數(shù)作為GMM模型初值，求出最終工況參數(shù)。

通過以上步驟即可完成對歷史工況的準確識別，下面通過第2節(jié)中的兩個例子對該方法進行驗證。

圖1

圖1   基于熱擴散核密度的工況識別方法流程圖

Fig.1   Flow chart of recognizing operating modes based on kernel density estimation of heat diffusion

2 方法驗證與結果分析

2.1 仿真數(shù)據

根據文獻[27]中的多模態(tài)仿真模型生成過程數(shù)據，然后分別用本文提出的方法、K-均值法和GMM（F-J）的方法進行工況識別，數(shù)據生成模型如式（17）所示：

?????x1=0.5768s1+0.3766s2+e1x2=0.7382s21+0.0566s2+e2x3=0.8291s1+0.4009s22+e3$\left(\begin{array}{c}{x}_1=0.5768s_1+0.3766s_2+e_1\\ x_2=0.7382s_1^2+0.0566s_2+e_2\\ x_3=0.8291s_1+0.4009s_2^2+e_3\end{array}\right)$(17)x?i=xi,j?xi,minxi,max?xi,min${\stackrel{?}{x}}_i=\frac{x_{i,j}-x_{i,min}}{x_{i,max}-x_{i,min}}$(18)

式（18）為數(shù)據標準化的方法，xi,j 表示第i個變量的第j個數(shù)據，x?i${\stackrel{?}{x}}_i$表示標準化處理后的第i個變量的向量數(shù)據，xi,min表示第i個變量的最小值，xi,max表示第i個變量的最大值。

其中,e1~e3是服從[0，0.01]的高斯白噪聲分布，通過調整s1和s2的參數(shù)，生成含3個變量（x1、x2和x3）的多模態(tài)過程數(shù)據。其中模態(tài)1是變量s1和s2分別服從高斯分布為[20，0.8]、[1，1.3]得到的300個數(shù)據；模態(tài)2是變量s1和s2分別服從高斯分布[5，0.6]、[20，0.7]得到的300個數(shù)據；模態(tài)3是變量s1和s2分別服從高斯分布[16，1.5]、[20，0.7]得到的300個數(shù)據；模態(tài)4是和模態(tài)2在相同參數(shù)（工況）下產生的300個數(shù)據，用于檢驗三種方法能否準確地獲取實際的工況狀態(tài)。

將生成數(shù)據用式（18）進行標準化處理，其分布情況如圖2所示，可以看出數(shù)據有四個階段，其中300~600和900~1200階段的狀態(tài)相同，然后分別對三種方法進行驗證。使用本文方法畫出關于密度和距離的決策圖，見圖3，圖中有三個中心點，表明本文方法根據過程數(shù)據識別出三種工況，每個中心點表示一種工況的數(shù)據中心。然后根據數(shù)據點到中心點的距離將其分類，將每一類的結果作為GMM模型的初值，從而得出不同工況的特征參數(shù)。使用K-均值法進行工況識別時，由于聚類數(shù)量是未知的，所以將聚類數(shù)量分別設置為3、4和5，其中K = 5時的結果與實際相差較大，其結果未在表1中列出。使用GMM（F-J）方法時，需要設置初始聚類數(shù)量（K）大于實際工況數(shù)量，本文分別設為4、5和6，GMM模型的初值設置為將過程數(shù)據平均分成K份，每份數(shù)據的特征參數(shù)[19]，先驗概率設為1/K，其中K = 6時的結果與K = 5時的結果幾乎相同，其結果未在表1中列出。三種方法的工況識別結果見表1?？梢钥闯霰疚姆椒ǐ@取的多模態(tài)過程的工況個數(shù)及其特征參數(shù)（變量的平均值和工況的先驗概率）與實際值一致。當K-均值法的聚類數(shù)量與實際工況數(shù)量一致時（K=3），得到的工況特征參數(shù)與實際值的相對偏差較小，當聚類數(shù)量大于實際工況數(shù)量時（K=4），得到的工況特征參數(shù)與實際值的相對偏差較大，由此看出該方法的工況識別效果取決于聚類數(shù)量的準確選擇。GMM（F-J）方法給定不同的初始聚類數(shù)量（K=4、5）時均將工況識別為4種，未能準確識別出實際工況的個數(shù)，但得到的工況特征參數(shù)與實際值的偏差在0.01~-20.39，其結果仍具有一定的參考價值。

圖2

圖2   仿真多模態(tài)過程數(shù)據標準化

Fig.2   Normalization of multi-modal process simulation data

圖3

圖3   仿真多模態(tài)過程數(shù)據的聚類中心決策圖

Fig.3   Clustering center decision diagram of process data for simulating multiple operating modes

表1   仿真多模態(tài)過程的工況識別結果

Table 1  Recognition results of simulation multiple operating modes

2.2 TE過程

Tennessee Eastman（TE）工業(yè)過程是由美國Eastman化學品公司開發(fā)的復雜工業(yè)過程的仿真平臺，它包括六種工作模態(tài)，每種模態(tài)具有不同的產品比例（G/H），該流程包含12個操作變量、22個連續(xù)過程測量變量和19個組成測量變量[28-30]。本文選取TE過程中模態(tài)1~模態(tài)4作為多模態(tài)過程，選取41個測量變量作為工況識別的變量，其中每種模態(tài)取300個數(shù)據為1組，第5組和第3組為相同模態(tài)下的數(shù)據，具體模態(tài)選取情況見表2。

表2   TE過程的模態(tài)選取情況

Table 2  Mode selection of the TE process

將過程數(shù)據用式（18）進行標準化處理，選取反應器溫度、A和C的混合進料量、產品分離器壓力和回收流量4個變量畫出分布圖，見圖4，可以看出5個階段中600~900和1200~1500的狀態(tài)相同，然后分別用三種方法進行驗證。本文方法得到的決策圖見圖5，圖中有4個中心點，表明本文方法識別出4種工況，然后使用GMM法獲取每種工況的特征參數(shù)。將K-均值法的聚類數(shù)量（K）分別設為4、5和6，其中K = 5和6的結果與實際相差較大，其結果未在表4中列出。將GMM（F-J）方法的初始聚類數(shù)量設為4、5和6時，該方法均無法獲取工況參數(shù)。三種方法工況識別結果的典型數(shù)據見表3、表4。

圖4

圖4   4個TE過程變量的標準化

Fig.4   Normalization of 4 TE process variables

圖5

圖5   TE多模態(tài)過程數(shù)據的聚類中心決策圖

Fig.5   Clustering center decision diagram of TE multi-modal process

表3   TE多模態(tài)過程的工況個數(shù)及先驗概率的識別結果

Table 3  Recognition results of the number and probability of the TE multi-modal process

表4   TE多模態(tài)過程變量的識別結果

Table 4  Recognition results of TE multi-modal process variables

從表3可以看出本文方法得到的歷史工況的個數(shù)和先驗概率與實際值一致；K-均值法的結果則取決于設定的聚類數(shù)量K，當K與實際一致（K = 4）時也可以較準確獲取歷史工況的個數(shù)及先驗概率，但是當K = 5和K = 6時，其結果與實際相差較大。GMM（F-J）法則無法獲取到工況的參數(shù)。

三種方法過程變量的識別結果見表4，可以看出本文提出的方法識別結果的平均相對偏差在 -0.0043~-1.3681，最大相對偏差為-4.748，最小相對偏差為-0.0014；K-均值法識別結果的平均相對偏差在-0.0043~-1.4371，最大相對偏差為-4.9847，最小相對偏差為-0.0014。結合表3、表4可以看出GMM（F-J）法不適合本案例的工況識別，本文方法和給定準確聚類數(shù)量的K-均值法都可以較準確地識別出工況特征，但K-均值法的準確性依賴于聚類數(shù)量的選擇，而本文方法則沒有這種約束。

3 結 論

針對目前工況識別方法的不足，提出將人工智能領域的CFSFDP-HD技術與GMM模型結合用于對多模態(tài)過程的歷史工況進行識別的方法，避免了K-均值法需要預先提供準確聚類數(shù)量的缺點，并利用案例對本文所提方法進行了驗證，結果表明：GMM（F-J）法不能保證準確地識別工況，K-均值法只有在給定正確工況數(shù)量的前提下才能獲得較好的結果，而本文方法則可方便、有效地對歷史多工況進行準確識別，具有更強的實用性。

符 號 說 明