曹 潔 李 釗 王進(jìn)花 余 萍

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

0 引 言

隨著以風(fēng)力發(fā)電為代表的清潔能源快速發(fā)展，風(fēng)電機(jī)組的裝機(jī)規(guī)?？焖僭鲩L(zhǎng)。如何保證風(fēng)電機(jī)組安全穩(wěn)定的運(yùn)行，已成為不可忽視的問(wèn)題。變槳距系統(tǒng)作為風(fēng)電機(jī)組運(yùn)行的關(guān)鍵，能夠保證高風(fēng)能利用系數(shù)、減輕干擾載荷和風(fēng)力制動(dòng)。因此，對(duì)變槳距系統(tǒng)進(jìn)行故障診斷可以有效地避免重大事故的發(fā)生，提高運(yùn)行可靠性，降低運(yùn)營(yíng)成本[1]。但受環(huán)境影響，變槳距系統(tǒng)存在各類信號(hào)擾動(dòng)和隨機(jī)噪聲，導(dǎo)致故障診斷的準(zhǔn)確性不高[2]。

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

針對(duì)上述問(wèn)題，本文設(shè)計(jì)并實(shí)現(xiàn)了一種基于CUDA的PSOPF并行算法，同時(shí)提出了一種合并訪問(wèn)全局內(nèi)存的拒絕重采樣算法，充分利用GPU的計(jì)算性能，提高算法的執(zhí)行速度。在此基礎(chǔ)上，利用PSOPF并行算法結(jié)合殘差平滑法對(duì)風(fēng)電機(jī)組變槳距系統(tǒng)進(jìn)行故障診斷。

1 粒子群優(yōu)化粒子濾波

1.1 粒子濾波原理

假設(shè)系統(tǒng)動(dòng)態(tài)空間模型描述如下：

狀態(tài)方程：

xk=f(xk-1)+wk

(1)

量測(cè)方程：

yk=h(xk)+vk

(2)

式中：xk為k時(shí)刻系統(tǒng)的狀態(tài)值；yk為k時(shí)刻的系統(tǒng)觀測(cè)值；f和h為非線性函數(shù)；wk和vk分別為系統(tǒng)噪聲和觀測(cè)噪聲。

粒子濾波算法是一種基于蒙特卡洛仿真的近似貝葉斯濾波算法，通過(guò)從后驗(yàn)概率中抽取的隨機(jī)粒子來(lái)表達(dá)其分布：

(3)

(4)

當(dāng)式(3)和式(4)隨著測(cè)量值遞推更新時(shí)，可通過(guò)粒子集的加權(quán)平均近似后驗(yàn)概率密度。

1.2 算法步驟

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

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為群體的粒子數(shù)；t為當(dāng)前迭代次數(shù)；w為慣性系數(shù)；c1和c2稱為學(xué)習(xí)因子，是非負(fù)常數(shù)；r1和r2為介于(0，1)區(qū)間的隨機(jī)數(shù)。

粒子群優(yōu)化粒子濾波算法利用PSO的尋優(yōu)能力對(duì)粒子采樣過(guò)程進(jìn)行優(yōu)化，使粒子集趨向于后驗(yàn)概率密度較高的區(qū)域運(yùn)動(dòng)，并將最新的觀測(cè)值引入采樣過(guò)程，定義適應(yīng)度函數(shù)為：

(7)

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

步驟2執(zhí)行PSO算法。

① 隨機(jī)初始化粒子的位置和速度。

② 初始化個(gè)體最優(yōu)XPbest(t)和全局最優(yōu)XGbest(t)。

③ 根據(jù)式(7)計(jì)算每個(gè)粒子的適應(yīng)度值f(xi)。

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

⑤ 根據(jù)式(5)和式(6)更新粒子每個(gè)粒子的狀態(tài)。

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

步驟3計(jì)算粒子的重要性權(quán)值并歸一化。

(8)

步驟4重采樣。

步驟6計(jì)算狀態(tài)估計(jì)值：

(9)

2 基于CUDA的PSOPF算法

2.1 CUDA并行計(jì)算架構(gòu)

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

圖1 CUDA編程模型

全局內(nèi)存的訪問(wèn)以32、64或128字節(jié)為單位，并且以32字節(jié)為基準(zhǔn)對(duì)齊。當(dāng)線程束中所有線程訪問(wèn)連續(xù)且對(duì)齊的內(nèi)存空間時(shí)，則所有訪問(wèn)可以合并為一次內(nèi)存讀取。如果內(nèi)存地址不在同一內(nèi)存區(qū)段中，則無(wú)法合并訪問(wèn)，內(nèi)存訪問(wèn)性能急劇下降。并且，當(dāng)同一線程束中的內(nèi)存訪問(wèn)不連續(xù)時(shí)，盡管線程不需要這些訪問(wèn)之間的數(shù)據(jù)，但這些數(shù)據(jù)仍然會(huì)被傳輸，浪費(fèi)了內(nèi)存帶寬。因此，為了獲得好的數(shù)據(jù)吞吐量，需要通過(guò)合并訪問(wèn)將大量的內(nèi)存獲取請(qǐng)求合并，減少內(nèi)存獲取次數(shù)。

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ī)實(shí)數(shù)；j是服從1到N之間均勻分布的隨機(jī)整數(shù)；j是在k時(shí)刻所選粒子的索引，該粒子是在k+1時(shí)刻的粒子i復(fù)制的祖先。

