計及日特性和季特性的風矢量仿真模型及應用

李 晗 鄒建華 張遠航 楊保華

(三峽大學 電氣與新能源學院，湖北 宜昌 443002)

風能是一種清潔的可再生能源，發電成本低且儲量巨大，因此風電規模在我國及全球范圍內發展迅速.2015年，我國新增風電裝機容量30.5 GW[1]，總容量和發電量分別達到145.1 GW和1.863×105GW·h[2].截至2021年2月，我國風電裝機總容量達到284.56 GW，其中2021年1～2月新增裝機容量達到3.57 GW[3].風速的隨機性和間歇性使得風力發電具有波動性和不可控性，隨著風電滲透率的升高，風電并網可能對電網的可靠性產生影響.研究含風能電力系統充裕度評估方法，評估風能對電力系統的可靠性效益，可以為電力系統運行與規劃提供幫助.

準確的風速模型是模擬風電場出力的重要條件，常用風速模型可分為概率分布模型和時間序列模型兩類.概率分布模型能夠較好地刻畫風速的分布特性[4-6]，但難以描述風速隨時間動態變化的特性；時間序列模型大多可分為兩類：一種基于回歸分析模型[7-9]，一種基于馬爾科夫鏈模型[10-11].風速具有復雜特征，既有一定的周期性，又同時隨季節更替而變化.為研究風速的日周期性，文獻[11]按時序將歷史風速數據分為24組數據集，分別建立24個馬爾科夫鏈模型，并且提出季節系數以描述風速的季節特性.除日周期特性外，風矢量也有一定的月周期特性，文獻[12]使用傅里葉變換將三角函數與國內多個地區歷史風速數據進行擬合，結果表明該地區風速具有明顯的月周期特性.當前多數模型專注于研究風速隨機性，較少刻畫風向隨機性，而風向是刻畫風電場尾流效應的重要變量.尾流效應對風電機組出力有顯著影響，導致風電機組發電量損失可達10%～15%[13].因此單一的風速模型難以計及尾流效應影響，進而導致對風電機組出力的樂觀預測.另一方面，風速和風向具有一定的季特性，文獻[14-15]分別統計國內和國外多個地區的歷史風速數據，結果顯示風速均值在不同季節存在較大差異.目前大多數風速模型并未考慮季特性，可能導致充裕度評估誤差.

為考慮尾流效應和風矢量的季特性，本文將風速和風向定義為風矢量，提出計及日特性和季特性的風矢量仿真模型.該模型的建模過程主要分為兩步：第一步，將實測數據分為24組數據集，分別建立馬爾科夫鏈模型得到計及日特性的風矢量模型，其中為劃分二維風矢量的馬爾科夫狀態空間，提出一種基于正交分解準則的聚類算法；第二步，引入最優季度系數，逐季度修正第一步得到的風矢量樣本，使其具有季節特性.另一方面，本文將提出的風矢量仿真模型與Jensen尾流模型和時序蒙特卡羅模擬法結合，提出計及日特性和季特性的含風能發電系統充裕度評估方法，分析風機布局和峰值負荷對系統充裕度的影響.

1 風矢量定義及其日特性和季特性

本文將風速和風向定義為二維笛卡爾坐標系中的方向性矢量：風矢量w=(v，θ)，如圖1所示.其中v為風速，表示風矢量長度；θ表示風向，為風矢量與x軸的夾角.定義正東方向為0°，并與x軸正方向重合，則當風矢量沿逆時針方向偏離x軸時，風向從0°增大至360°.此外，vx和vy分別為風矢量在x軸和y軸的分量.

圖1 風矢量示意圖

為說明風矢量的日特性和季特性，通過訪問開放氣象數據庫North Dakota Agricultural Weather Network(NDAWN)[16]，得到美國North Dakota州Berthold觀測站2006年至2010年間的實測每小時風矢量數據.圖2所示為24 h風速和風向均值與標準差.可以看到，在下午時段風速均值較大，同時風速和風向的波動性略高于平均水平.這表明由風速和風向組成的風矢量具有明顯的日特性，應以24 h為周期，研究風矢量的日特性.

