基于聯合概率分布和置信邊界模型的風力發電機功率曲線自適應修正算法

徐姣新 楊 召

1(商丘學院 河南 商丘 476000)2(商丘工學院 河南 商丘 476000)

0 引 言

近十年來，風電裝機容量快速增長，風電裝機規模在世界能源結構中發揮著越來越重要的作用[1]。風能利用的質量也隨之受到了風電市場的重視。隨著風電場SCADA(Supervisory Control And Data Acquisition)系統歷史運行數據的積累，海量的數據源已成為提高風電經濟性可行的解決方案[2]。

由于風電機組功率曲線(Wind Turbine Power Curve，WTPC)經常隨時間變化，將實測輸出功率與風速聯系起來的數據驅動型WTPC已應用于較多領域。它是電力預測、性能評估和狀態監測等不同場景下的一種有效的技術手段，受到了極大關注[3-4]。WTPC建模主要有參數化建模和非參數化建模兩種技術路線。前者包括線性分段模型、多項式冪曲線、最大值法、Logistic函數和概率模型等，后者包括Copula冪曲線、三次樣條插值、神經網絡、模糊邏輯方法和數據挖掘算法等。但是，基于數據的風力發電應用程序的性能在很大程度上取決于數據質量[5-6]以及WTPC建模[7-16]SCADA系統可能存在傳感器故障、數據傳輸噪聲和損耗、有限的功率輸出或設備異常等原因，由此產生大量的異常數據和缺失數據。所以通常需要消除異常數據，恢復丟失數據。文獻[7]將加權歐氏距離定義為相似性度量，并采用局部離群因子算法來識別離群值，但若含有高比率無效數據的數據源，這種方法可能不準確。文獻[8]簡單地用物理極限作為識別邊界。文獻[9]則利用控制圖給出了過濾高度依賴于標準非線性WTPC的異常值的上/下控制限。文獻[10]對標準WTPC的左/右運動進行了優化，得到了上下限，但其中的PDL指數對優化不夠合理。

在假設每個風速區間輸出功率的概率分布為正態分布的情況下，可以利用平均值附近的多個標準差來確定WTPC的置信邊界，以消除異常值[11]。然而，正態分布的假設并不總是適用的。文獻[12]采用四分位算法、文獻[13]使用基于密度的噪聲應用空間聚類法(Density-Based Spatial Clustering of Applications with Noise，DBSCAN)給出異常值的識別邊界，而它們的性能對初始設定參數敏感。鑒于文獻[11-13]中基于每個風速區間邊界集合的方法總是忽略輸出功率對風速的總體分布特征，文獻[14]使用Copula工具將WTPC作為二元聯合分布。文獻[15]則遵循這一思想，將概率WTPC應用于風電機組狀態監測，采用經驗Copula方法。在其基礎上，文獻[16]采用GMCM(Gaussian mixed Copula)模型建立風速和輸出功率的聯合概率分布，并給出識別異常值的置信度曲線。上述Copula方法考慮了風電機組整個運行范圍內的整體聯合分布特性。然而，自然風速和產生的輸出功率通常對它們的聯合分布特性有不同的影響。此外，在風力機的分區運行區域中，多個Copula邊界模型的集合尚未得到充分的研究。同時，不同區域關節分布的尾部特征是不規則且多樣的。

特別是SCADA系統中的記錄也反映出棄風現象較為嚴重。如果消除了由棄風引起的堆積異常值，則不可避免地存在大量缺失數據。它們通常可以通過基于時間序列回歸、插值和統計方法的時間或空間方法來恢復[17]。然而，回歸方法通常需要大量連續的數據點。當連續缺失的數據點達到一定數量時，牛頓、三次樣條等插值方法的累積插值誤差往往會增大，此外，考慮到數據量有限，統計方法可能還不夠。因此，需要探索一種有效的插值方法。

