基于粒子群優化粒子濾波和CUDA加速的故障診斷方法

曹 潔 李 釗 王進花 余 萍

1(蘭州理工大學電氣工程和信息工程學院 甘肅 蘭州 730050)2(甘肅省制造業信息化工程研究中心 甘肅 蘭州 730050)

0 引 言

隨著以風力發電為代表的清潔能源快速發展，風電機組的裝機規?？焖僭鲩L。如何保證風電機組安全穩定的運行，已成為不可忽視的問題。變槳距系統作為風電機組運行的關鍵，能夠保證高風能利用系數、減輕干擾載荷和風力制動。因此，對變槳距系統進行故障診斷可以有效地避免重大事故的發生，提高運行可靠性，降低運營成本[1]。但受環境影響，變槳距系統存在各類信號擾動和隨機噪聲，導致故障診斷的準確性不高[2]。

粒子濾波(Particle Filter，PF)由于其適用于非線性非高斯系統，在故障診斷領域得到了廣泛的研究和應用[3]。但是算法存在粒子退化和實時性問題，不利于故障診斷的快速響應。針對粒子退化問題，將粒子群優化算法(Particle Swarm Optimization，PSO)引入粒子濾波算法中得到粒子群優化粒子濾波算法，使得粒子向狀態后驗概率密度高值區移動，改善粒子退化問題[4]。針對實時性問題，文獻[5-6]提出了基于圖像處理單元(GPU)加速的并行粒子濾波算法，提高了算法的執行速度。文獻[7]基于GPU存儲架構提出一種共享內存系統重采樣，顯著提升了粒子濾波的運行速度。文獻[8]提出了基于GPU的粒子濾波并行算法，并應用于目標跟蹤，在不降低跟蹤準確性的同時，提升了跟蹤算法的速度。文獻[9]提出了GPU加速的粒子群優化粒子濾波算法，提升了算法的速度。文獻[10]為解決移動機器人故障診斷實時性問題，提出一種并行化粒子濾波器，減少了故障診斷的反應時間。文獻[11]使用CUDA(compute unified device architecture)在GPU上加速粒子濾波和輔助粒子濾波，并應用于工業控制，獲得了相比物理過程兩個數量級的加速。由于粒子間的相互作用，重采樣方法的有效并行實現是一個挑戰。文獻[12-13]提出了局部重采樣，將全部的粒子劃分為多個子集，減少重采樣對權值的全局操作，從而提升算法的并行性。文獻[14-15]分析了Metropolis重采樣和拒絕重采樣在GPU上的并行實現，并將這兩種方法與多項式、分層和系統重采樣進行比較，獲得了更快的執行速度。但是以上方法存在全局內存的非合并訪問問題，阻礙了重采樣有效并行化，不利于算法執行效率。

針對上述問題，本文設計并實現了一種基于CUDA的PSOPF并行算法，同時提出了一種合并訪問全局內存的拒絕重采樣算法，充分利用GPU的計算性能，提高算法的執行速度。在此基礎上，利用PSOPF并行算法結合殘差平滑法對風電機組變槳距系統進行故障診斷。

1 粒子群優化粒子濾波

1.1 粒子濾波原理

假設系統動態空間模型描述如下：

狀態方程：

xk=f(xk-1)+wk

(1)

量測方程：

yk=h(xk)+vk

(2)

式中：xk為k時刻系統的狀態值；yk為k時刻的系統觀測值；f和h為非線性函數；wk和vk分別為系統噪聲和觀測噪聲。

粒子濾波算法是一種基于蒙特卡洛仿真的近似貝葉斯濾波算法，通過從后驗概率中抽取的隨機粒子來表達其分布：

(3)

(4)

當式(3)和式(4)隨著測量值遞推更新時，可通過粒子集的加權平均近似后驗概率密度。

1.2 算法步驟