圖2 各時刻實測風矢量柱狀圖

圖3為2006年至2010年間各季節實測風矢量均值和標準差.可以看到，在相同年份，四季風速風向均值和標準差不盡相同；在不同年份，同一季節的風速風向均值和標準差也存在差異.總的來說，風速的均值和標準差在春季和冬季較高，在夏季和秋季較低；風向的均值在秋季和冬季較高，標準差在春季最高.所以，由風速和風向組成的風矢量的季特性，體現在各季的風速風向均值和標準差隨著季度和年份的更替而波動變化.

圖3 各季節實測風矢量對比圖

2 計及日特性和季特性的風矢量仿真模型

2.1 計及日特性的風矢量馬爾科夫鏈模型

將每天同一時刻的歷史風矢量數據劃分至同組集合，則可得到24組風矢量數據集.以該24組風矢量數據集，建立24個馬爾科夫鏈模型，進而抽樣產生計及日特性的風矢量仿真樣本，以下為詳細建模過程.建立馬爾科夫模型，首先要將風矢量劃分為多個馬爾科夫狀態空間.由于風矢量中的風向為圓周變量，對風矢量狀態劃分造成困難，故本文采用正交分解方法將風矢量轉換為圖1二維笛卡爾坐標系中的x、y軸線性分量vx和vy，而后通過聚類算法進行風矢量狀態的比較和劃分.

以第1小時為例，該小時的風速馬爾科夫鏈的建模步驟如下：

步驟1：將第一時刻所有風矢量數據轉換為x，y分量形式：vx=vcosθ，vy=vsinθ；

(1)

步驟4：應用每組風矢量聚類中的風矢量x軸和y軸分量的均值更新該組風矢量聚類中心，進而得到5組更新后的風矢量聚類中心；

步驟5：重復步驟3和4，不斷更新風矢量聚類中心，直到當前聚類中心與前一聚類中心的歐式距離低于δ(δ=0.001)為止，得到第1小時中5組風矢量聚類中心和各聚類中的風矢量數據.

其他時刻的風矢量數據聚類算法和第1時刻相同，至此完成對一天24個時刻的風矢量聚類劃分.

根據上述步驟的劃分結果，可將每時刻風矢量均劃分為5組聚類空間.對每時刻的風矢量建立馬爾科夫鏈模型，以每組聚類空間作為馬爾科夫狀態空間，建立每個時刻的風矢量所屬聚類空間的轉移率矩陣，見式(2)：

(2)

式中，轉移率矩陣G的每個元素表示一種轉移關系，行數和列數分別表示轉移前后的狀態.例如，λ[1，1]和λ[1，5]分別表示風矢量狀態S1向自身轉移的轉移率和向狀態S5轉移的轉移率.

根據馬爾科夫鏈方法，轉移概率的表達式如下：

(3)

式中：fi，j為風矢量從第i個狀態轉移到第j個狀態的轉移次數；K為總狀態數.

根據前面的每時刻風矢量聚類狀態劃分結果，需對每個風矢量聚類狀態建立聯合概率分布模型，以抽樣每個馬爾科夫狀態空間的風矢量具體數值.本文應用copula方法建立風速與風向的聯合概率分布，以此表示風矢量概率分布.根據copula方法，風矢量的概率密度函數可表示為[17]：

f(v，θ)=c(F1(v)，F2(θ))f1(v)f2(θ)

(4)

式中：c(v，θ)為copula概率密度函數；f1(v)和f2(θ)分別為風速和風向邊緣概率密度函數；F1(v)和F2(θ)分別為風速和風向邊緣累積概率密度函數.copula函數種類較多，具體copula函數類型的選擇過程在4.1節中敘述.

本文選擇風速模型中應用最為廣泛的Weibull模型為風速變量的邊緣概率分布模型.Weibull分布的概率密度分布和累積概率分布參見文獻[18].

采用Von Mises模型為風向的邊緣概率分布模型，Von Mises模型可描述隨機變量的多峰分布特征，其分布模型參見文獻[19].