本文在三維經驗Copula空間中，以轉子速度作為輔助辨識因子進行初步數據過濾，并對WTPC的離群點識別邊界進行了初步改進，提出一種基于混合阿基米德Copula函數的自適應置信邊界模型。整個邊界是變速變槳控制(Variable Speed Variable Pitch，VSVP)風力機分區運行區域邊界的集合。在每個操作區域中，利用梯度下降進行期望最大化(Expectation Maximization，EM)。然后，將推導出的Copula條件概率與各風速下的條件核密度估計(Conditional Kernel Density Estimation，CKDE)相結合，建立了置信邊界模型，并提出保證性能的自適應建模評價系統。在邊界模型剔除異常點后，提出一種雙向馬爾可夫鏈插值(BMCI)方法，用優化的前向和后向權值恢復缺失的數據點。因此，所提方法有助于獲得更準確的WTPC用于監測、故障檢測和功率預測。

1 基于置信邊界的異常值剔除

1.1 初步清理

通常，具有雙饋感應發電機(Doubly Fed Induction Generator，DFIG)的主流三葉片水平軸風力發電機具有VSVP能力，其整個運行范圍可分為五個區域，如圖1所示。風速(V)、轉子轉速(ωrot)和輸出功率(P)之間的關系能有效地反映各區域的主要運行特征。帶有“rated”和“min”的變量分別表示相關變量的額定值和最小值；Vcut-in和Vcut-out分別表示切入風速和切出風速。

圖1 風機三維運行區域

在I區，DFIG空載空轉，不并網，不等于零的輸出功率值需要清除。

在II區，DFIG通過雙向PWM變換器連接電網。通過核密度估計(Kernel Density Estimation，KDE)計算風速、轉速和輸出功率的累積邊際概率分布，并將其范圍劃分為若干個小的區間。采用經驗Copula方法建立聯合概率分布。然后，在三維Copula空間中，聯合概率點分布在對稱線交叉點(0，0，0)和(1，1，1)周圍。點與線之間的距離表示點的偏差程度。同時，最偏離點的概率往往較低。為了消除一定置信水平下的異常值，在對稱直線上設置距離閾值，對偏差程度明顯的點進行隔離。然后，在一定的置信水平下，優先清除其中概率較低的點。在Copula空間中，一個bin對應于許多實際的數據點。考慮到小區間劃分間隔對建立經驗Copula模型的影響，利用PDL索引來最大化剩余數據量，并保證滿足置信度要求。

在III區和IV區中，可以執行相同的過程來清除原始數據。當然，應該考慮轉子轉速和輸出功率的最大物理極限。

在V區，風輪葉片通常變成全順槳，渦輪轉子以非常低的速度旋轉或剎車。同時，電氣系統不并網，輸出功率保持恒定為零。

1.2 WTPC的置信邊界建模

經過初步的數據清理，風速和輸出功率的散點圖主要表現為II、III、IV區的帶狀分布，這些區域的控制策略和算法不同，因此最好能準確地剔除不同區域的異常值。

首先，估計風速和輸出功率的邊際概率分布。隨后，需要在每個區域選擇合適的Copula函數。正態Copula、t-Copula和Frank-Copula等具有對稱尾特征的Copula函數不能捕捉隨機變量的非對稱尾特征。考慮到Gumbel Copula和Clayton Copula能夠捕捉隨機變量的非對稱尾特征，因此采用了Gumbel Copula、Clayton Copula和Frank Copula等混合函數。

具體的建模步驟如下。

步驟1采用KDE方法分別估計II、III、IV區的累積邊際分布FP(P)和FV(V)。

步驟2混合阿基米德Copula函數如下：

C(u,v)=φGCG(u,v;θG)+φCCC(u,v;θC)+

φFCF(u,v;θv)

(1)

式中：u=FP(P)和v=FV(V)；φ*是加權系數；θ*是隨機變量間的關聯關系；C*(·)是Copula函數。此外，Gumbel函數、Clayton函數和Frank Copula函數如下：

CG(u,v;θG)=exp{-[(-lnu)θG+(-lnv)θG]1/θG}

(2)

CC(u,v;θC)=max[(u-θC+v-θC-1)1/θC,0]

(3)

(4)

它們分別對應于聯合概率分布的UT(Upper Tail)特征、LT(Lower Tail)特征和ST(Symmetric Tail)特征。