粒子群優化算法由Eberhart和Kennedy于1995年提出，其原理是利用群體中個體的運動來探尋問題求解空間。群體中的每個個體都將自己的經驗貢獻給群體的發展，并利用群體的全局經驗進行自身的發展。因此，信息在兩個方向進行傳遞，分別為從群體到個體和從個體到群體。算法以隨機初始化的粒子群作為起始，通過迭代尋優使粒子群找到最優解。在每次迭代中，每個粒子的速度Vi(t)和位置Xi(t)根據先前個體的最優解XPbest(t)和群體的最優解XGbest(t)進行更新：

Vi(t+1)=wVi(t)+c1r1(XPbest(t)-Xi(t))+

c2r2(XGbest(t)-Xi(t))

(5)

Xi(t+1)=Xi(t)+Vi(t)

(6)

式中：i=1,2,…,N，N為群體的粒子數；t為當前迭代次數；w為慣性系數；c1和c2稱為學習因子，是非負常數；r1和r2為介于(0，1)區間的隨機數。

粒子群優化粒子濾波算法利用PSO的尋優能力對粒子采樣過程進行優化，使粒子集趨向于后驗概率密度較高的區域運動，并將最新的觀測值引入采樣過程，定義適應度函數為：

(7)

式中：R為觀測噪聲；ynew為最新測量值；ypre為預測值。具體步驟如下：

步驟2執行PSO算法。

① 隨機初始化粒子的位置和速度。

② 初始化個體最優XPbest(t)和全局最優XGbest(t)。

③ 根據式(7)計算每個粒子的適應度值f(xi)。

④ 對比和更新個體最優XPbest(t)和全局最優XGbest(t)。對每個粒子，將其適應度值與個體最優XPbest(t)相比較，若較好，則替換個體最優。再將每個粒子的個體最優XPbest(t)與群體的全局最優XGbest(t)相比較，若較好，則替換全局最優。

⑤ 根據式(5)和式(6)更新粒子每個粒子的狀態。

⑥ 判斷是否滿足結束條件，若不滿足，則返回到③。

步驟3計算粒子的重要性權值并歸一化。

(8)

步驟4重采樣。

步驟6計算狀態估計值：

(9)

2 基于CUDA的PSOPF算法

2.1 CUDA并行計算架構

CUDA是采用基于單指令多線程(Single Instruction Multiple Thread，SIMT)的細粒度并行模型[16]，其基本編程模型如圖1所示。其中：CPU作為主機，負責處理串行任務和流程控制；GPU作為協處理器，負責執行并行化任務。CUDA通過調用核函數(kernel)執行并行計算任務，并將線程組織成線程塊(block)和線程網格(grid)形式。線程通過線程束(warp)調度執行，一個線程束由連續的32個線程組成，每個線程束中的線程同時執行。線程可以從多個內存空間訪問數據：私有寄存器，一個線程塊內所有線程擁有的共享內存，以及GPU上的全局內存。此外，所有線程都可以訪問另外兩個只讀存儲空間：常量內存和紋理內存。所有這些內存類型都有不同級別的延遲和帶寬，因此需要正確選擇要使用的內存級別。

圖1 CUDA編程模型

全局內存的訪問以32、64或128字節為單位，并且以32字節為基準對齊。當線程束中所有線程訪問連續且對齊的內存空間時，則所有訪問可以合并為一次內存讀取。如果內存地址不在同一內存區段中，則無法合并訪問，內存訪問性能急劇下降。并且，當同一線程束中的內存訪問不連續時，盡管線程不需要這些訪問之間的數據，但這些數據仍然會被傳輸，浪費了內存帶寬。因此，為了獲得好的數據吞吐量，需要通過合并訪問將大量的內存獲取請求合并，減少內存獲取次數。

2.2 基于CUDA的拒絕重采樣

算法1拒絕重采樣算法

for i=1:N

j=i

u～U[0,1]

while u>wj/sup w

j～U[1,…,N]

u～U[0,1]

end for