拒絕重采樣可以作為單個(gè)CUDA內(nèi)核函數(shù)調(diào)用實(shí)現(xiàn)，使每一個(gè)線程對(duì)應(yīng)一個(gè)粒子，線程為其代表的粒子選取和復(fù)制祖先。拒絕重采樣在隨機(jī)抽選粒子時(shí)，存在非合并的全局內(nèi)存訪問(wèn)。對(duì)于較大的粒子數(shù)，非合并訪問(wèn)次數(shù)增多，導(dǎo)致算法執(zhí)行效率降低。考慮到線程束中線程在同一內(nèi)存區(qū)段上執(zhí)行讀/寫(xiě)操作能夠避免非合并的內(nèi)存訪問(wèn)，提出改進(jìn)的拒絕重采樣算法，其偽代碼如下：

算法2改進(jìn)的拒絕重采樣算法

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為全局內(nèi)存中一組固定數(shù)量的連續(xù)區(qū)段，最小的s為單個(gè)內(nèi)存區(qū)段，最大的區(qū)段為跨越權(quán)值數(shù)組的所有區(qū)段的集合。設(shè)Esize是一個(gè)權(quán)值所占字節(jié)大小，Ssize為一個(gè)s區(qū)段所占字節(jié)大小；Scount是s區(qū)段的數(shù)量，則Scount=Esize×N/Ssize。Ecount是一個(gè)s區(qū)段中權(quán)值的數(shù)量，則Ecount=Ssize/Esize。j是服從的s段中第一到最后一個(gè)元素索引之間均勻分布的隨機(jī)整數(shù)；s服從第一個(gè)到最后一個(gè)s區(qū)段索引之間均勻分布的隨機(jī)整數(shù)。為方便線程和全局內(nèi)存對(duì)齊，粒子數(shù)N的值設(shè)為2的次冪。

2.3 基于CUDA的PSOPF并行算法

通過(guò)對(duì)PSOPF算法的并行性分析，采取線程與粒子一一對(duì)應(yīng)的思路，實(shí)現(xiàn)基于CUDA的PSOPF并行算法。具體實(shí)施思路為：創(chuàng)建與粒子數(shù)相同的線程數(shù)，并分配相應(yīng)的存儲(chǔ)空間，每個(gè)線程執(zhí)行一個(gè)粒子的運(yùn)動(dòng)計(jì)算。PSOPF算法流程中更新全局最優(yōu)和計(jì)算狀態(tài)估計(jì)值會(huì)產(chǎn)生粒子間的數(shù)據(jù)關(guān)聯(lián)，無(wú)法完全并行化。本文采用并行規(guī)約算法尋找極值和求和運(yùn)算，每個(gè)線程處理一個(gè)粒子，提高并行度，以充分發(fā)揮GPU并行計(jì)算性能。根據(jù)PSOPF算法的計(jì)算流程，制定的基于CUDA的PSOPF并行算法流程如圖2所示。

圖2 基于CUDA架構(gòu)的PSOPF并行算法流程

3 實(shí)驗(yàn)結(jié)果及分析

3.1 算法性能測(cè)試

為驗(yàn)證本文算法的精度和運(yùn)行時(shí)間等性能，采用一維非線性系統(tǒng)模型在CPU和GPU平臺(tái)上進(jìn)行實(shí)驗(yàn)仿真，并與PF以及PSOPF算法進(jìn)行對(duì)比。實(shí)驗(yàn)所采用的計(jì)算平臺(tái)如表1所示。

表1 計(jì)算平臺(tái)

狀態(tài)模型和量測(cè)模型如下：

(10)

(11)

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

(12)

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

圖3 不同算法的狀態(tài)估計(jì)結(jié)果

圖4 不同算法的RMSE對(duì)比結(jié)果

表2 不同算法RMSE對(duì)比

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

在不同粒子數(shù)時(shí)本文算法的運(yùn)行時(shí)間如表3所示。

表3 本文算法運(yùn)行時(shí)間對(duì)比

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

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

圖5 不同算法的加速比對(duì)比

由圖5可見(jiàn)，隨著粒子數(shù)從28增加到214，本文算法的加速比增長(zhǎng)快于PSOPF算法，這是由于隨著粒子數(shù)的增加，算法計(jì)算量增加，采用改進(jìn)拒絕重采樣的本文算法避免了非合并訪問(wèn)，提高了GPU內(nèi)存使用效率，顯著地提高了算法的實(shí)時(shí)性。