在FP(P)和FV(V)的基礎上，構建聯合頻率直方圖，并擬合各操作區域的混合阿基米德Copula函數。采用梯度下降優化的EM估計式(1)未知參數。

步驟3利用混合Copula函數C(u,V)可以建立聯合概率分布C(FP(P)，FV(V))。然后，得到FV(V)下的FP(P)的條件概率為：

(5)

在一定條件下，將V和P的聯合概率分布與局部概率分布聯系起來。然后，在每一個V值下，計算其條件概率。

步驟4對于II、III和IV之間的轉換，由于不確定性，很容易產生異常數據。如果設置了各區域FP(P)的置信水平，則異常數據量可以降低到一定的水平。FP(P)在一個操作區域內具有置信度1-β。考慮到FP(P)的不對稱性，將不對稱系數設為k。那么，置信度1-β的概率分位數可以計算如下：

βlow=κβ

(6)

βup=1-(1-κ)β

(7)

式中：βup和βlow是上概率分位數和下概率分位數。

步驟5對于風速值V=V0，FV(V0)可以通過其KDE得出。通過式(5)，F(FP(P)|FV(V0))成為單變量概率分布。對應于βup和βlow，F(FP(P)|FV(V0))的輸出分別為γup和γlow，這是P在值V0下的條件概率。實際上，它們相當于條件概率分布H(P|V0)的輸出。那么，V0下P的分位數如下：

Pup=H-1(γup|V0)

(8)

Plow=H-1(γlow|V0)

(9)

式中：Pup和Plow是在一定值V0下WTPC的上下置信邊界。

FP(P)和H(P|V)通過式(5)聯系起來。然后綜合考慮聯合概率分布和局部概率分布，得到期望置信水平。

步驟6在Ⅱ、Ⅲ和Ⅳ區，按一定間隔計算各風速值的上下邊界。然后，對所有區域的點進行描述，得到WTPC的等價置信邊界。利用上下邊界，可以直接識別和消除異常數據。

利用上述方法，可以對WTPC的置信邊界進行建模。然而，它是基于一定數量的SCADA歷史數據的統計模型。直觀地說，風力發電機的動力會隨著時間而變化。因此，在評估數據清理效果時，需要更新置信邊界模型。

1.3 模型評估與自適應更新

利用WTPC的上述置信邊界模型，直接從SCADA系統中清除原始數據。比較原始數據和清理數據的概率分布，清理數據的概率分布變得更加規則。為了量化這些變化，本文給出了幾種統計指標來表示清理數據在概率分布上的變化趨勢。

對于清潔能源輸出數據的分布，定義其顯著性水平為α。使用式(6)和式(7)確定置信度1-α下的上下分位數。然后，指數置信帶寬比(Confidence Bandwidth Ratio，CBR)可以由式(10)計算。

(10)

式中：χCBR代表CBR指數；ΔP是上分位數和下分位數之間的帶寬。另外，其他兩個參數定義如下：

(11)

(12)

式中：μSke是偏度指數；δKur是峰度指數；D2、D3和D4分別是2階、3階和4階中心距。

CBR指數表示有效數據的平均概率分布水平，較大的CBR意味著在一定的置信度下，更多的輸出功率數據分布在相對較少的帶寬上。偏度指數表明對稱性，如果它接近于零，則概率分布更為對稱。峰度指數反映了其峰值附近概率分布的集中程度。峰度指數越大意味著數據在峰值附近越集中。

基于以上評價指標，可以看出數據清洗的變化趨勢。然而，它們都是描述清理數據的概率特性的統計指標。在一定程度上，它們受數據樣本大小的影響。凈化數據的評價指標越穩定，反映了置信邊界WTPC模型對數據凈化的穩定效果。本文采用隨機抽樣和交叉驗證(Cross Validation，k-CV)方法進行研究，從而獲得一個穩定且普遍適用的模型