風速和風向邊緣分布中的參數均通過最大似然估計法求出.至此，建立風矢量馬爾科夫鏈模型和每個馬爾科夫狀態內的風速風向聯合分布模型.通過先抽樣馬爾科夫狀態空間序列，再抽樣風速風向的累積分布值，最后逆抽樣得到風速風向的具體數值，流程如圖4所示.

圖4 計及日特性的風矢量模型建立和抽樣流程圖

2.2 計及季特性的最優季度系數

w+=(ωvv*，ωθθ*)

(5)

式中：ωv為風速修正系數；ωθ為風向修正系數；w+為修正后的風矢量；v*和θ*為計及日特性的風速樣本和風向樣本.

以風速修正為例，根據樣本均值和標準差計算公式，可知修正前后風速均值和標準差滿足以下條件：

(6)

(7)

根據式(6)和(7)，修正后風速數據的第t年第j季度的均值和標準差與實測風速均值和標準差的偏差可表示為以下偏差函數：

(8)

式中：ωv，μ和ωv，σ為風速的均值和標準差在偏差函數中的權重系數，可根據工程人員的實際需要進行設定和調整.該式表明，修正風速的各季度的風速均值和標準差與實測風速均值和標準差間的偏差可表示為關于修正系數ωv的非線性函數.

對式(8)中偏差函數fv(ωv)求取關于修正系數ωv的偏導函數hv(ωv)，其表達式如下：

(9)

根據經典優化算法可知，當選取最優修正系數ωv使式(9)中偏導數為零，可使得式(8)偏差函數達到最小值.即令式(9)為零，可得到最優修正系數的表達式如下：

(10)

類似地，可推導得到第t年第j季度風向修正系數為：

(11)

式中：ωθ，μ和ωθ，σ分別為風向的均值和標準差在風向修正過程中的權重系數.

圖5為引入風矢量季特性的流程圖，根據式(10)和(11)，由每季度計及日特性的模擬風速樣本和實測風速樣本的均值和標準差，可計算得到每季度的最優季度系數，并使用該系數使修正后的風矢量具備季特性.

圖5 風矢量季特性引入流程圖

2.3 計及日特性和季特性的風矢量仿真算法

根據2.1節和2.2節內容，可給出計及日特性和季特性的風矢量仿真算法，以一個季度的抽樣過程為例，其詳細步驟如下.

步驟1：按照時刻劃分，分別抽樣24個馬爾科夫狀態空間序列；

步驟2：根據copula模型，抽樣每個馬爾科夫狀態空間的風速和風向累積概率分布值，使用逆抽樣法得到計及日特性的風速和風向樣本；

步驟3：計算最優季度系數，得到計及日特性和季特性的仿真樣本.

以上是一個季度的風矢量抽樣過程，改變年份和季度，可得到更多的修正風矢量仿真樣本.由于統計樣本有限，當仿真年限超過實測樣本長度時，將實測樣本重置為第一年，進行循環抽樣.

3 計及日特性和季特性含風能電力系統充裕度評估方法

3.1 發電機組兩狀態運行模型

為驗證計及風矢量日特性和季特性的風矢量仿真模型，擬將風電場并入IEEE-RTS79節點系統，使用時序蒙特卡羅法計算系統的充裕度指標.其中發電機停運模型為兩狀態馬爾科夫模型，通過“運行-停運-運行”的循環過程來模擬發電機的運行狀態[20].

圖6為發電機狀態轉移圖，其中λ和μ分別為機組的失效率和修復率，根據λ和μ可循環抽樣每臺機組的運行時間TF和修復時間TR，來模擬發電機的運行狀態.

圖6 發電機狀態轉移圖

(12)

(13)

式中：U和U′為(0，1)之間的均勻分布隨機數.

3.2 計及日特性和季特性的含風能發電系統充裕度評估流程

本文將風矢量仿真模型與時序蒙特卡羅模擬法結合，提出計及日特性和季特性含風能電力系統充裕度評估方法，使用缺電不足期望ELOE和缺電時間期望ELOL描述系統充裕度水平，詳細步驟如下.

步驟1：根據2.3節的方法得到計及日特性和季特性的修正風矢量仿真樣本；