其中：xnew是重新采樣后的新粒子集；u是服從0到1之間均勻分布的隨機實數；j是服從1到N之間均勻分布的隨機整數；j是在k時刻所選粒子的索引，該粒子是在k+1時刻的粒子i復制的祖先。

拒絕重采樣可以作為單個CUDA內核函數調用實現，使每一個線程對應一個粒子，線程為其代表的粒子選取和復制祖先。拒絕重采樣在隨機抽選粒子時，存在非合并的全局內存訪問。對于較大的粒子數，非合并訪問次數增多，導致算法執行效率降低?？紤]到線程束中線程在同一內存區段上執行讀/寫操作能夠避免非合并的內存訪問，提出改進的拒絕重采樣算法，其偽代碼如下：

算法2改進的拒絕重采樣算法

for i=1:N

j=i

s～U[1,…,Scount]

while u>wj/sup w

j～U[(s-1)*Ecount+1,…,s*Ecount]

u～U[0,1]

end for

其中：s為全局內存中一組固定數量的連續區段，最小的s為單個內存區段，最大的區段為跨越權值數組的所有區段的集合。設Esize是一個權值所占字節大小，Ssize為一個s區段所占字節大?。籗count是s區段的數量，則Scount=Esize×N/Ssize。Ecount是一個s區段中權值的數量，則Ecount=Ssize/Esize。j是服從的s段中第一到最后一個元素索引之間均勻分布的隨機整數；s服從第一個到最后一個s區段索引之間均勻分布的隨機整數。為方便線程和全局內存對齊，粒子數N的值設為2的次冪。

2.3 基于CUDA的PSOPF并行算法

通過對PSOPF算法的并行性分析，采取線程與粒子一一對應的思路，實現基于CUDA的PSOPF并行算法。具體實施思路為：創建與粒子數相同的線程數，并分配相應的存儲空間，每個線程執行一個粒子的運動計算。PSOPF算法流程中更新全局最優和計算狀態估計值會產生粒子間的數據關聯，無法完全并行化。本文采用并行規約算法尋找極值和求和運算，每個線程處理一個粒子，提高并行度，以充分發揮GPU并行計算性能。根據PSOPF算法的計算流程，制定的基于CUDA的PSOPF并行算法流程如圖2所示。

圖2 基于CUDA架構的PSOPF并行算法流程

3 實驗結果及分析

3.1 算法性能測試

為驗證本文算法的精度和運行時間等性能，采用一維非線性系統模型在CPU和GPU平臺上進行實驗仿真，并與PF以及PSOPF算法進行對比。實驗所采用的計算平臺如表1所示。

表1 計算平臺

狀態模型和量測模型如下：

(10)

(11)

式中：系統噪聲wk為符合N(0,10)分布的高斯噪聲；測量噪聲vk為符合N(0,1)分布的高斯噪聲。均方根誤差公式為：

(12)

令粒子數N=256，單次仿真結果如圖3和圖4所示。分別令粒子數N=128和N=256，進行50次蒙特卡洛仿真并取平均值，得到三種算法的均方根誤差對比如表2所示。

圖3 不同算法的狀態估計結果

圖4 不同算法的RMSE對比結果

表2 不同算法RMSE對比

從圖3、圖4和表2可以看出，由于增加了PSO算法的尋優過程，引入粒子群優化的粒子濾波算法的均方根誤差明顯小于標準PF算法。而采用改進拒絕重采樣的本文算法跟蹤效果最好，均方根誤差小于PSOPF算法，說明其估計精度最高。

在不同粒子數時本文算法的運行時間如表3所示。

表3 本文算法運行時間對比

由表3可以看出，GPU的運行時間明顯小于CPU。算法在GPU上的運行時間始終小于1 s，而CPU的運行時間可達數秒，因此基于CUDA執行的本文算法能夠顯著地提高運行速度。并且隨著粒子數的增加，運行時間逐漸增加，GPU的增加幅度明顯小于CPU。在粒子數較多時，能夠獲得較好的計算性能。