對于II、III或IV區，需要選擇適當的時間尺度以獲得足夠的原始數據樣本量。然后，在對原始數據進行初步數據清理后，進行k-CV法。對于數據樣本，它們被隨機分成k個部分。每個部分都用作測試數據集，其余的k-1部分用作訓練數據集。然后，對數據樣本進行隨機劃分，進行k次訓練和測試，得到k個CBR指數、k個偏態指數和k個峰度指數。如果它們是收斂的，則使用k組的平均值來評估數據清理效果。為了保證統計指標的穩定性，可以將數據樣本的隨機劃分和非參數模型的k-CV重復n次以驗證其性能。

如果n組平均指標在常數附近收斂到要求的精度，則認為WTPC的置信邊界模型是可靠的。由于風電機組運行特性隨時間的變化，需要對SCADA數據以及數據過濾模型進行更新。本文采用基于滑動時間窗(Sliding Time Window，STW)機制的時間驅動模型更新方法。合理選擇時間窗T，對采集到的數據進行增量更新，并保留足夠的數據樣本，以獲得穩定的模型，從而構建了WTPC置信邊界模型的自適應更新流程。

2 雙向馬爾可夫鏈插值

處理數據中不可避免地存在數據丟失。如果連續丟失大量數據點，即使重復執行數據恢復，也會造成很大的累積誤差。此外，剔除異常值后，剩余數據量也會減少。考慮到馬爾可夫鏈無后遺癥，基于隨機切換統計的馬爾可夫鏈可以用于不連續數據，且需要較少的數據量。此外，馬爾可夫鏈還可以利用時間序列的時變特性。

為了充分利用剔除異常值后的剩余數據量，本文提出一種前向和后向馬爾可夫鏈相結合的雙向機制。首先，在一定的區間內對輸出功率范圍進行均勻劃分，得到連續的離散狀態。然后，選擇缺失數據段前后的剩余數據樣本。在此基礎上，分別建立了前向馬爾可夫鏈和后向馬爾可夫鏈。前者表示為：

…,Pt-Np)

(13)

式中：Pt是時間t時輸出功率的離散狀態；0，1，…，t是正向時間序列號；NP是馬爾可夫鏈的階。雖然風電隨時間的變化而變化，但在變化中存在一定的連續性，從而提供了變化的向后可追溯性。然后，后向馬爾可夫鏈如下：

…,Pt+Np)

(14)

式中：inf，1，…，t是后向時間序列號。使用式(15)和式(16)，在t時恢復的離散狀態分別為前Pt和后Pt。然后，使用隨機數生成器如下：

(15)

(16)

式中：Pi是P的第i個狀態；“upp”和“low”表示P的狀態劃分區間的上下限；εi是均勻分布在[0,1]中的隨機數；ω1和ω2是加權系數。

為了提高丟失數據恢復的適應性，采用均方根誤差(Root Mean Square Error，RMSE)指標和梯度下降算法對加權系數進行優化。

在剔除異常值和恢復丟失數據后，對SCADA系統的運行數據進行了充分的預處理，得到的數據對許多應用都是有效的。

3 算例分析

3.1 異常值過濾

為了驗證本文方法的有效性，本節對其進行了仿真驗證和應用分析。原始數據來自華北某風電場的SCADA系統。時間為2017年1月1日至12月31日，采樣周期為10 min。選擇了主流的1.5 MW DFIG。它可以在圖1所示的模式下工作。切入、切出、額定轉速和轉子功率的風速分別為3 m/s、25 m/s、8.87 m/s和10.8 m/s。最小、額定和最大風力渦輪機轉子轉速分別為9 r/min、17.3 r/min和18 r/min。額定功率為1.5 MW。額定轉速下的有功功率為0.99 MW。

使用k-CV過程檢查不同操作區域中數據樣本的充分性和一致性。通過測試，三個月左右的時間窗可以提供足夠的數據量建立一個可靠的WTPC置信邊界模型。當然，如果不能滿足要求，可以延長時間窗口，增加數據量。選擇2017年WT1的1月-3月和WT2的7月-9月這兩個時間段，顯示完整的驗證過程。如表1所示，隨機分割重復5次，3倍交叉驗證的測試數據集的平均指標均穩定。結果表明，三個月的時間窗可以得到一個穩定的WTPC置信邊界模型。利用時間窗T=3個月的數據，將WTPC的置信邊界建模過程顯示如下。