3.2 變槳距系統(tǒng)的故障診斷仿真與分析

風(fēng)電機(jī)組變槳距系統(tǒng)由三個(gè)相同且獨(dú)立的槳距控制器組成，每個(gè)槳距角控制器被建模為實(shí)際的槳距角β與其參考值βref間的閉環(huán)傳遞函數(shù)。βref為閉環(huán)傳遞函數(shù)的輸入，由系統(tǒng)的槳距角控制器提供；β為傳遞函數(shù)的輸出，作為槳距角的測(cè)量值。實(shí)質(zhì)上，變槳距系統(tǒng)是一個(gè)活塞伺服系統(tǒng)，可以通過(guò)二階傳遞函數(shù)建模[17]：

(13)

(14)

表4 故障模式表

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

(15)

圖6 變槳距系統(tǒng)故障診斷原理框圖

首先進(jìn)行故障檢測(cè)，設(shè)置殘差閾值δ(δ>0)，當(dāng)殘差平滑值d(0)大于閾值δ時(shí)，則判斷系統(tǒng)發(fā)生故障，否則判斷系統(tǒng)無(wú)故障。檢測(cè)到故障發(fā)生后，通過(guò)比較每個(gè)故障模式的殘差平滑值d(m)(m=1,2,3)的大小來(lái)進(jìn)行故障隔離，若其中某一個(gè)模式的殘差平滑值最小，則判斷系統(tǒng)發(fā)生了該模式故障。

根據(jù)正常模式的系統(tǒng)模型計(jì)算當(dāng)前時(shí)刻的殘差平滑值d(0)，當(dāng)d(0)大于閾值時(shí)，則檢測(cè)到系統(tǒng)發(fā)生故障。系統(tǒng)正常運(yùn)行和發(fā)生故障M1-M3時(shí)的殘差曲線如圖7所示。

(a) 系統(tǒng)正常運(yùn)行時(shí)的槳距角殘差

(b) 系統(tǒng)發(fā)生故障M1時(shí)的檢測(cè)結(jié)果

(c) 系統(tǒng)發(fā)生故障M2時(shí)的檢測(cè)結(jié)果

(d) 系統(tǒng)發(fā)生故障M3時(shí)的檢測(cè)結(jié)果圖7 不同故障的診斷結(jié)果

從圖7(a)可見(jiàn)，當(dāng)系統(tǒng)正常運(yùn)行時(shí)，槳距角殘差基本為零。從圖7(b)-(d)可以看出，對(duì)于故障M1-M3，系統(tǒng)無(wú)故障時(shí)的統(tǒng)槳距角殘差幅值較小，且變化比較平穩(wěn)；當(dāng)故障發(fā)生后，殘差變化劇烈且幅度高于閾值，因此判斷故障發(fā)生。

檢測(cè)到故障后，啟動(dòng)故障M1-M3模型進(jìn)行故障隔離，結(jié)果如圖8所示。

(a) 系統(tǒng)發(fā)生故障M1時(shí)的隔離結(jié)果

(b) 系統(tǒng)發(fā)生故障M2時(shí)的隔離結(jié)果

(c) 系統(tǒng)發(fā)生故障M3時(shí)的隔離結(jié)果圖8 不同故障的隔離結(jié)果

由圖8(a)可見(jiàn)，當(dāng)檢測(cè)到故障M1發(fā)生后，由于M3為槳距角傳感器偏差故障，起始?xì)埐钤?.3左右。M1和M2的殘差曲線較為相似，但M1的殘差是最小的，故判斷為M1故障。同理分析圖8(b)和(c)可見(jiàn)，當(dāng)系統(tǒng)實(shí)際狀態(tài)與某故障模型匹配時(shí)，其殘差平滑值相比其他故障模型較小，則判斷該故障模式發(fā)生。

4 結(jié) 語(yǔ)

本文基于CUDA架構(gòu)設(shè)計(jì)實(shí)現(xiàn)了一種GPU加速的PSOPF并行算法，通過(guò)分析PSOPF算法的并行性，采取線程和粒子一一對(duì)應(yīng)的方式實(shí)現(xiàn)并行化算法，運(yùn)行速度得到顯著提升。利用改進(jìn)的拒絕重采樣并行算法，解決全局內(nèi)存讀取時(shí)出現(xiàn)的非合并訪問(wèn)問(wèn)題，提高了PSOPF并行算法在GPU上的執(zhí)行效率，并將其應(yīng)用到風(fēng)電機(jī)組變槳距系統(tǒng)故障診斷中。實(shí)驗(yàn)結(jié)果表明，所提出的改進(jìn)的PSOPF并行算法有效地提高了實(shí)時(shí)性，是一種進(jìn)行故障診斷有效的方法。

下一步工作將設(shè)計(jì)并實(shí)現(xiàn)基于CUDA的多模型粒子濾波故障診斷方法，通過(guò)故障診斷的并行化，以提高故障診斷的實(shí)時(shí)性能。