步驟2：根據公式(12)和(13)抽樣得到每臺機組的運行狀態曲線.常規機組處于運行狀態時，其輸出功率為額定功率，風電機組的功率計算按照Jensen尾流模型[21]進行計算，當風機流入風速為vw時，其輸出功率為：

(14)

式中：vci、vr和vco分別為風機切入風速、額定風速和切出風速，查詢機組運行參數可知.a、b和c為恒定值，可由vci、vr和vco計算得到[21].

步驟3：將所有機組的輸出功率相加，得到含風電系統的發電功率曲線；

步驟4：逐小時比較發電功率曲線和負荷曲線，得到每小時缺電量ENS曲線和缺電持續時間TLD曲線；

步驟5：根據ENS和TLD計算系統充裕度指標ELOE和ELOL.

(15)

(16)

為確保仿真結果的準確性，以ELOE的方差系數作為收斂條件，當其方差系數小于0.05或仿真年限達到8 000年時結束仿真.

按照以上步驟可計算得到系統的充裕度指標，應用于具體算例中.

4 算例分析

4.1 計及日特性和季特性的風矢量仿真模型驗證

本文從NDAWN網站收集到美國North Dakota州Berthold觀測站2006年至2010年間的風矢量數據[16].應用本文收集到的Berthold觀測站的實測風矢量數據驗證本文模型，結合2.1節和2.2節的模型參數，采用最大似然估計法，分別計算Gaussian copula函數、Gumbel copula函數、Frank copula函數和Clayton copula函數的未知參數，抽樣得到計及日特性的模擬風矢量仿真樣本和計及日特性和季特性的修正風矢量仿真樣本.將仿真風矢量樣本與實測風矢量樣本的均值進行比較，見表1.

根據表1中的的對比結果，發現應用Gumbel copula函數后，可使得仿真風矢樣本的均值最接近實測風矢量樣本的均值，故選取Gumbel copula函數為風速和風向的聯合概率分布模型.

表1 不同copula函數類型的仿真風矢量樣本與實測風矢量樣本的均值對照表

由于數據分類較多，以時刻1：00的風矢量模型為例，該數據集的風矢量分為5個聚類空間，聚類空間之間的轉移概率見表2，每個聚類空間內風矢量的分布參數見表3.

表2 馬爾科夫轉移概率表

表3 模型參數表

本文以ARMA(自回歸滑動平均模型)[22]和一階馬爾科夫鏈模型[23]作為對照以驗證模型的日特性，由于數據分類較多，選取其中6個時刻的數據對比見表4和表5.

表4 不同風速時間序列的平均值比較(單位：m·s-1)

表5 不同風速時間序列的標準差比較(單位：m·s-1)

通過表4和表5可以看出，本文提出的計及日特性和季特性的風矢量修正模型在各時刻風速均值方面與歷史風速序列最為接近，在各時刻標準差方面較為接近，體現了風矢量的日特性.

以一年4個季度的模擬風矢量和仿真風矢量為樣本，前者僅計及風矢量的日特性，后者計及風矢量的日特性和季特性.將這兩組風矢量與實測數據進行對比，當公式(10)和公式(11)的權重系數均取值為0.5時，修正前后的風速數據對比如圖7所示.圖7中藍色數據為實測風矢量，棕色數據為只計及日特性的模擬風矢量，黃色數據為計及日特性和季特性的修正風矢量.通過表1可以看出，僅計及日特性的模擬風矢量數據的全年均值與歷史風矢量數據很接近.但觀察圖7可發現，模擬風矢量的每季度數據與實測數據存在一定差異，并不具備季節特征.當風速風向各自的均值和標準差權重系數均取0.5時，修正系數對每季度風速和風向均值的修正效果較為明顯，對風向標準差的修正效果較小.結果表明，最優季度修正系數的引入，使風矢量模型具備一定的季特性.

圖7 不同風矢量各季度對比圖

4.2 日特性和季特性對系統充裕度的影響