表1 評價指標驗證過程

續表1

1) 從SCADA系統采集原始數據，用本文方法對數據進行初步清理，最終劃分區間到邊際概率為0.023。WT1和WT2的PDL指數分別為90.3%和87.6%。Copula空間中的聯合概率分布如圖2所示，其中顯著偏離的異常值具有較深的顏色且概率較小。相對于風速，轉速對輸出功率的不確定性影響較小。垂直于直線，計算從對稱直線到(0,0,0)和(1,1,1)點的距離。

(a) WT1：1月-3月 (b) WT2:7月-9月圖2 Copula空間中的聯合概率分布

Copula空間中距離的概率密度如圖3所示。在0附近的第一個峰值附近，風速和轉子速度的影響相互重疊。由于轉子速度的最大點在這里，所以得到了轉子速度的最大概率密度。然后，隨著距離的增加，受轉速影響的概率密度不斷減小到谷點。在第二個峰值，距離的概率密度主要是由風速引起的。因此，選擇第一個峰值處的距離作為沿轉子速度維度的距離閾值。選擇第二個峰值處的距離作為沿風速維度的距離閾值。它們被用來識別顯著偏離的異常值，并避免僅僅通過置信度消除錯誤的異常值。在此基礎上，將置信度設置為0.95。過濾后的數據量對置信水平敏感，應保證有更多的剩余數據量。原始數據和初步清理數據的散點圖如圖4所示。這表明顯著的異常值已經被清除，剩余的數據變得更加集中。

(a) WT1：1月-3月 (b) WT2:7月-9月圖3 Copula空間距離概率密度

(a) WT1：1月-3月 (b) WT2:7月-9月圖4 數據清理散點圖

2) 根據初步清理的數據，用KDE計算風速和輸出功率的累積邊際概率分布。然后，用梯度下降優化的EM算法估計混合阿基米德Copula函數的權值和參數，如表2所示。

表2 Copula函數估計及1-β置信下概率分位數值

3) 推導出式(5)F(FP(P)|FV(V))中的條件概率分布。在每個區域，將1-β設置為0.95。在已知值V0下，用F(FP(P)|FV(V0))分別計算γup和γlow，如圖5所示。結果表明，不同風速下，輸出功率的不確定性隨湍流強度的變化而變化。它符合實際輸出功率的特點。

(a) WT1：1月-3月

(b) WT2:7月-9月圖5 不同風速下的γup和γlow

4) 用Parzen-Rosenblatt-CKDE計算了各區域的條件概率分布H(P|V)。對應于V0下的每一對γ，反求H(P|V0)，可得到分位數Pup和Plow。計算每個風速值對應的Pup和Plow，并描繪所有區域的所有邊界點。然后，可以得到WTPC的等效置信邊界，如圖6所示。利用WTPC的上下置信邊界，可以直接識別和消除異常數據。

(a) WT1：1月-3月 (b) WT2:7月-9月圖6 WTPC的散點以及置信邊界

5) 為了顯示最終清理數據的效果，比較不同階段的輸出功率fP(P)的概率密度曲線，如圖7所示。評價指標比較見表3。結果表明，最終清洗后的數據更加集中，置信邊界模型在各個區域都具有可靠的數據清洗性能。

(a) II區兩時段

(b) III區兩時段

(c) IV區兩時段圖7 不同區域概率分布對比

表3 不同階段數據評價指標

此外，為了進一步證明本文方法的有效性，還將其與經驗Copula和GMCM方法進行了比較。使用經驗Copula方法，Copula空間和散點圖中的聯合概率分布如圖8所示。圖8(b)和圖8(d)中的輪廓是沿著剩余數據點的邊緣描繪的。圖8(b)和圖8(d)中等高線外的數據點是消除的異常值。結果表明，在整個工作區域內，輪廓線非常不規則。

(a) WT1聯合概率分布 (b) WT1散點圖