在不同粒子數時分別在CPU和GPU上執行PSOPF算法和本文算法，獲得加速比曲線如圖5所示。

圖5 不同算法的加速比對比

由圖5可見，隨著粒子數從28增加到214，本文算法的加速比增長快于PSOPF算法，這是由于隨著粒子數的增加，算法計算量增加，采用改進拒絕重采樣的本文算法避免了非合并訪問，提高了GPU內存使用效率，顯著地提高了算法的實時性。

3.2 變槳距系統的故障診斷仿真與分析

風電機組變槳距系統由三個相同且獨立的槳距控制器組成，每個槳距角控制器被建模為實際的槳距角β與其參考值βref間的閉環傳遞函數。βref為閉環傳遞函數的輸入，由系統的槳距角控制器提供；β為傳遞函數的輸出，作為槳距角的測量值。實質上，變槳距系統是一個活塞伺服系統，可以通過二階傳遞函數建模[17]：

(13)

(14)

表4 故障模式表

變槳距系統的故障診斷原理如圖6所示，采用本文提出的改進粒子濾波算法，對應模式M0-M3，分別為每個模式配置一個粒子濾波器進行狀態估計。計算各個模式的槳距角估計值y(m)(m=0,1,2,3)，對比系統的實際測量值y獲得殘差e(m)，并根據式(15)計算最近M時刻內的殘差平滑值d(m)。

(15)

圖6 變槳距系統故障診斷原理框圖

首先進行故障檢測，設置殘差閾值δ(δ>0)，當殘差平滑值d(0)大于閾值δ時，則判斷系統發生故障，否則判斷系統無故障。檢測到故障發生后，通過比較每個故障模式的殘差平滑值d(m)(m=1,2,3)的大小來進行故障隔離，若其中某一個模式的殘差平滑值最小，則判斷系統發生了該模式故障。

根據正常模式的系統模型計算當前時刻的殘差平滑值d(0)，當d(0)大于閾值時，則檢測到系統發生故障。系統正常運行和發生故障M1-M3時的殘差曲線如圖7所示。

(a) 系統正常運行時的槳距角殘差

(b) 系統發生故障M1時的檢測結果

(c) 系統發生故障M2時的檢測結果

(d) 系統發生故障M3時的檢測結果圖7 不同故障的診斷結果

從圖7(a)可見，當系統正常運行時，槳距角殘差基本為零。從圖7(b)-(d)可以看出，對于故障M1-M3，系統無故障時的統槳距角殘差幅值較小，且變化比較平穩；當故障發生后，殘差變化劇烈且幅度高于閾值，因此判斷故障發生。

檢測到故障后，啟動故障M1-M3模型進行故障隔離，結果如圖8所示。

(a) 系統發生故障M1時的隔離結果

(b) 系統發生故障M2時的隔離結果

(c) 系統發生故障M3時的隔離結果圖8 不同故障的隔離結果

由圖8(a)可見，當檢測到故障M1發生后，由于M3為槳距角傳感器偏差故障，起始殘差在0.3左右。M1和M2的殘差曲線較為相似，但M1的殘差是最小的，故判斷為M1故障。同理分析圖8(b)和(c)可見，當系統實際狀態與某故障模型匹配時，其殘差平滑值相比其他故障模型較小，則判斷該故障模式發生。

4 結 語

本文基于CUDA架構設計實現了一種GPU加速的PSOPF并行算法，通過分析PSOPF算法的并行性，采取線程和粒子一一對應的方式實現并行化算法，運行速度得到顯著提升。利用改進的拒絕重采樣并行算法，解決全局內存讀取時出現的非合并訪問問題，提高了PSOPF并行算法在GPU上的執行效率，并將其應用到風電機組變槳距系統故障診斷中。實驗結果表明，所提出的改進的PSOPF并行算法有效地提高了實時性，是一種進行故障診斷有效的方法。

下一步工作將設計并實現基于CUDA的多模型粒子濾波故障診斷方法，通過故障診斷的并行化，以提高故障診斷的實時性能。