為評估季特性對系統充裕度的影響，本文假設North Dakota州Berthold觀測站的某座風電場由100臺相同型號風電機組組成，且該風電場接入IEEE-RTS79發電系統，該系統由32臺常規機組組成，總裝機容量3 405 MW，峰值負荷為2 850 MW，系統具體參數參見文獻[24].風電場內100臺風電機組按照10行10列的方式均勻布置，且相鄰風電機組間的間距取值12r，其中r為風電機組葉輪半徑40 m.風力發電機的額定功率為2 MW，其切入、額定和切出風速分別為3、12和25 m/s，風機的故障率和修復率分別為0.02和0.002 2.

風電場接入前，原始IEEE-RTS系統的ELOL和ELOE分別為9.37 h/a和1 197.4 MW·h/a.應用本文3.3節充裕度評估流程，評估以上含風能IEEE-RTS發電系統充裕度.為進行對比，本文也使用一階馬爾科夫鏈模型和ARMA模型作為風速仿真模型，評估該風能IEEE-RTS發電系統充裕度.三者對比結果見表6.

表6 不同風矢量模型對應的系統充裕度指標

可以看到：風電場接入后，系統的ELOL和ELOE指標均降低，表明風電接入對系統充裕度具有有益貢獻；如忽略風矢量季特性和尾流效應，將導致充裕度結果評估過于樂觀.在模型中加入風矢量的季特性，將對結果產生一定影響.

4.3 尾流效應對系統充裕度的影響

為評估尾流效應對系統充裕度的影響，本文選取4種風電場布局情形，每種情形下的風電場仍按照10行10列的方式均勻布置，但相鄰風電機組間距分別取值5r，10r，12r和15r.應用本文充裕度評估方法，評估不同布局情形下的系統充裕度結果見表7.

表7 不同風機間距下的系統充裕度指標

可以看到：風機間的尾流效應對系統充裕度有負面影響，且風機間距越小，尾流效應越強，對系統的充裕度影響越大.在不同風機間距條件下，應用修正風矢量進行系統充裕度評估得到的指標均低于模擬風矢量.結果表明在不同尾流條件下，忽略風矢量的季特性均可能導致過于樂觀的充裕度評估結果.

4.4 峰值負荷對系統充裕度的影響

以2 850 MW作為系統峰值負荷的基準情形，改變峰值負荷分別取2 750、2 950和3 050 MW，不同峰值負荷參數下的系統充裕度指標見表8.

表8 不同峰值負荷下的系統充裕度指標

從表8可以看出：風電場的接入顯著提高了電力系統的充裕度；系統峰值負荷是影響系統充裕度的重要因素，峰值負荷增大導致系統充裕度降低.在系統不同峰值負荷的條件下，兼顧日特性和季特性的修正風矢量模型與只計及日特性的模擬風矢量模型相比，其系統充裕度指標較低.結果表明在不同峰值負荷的條件下，忽略風矢量的季特性均可能使充裕度評估結果過于樂觀.

5 結 論

本文提出了一種計及日特性和季特性的風矢量仿真模型，該模型基于正交分解法對風矢量進行聚類分析以劃分馬爾科夫狀態空間，建立了計及日特性的風矢量馬爾科夫鏈模型，并提出兩組修正系數來描述風矢量的季節特征，對只計及日特性的仿真風矢量樣本加以修正.本文的風矢量模型在考慮風速風向相關性的基礎上，計及日周期特征與季節特征.

計及日特性的風矢量馬爾科夫鏈模型，在統計特征上較為接近歷史風矢量數據.提出一種同時修正仿真樣本均值和標準差的季節系數，使風矢量仿真樣本在保持統計特征的前提下呈現季節特征，從而提高了風矢量仿真模型的準確性.使用時序蒙特卡羅法對風矢量仿真模型加以應用，進一步證明了該模型的實用性和普適性.如果忽略風矢量的季節特征，使用未經修正的風矢量樣本應用于含風電系統的時序蒙特卡羅法，會使充裕度評估結果過于樂觀.本文驗證了季特性對系統充裕度評估的結果，但仍存在不足之處，地形條件和風機高度也可能對風速分布產生影響，需要進一步驗證以排除相關因素的干擾.