(c) WT2聯合概率分布 (d) WT2散點圖圖8 聯合概率分布及功率散點圖

不同區域聯合概率密度的GMCM擬合和散點圖如圖9所示。不同區域的概率密度差異很大，如圖9(a)和圖9(c)所示，得到的輪廓也非常不規則，如圖9(b)和圖9(d)所示。結果表明，僅用概率密度剔除異常值，不能平衡WTPC的整體和局部特征。它們不足以消除異常值。

(a) WT1不同地區的GMCM (b) WT1功率散點圖

(c) WT2不同地區的GMCM (d) WT2功率散點圖圖9 不同地區的GMCM及功率散點圖

(17)

使用式(17)中的R2指數比較不同方法的數據清理效果。證明了Copula空間中從直線到(0,0)和(1,1)數據概率分布的中心性。R2值越大，數據的概率分布越集中。這也意味著數據清理效果更好。計算結果見表4。

表4 不同方法的評價指標

本文方法除了在II區中有相似的性能外，在III區和IV區中也有較好的性能。此外，在本文的初步數據清理階段，風速、轉子轉速和輸出功率的三維空間具有較強的識別和剔除異常值的能力。然而，這一優勢并沒有得到R2指標的正確評價，需要定義一個更適合三維空間的指標。

綜上所述，本文提出的異常值剔除方法包括初步數據清理和置信邊界建模是非常有效的。與經驗Copula和GMCM方法相比，該方法在平衡WTPC上運行數據的整體和局部分布特性的同時，具有更好的性能。

3.2 丟失數據恢復

本文利用二階馬爾可夫鏈建立輸出功率模型。優化后的權重ω1和ω2分別為0.609 8和0.390 2。當連續缺失10個數據點時，比較BMCI法與牛頓法和分段三次Hermite插值法，結果如圖10和表5所示。隨著連續缺失點數的增加，BMCI方法的精度不斷提高。

圖10 不同插值方法對比

表5 不同數據恢復方法的RMSE

3.3 基于處理數據的WTPC建模應用

處理后的數據被用于生成代表風機實際運行特性的WTPC。為了更好地刻畫風速與輸出功率之間的非線性關系，采用了4參數和5參數Logistic函數等性能較好的參數化方法對處理后的數據進行擬合。此外，還采用傅里葉函數進行比較。曲線擬合函數表示為：

(18)

式中：PTheo是WTPC的功率；a*、b*、c*、m*、n*、g*和τ*是擬合參數。

擬合曲線如圖11所示。結果表明，所安裝的WTPC與制造商提供的WTPC完全不同。基于處理數據的WTPC位于置信區間之間。另外，雖然四參數和五參數logistic函數的參數較多，但三種方法的擬合性能接近。它反映了使用最終清理數據的穩定擬合性能。為了驗證數據清理程序的有效性，使用擬合的風電機組來計算風電機組的理論功率。使用以下指標評估計算出的理論功率的性能。

(a) WT1：1月-3月 (b) WT2:7月-9月圖11 不同方法的WTPC比較

(19)

理論和實際功率曲線如圖12所示。結果表明，擬合后的WTPC理論功率曲線在上下置信邊界之間，有較好的位置。評價指標見表6。結果表明，與制造商的計算結果相比，所擬合的WTPC的計算結果具有更大的相關性和更小的誤差。

(a) WT1：1月-3月

(b) WT2:7月-9月圖12 不同方法功率對比

表6 評價指標比較

4 結 語

本文利用轉子速度的狀態變量作為輔助辨識因子，在三維Copula空間中合理地消除顯著異常值，并基于混合阿基米德-Copula函數的EM估計，將各區域輸出功率的總體概率分布和各風速下輸出功率的條件概率相結合。總體而言，本文方法綜合考慮了運行數據的機理和概率分布。它不同于只考慮概率密度或偏差距離的方法。而且，所用的評估系統對于置信邊界模型的自適應更新也極有效。此外，所提出的BMCI方法能有效地填補連續缺失數據點，減少累積誤差。最后，算例證明了本文方法優越的性能，所用數據處理過程具有巨大的應用潛力。