基于交替方向隱式平衡截斷法的直驅(qū)風電場次同步振蕩分析的模型降階研究

王進釗，嚴干貴，劉侃

（現(xiàn)代電力系統(tǒng)仿真控制與綠色電能新技術教育部重點實驗室(東北電力大學），吉林省 吉林市 132000）

0 引言

近年來，構建以新能源為主體的新型電力系統(tǒng)被作為能源電力行業(yè)實現(xiàn)碳達峰、碳中和的主要措施[1-4]。我國風電發(fā)展迅猛，截至2021年11月14日，我國風電并網(wǎng)裝機容量達到30 015萬kW，突破3 億kW 大關，15 年間裝機容量增加了近百倍，年均增速30%左右。據(jù)行業(yè)權威機構預測：風電并網(wǎng)裝機容量到2030 年將達8 億kW，到2060 年將達30 億kW。2013 年以來，我國河北沽源風電系統(tǒng)屢次發(fā)生3～10 Hz 的次同步振蕩(subsynchronous oscillation，SSO)，曾造成上千臺風電機組異常脫網(wǎng)。2015 年7 月1 日，我國新疆哈密地區(qū)風電集群發(fā)生跨越5 個電壓等級的次/超同步振蕩事故，造成花園電廠3臺660 MW機組全切的嚴重后果。近年來，風力發(fā)電的滲透率占比持續(xù)攀升，使得風電并網(wǎng)引發(fā)的新型次同步現(xiàn)象不斷危害電網(wǎng)的安全穩(wěn)定運行[5-10]。

在一個由直驅(qū)風機構成的風電場中，若對風電機組構建詳細的數(shù)學模型需要上千個微分方程，“自下而上”堆積木式的建模方法面臨“維數(shù)災”問題，仿真計算代價大，并且海量仿真結果難以分析，尤其是機群在功率振蕩中的作用難以估量，使得搭建模型的有效性難以檢驗。因此，建立一個能夠捕獲高階系統(tǒng)重要動態(tài)的風電場降階模型對系統(tǒng)研究具有重要意義。模型降階方法在滿足問題分析需要的同時，將一個較大系統(tǒng)轉(zhuǎn)化為一個近似的較小系統(tǒng)，使得降階后的系統(tǒng)仿真時長和方程維數(shù)大大降低[11]。常見的模型降階方法有如下3種[12]：

1）基于奇異值分解(singular value decomposition，SVD)的方法。該方法主要包括平衡截斷(balanced truncation，BT)法和基于交替方向隱式(alternating direction implicit，ADI)的平衡截斷法。

2）Krylov子空間類方法。該方法主要包括有理克雷洛夫(rational Krylov，RK)算法和迭代有理克雷洛夫算法(iterative RK algorithm，IRKA)。

3）模態(tài)截斷方法。該方法識別系統(tǒng)特定的模式，并且將模式保留下來。

針對實際出現(xiàn)的振蕩問題，國內(nèi)外開展了廣泛研究。在風電場建模方面，主要有聚合建模和降階2 種。聚合建模又分為單機等值模型和多機等值模型，多用于穩(wěn)態(tài)潮流計算和暫態(tài)計算，較早一些用于振蕩分析的等值模型采用的是單機等值模型。由于次同步振蕩現(xiàn)象是“場-網(wǎng)”作用結果，采用單機等值模型后，對于源網(wǎng)耦合作用考慮不足，尤其當電網(wǎng)強度較弱時，等值模型有效性值得商榷；另外，等值理論依據(jù)不足，如何證明風電場等值前后主導振蕩特性一致也有待分析。在風電場的降階模型中，最初采用基于聚合的方法，該方法可進一步分為單機聚合和多機聚合[13-15]。此外，文獻[16]應用模態(tài)分析法，通過分析并保留特定需要的模式，模擬風電機組輸出功率的降階。文獻[17]應用奇異值攝動的方法對復雜電力系統(tǒng)進行模型降階研究。文獻[18]應用Krylov 子空間類方法對系統(tǒng)進行模型降階研究。但是目前對于面向直驅(qū)風電場次同步振蕩分析的模型降階研究鮮見報道，值得更加深入的研究。

本文建立一種面向直驅(qū)風電場次同步振蕩分析的降階等值模型，采用基于ADI 的平衡截斷方法對風電場模型進行降階分析，在直驅(qū)風電場次同步振蕩穩(wěn)定場景下驗證該降階方法的有效性，從降階模型階數(shù)、計算耗時等方面對比該降階方法的應用效果。

1 直驅(qū)永磁同步風電場模型的建立

1.1 直驅(qū)式永磁同步風電機組結構

直驅(qū)式永磁同步風電機組結構如圖1 所示，包括風輪機、永磁同步發(fā)電機(permanent magnet synchronous generator，PMSG)、四象限型變流器等[14]。其中：AC 表示交流電(alternating current)；DC 表示直流電(direct current)；PWM 表示脈沖寬度調(diào)制(pulse width modulation)；u1和i1分別為發(fā)電機側(cè)的電壓和電流；u2和i2分別為電網(wǎng)側(cè)的電壓和電流；Udc和Udcref分別為直流側(cè)兩端的電壓及其參考值；P1，ref和Q1，ref分別為發(fā)電機側(cè)輸出的有功功率、無功功率參考值；Q2，ref為電網(wǎng)側(cè)輸出的無功功率參考值；ωr為發(fā)電機轉(zhuǎn)子轉(zhuǎn)速。

圖1 直驅(qū)永磁同步風電系統(tǒng)示意圖Fig.1 Diagram of direct drive permanent magnet synchronous wind power system

1.2 風力機的風能捕獲模型

風力機通過葉片捕獲風能，葉片固定在風力機輪轂處，假設輪轂處風速為v0，則風力機所捕獲的風功率為

式中：Cp為風能利用系數(shù)，是風力機葉片從自然風能中所吸收能量與葉片掃過面積內(nèi)未擾動氣流所具有動能之比，是表征風力機捕獲風能大小的一個量；Pw=，為風功率，其中m為風質(zhì)量；ρ為空氣密度；A1=πR2，為風力機葉片掃掠面積，其中R為風力機葉片半徑。

在風速恒定的情況下，風機捕獲的機械功率取決于Cp的大小。對于變速風力發(fā)電機組，Cp直接與葉尖速比λ、槳距角β有關，表達式如下：

槳距角控制系統(tǒng)數(shù)學模型可用以下方程表示：

式中：β0為初始槳距角；Tsq為慣性時間常數(shù)。

1.3 直驅(qū)式永磁同步風電機組模型

為建立永磁同步發(fā)電機數(shù)學模型，對定子繞組回路電壓電流、磁鏈正方向進行以下規(guī)定：定子繞組電壓、電流正方向為非關聯(lián)參考方向；定子負值電流產(chǎn)生正值磁鏈。

兩相旋轉(zhuǎn)坐標系下電壓方程為：

式中：U1d、U1q分別為d、q軸定子繞組電壓；Rs為定子繞組電阻；i1d、i1q分別為d、q軸定子繞組電流；ψsd、ψsq分別為d、q軸定子繞組磁鏈。

將旋轉(zhuǎn)坐標系d軸定在轉(zhuǎn)子磁鏈方向上，定子磁鏈可表示為：

式中：Lsd、Lsq分別為d、q軸定子繞組電感；φ為轉(zhuǎn)子磁鏈。

將式(6)代入式(5)可得：

根據(jù)牛頓運動定律，轉(zhuǎn)子運動方程為

式中：Tm為原動機加于電機軸的機械力矩；Te為發(fā)電機電磁力矩；Ωm為轉(zhuǎn)子機械角速度；J為轉(zhuǎn)子的轉(zhuǎn)動慣量。

在直驅(qū)式永磁同步風電機組拓撲結構中，發(fā)電機側(cè)變流器控制結構如圖2 所示，一般采用基于轉(zhuǎn)子磁鏈定控制策略實現(xiàn)d、q軸解耦控制，由發(fā)電機側(cè)輸出的有功功率、無功功率參考值經(jīng)過PI控制器之后，形成d、q軸電流分量參考值i1dref、i1qref，d、q軸電流分量及其參考值經(jīng)過PI 控制器之后，得到d、q軸機端電壓分量。

圖2 發(fā)電機側(cè)變流器控制框圖Fig.2 Control block diagram of generator side converter

發(fā)電機側(cè)變流器的數(shù)學模型表示如下：

式中：x1、x2、x3和x4為發(fā)電機側(cè)變流器的控制參數(shù)；P1和Q1分別為發(fā)電機側(cè)變流器輸出的有功功率和無功功率；ΔPe和ΔQe分別為發(fā)電機側(cè)變流器輸出的有功功率增量和無功功率增量；Kp1,Ki1,Kp2,Ki2,Kp3,Ki3,Kp4,Ki4為PI控制器的控制參數(shù)；Ls為發(fā)電機定子電感。

在直驅(qū)式永磁同步風電機組拓撲結構中，電網(wǎng)側(cè)變流器控制結構如圖3 所示，采用并網(wǎng)點電壓矢量控制，直流電壓及其參考值經(jīng)過PI控制器后，形成d軸電流分量參考值i2dref，d軸電流分量及其參考值經(jīng)過PI 控制器后，得到d軸網(wǎng)側(cè)電壓分量；網(wǎng)側(cè)無功功率及其參考值經(jīng)過PI控制器后，形成q軸電流分量參考值i2qref，q軸電流分量及其參考值經(jīng)過PI控制器后，得到q軸網(wǎng)側(cè)電壓分量。

圖3 電網(wǎng)側(cè)變流器控制框圖Fig.3 Control block diagram of grid side converter

電網(wǎng)側(cè)變流器的數(shù)學模型表示如下：

式中：U2d和U2q分別為電網(wǎng)側(cè)變流器的d、q軸電壓；i2d和i2q分別為電網(wǎng)側(cè)變流器的d、q軸電流；x5、x6、x7和x8為電網(wǎng)側(cè)變流器的控制參數(shù)；P2和Q2分別為電網(wǎng)側(cè)變流器輸出的有功功率和無功功率；ΔQg為電網(wǎng)側(cè)變流器輸出的無功功率增量；ΔUdc為直流側(cè)兩端的電壓增量；Ug為電網(wǎng)電壓；Kp5,Ki5,Kp6,Ki6,Kp7,Ki7,Kp8,Ki8為PI控制器的控制參數(shù)；Lg為電網(wǎng)側(cè)變流器進線電抗器的電感；ωg為電網(wǎng)同步角轉(zhuǎn)速。

在直驅(qū)式永磁同步風電機組拓撲結構中，直流側(cè)電容C的數(shù)學模型表示如下：

1.4 并網(wǎng)風電場線性化狀態(tài)空間模型

i臺直驅(qū)風電機組的線性化狀態(tài)空間可表示為：

式中xi、ui、yi分別為第i臺直驅(qū)風電機組的狀態(tài)變量、并網(wǎng)點電壓和輸出電流。狀態(tài)向量x=(β，ωr，x1，x2，x3，x4，Udc，x5，x6，x7，x8)T；輸入向量u=(U1d，U1q，U2d，U2q)T；輸出向量y=(i1d，i1q，i2d，i2q)T。

忽略集電線路的電阻，風電場網(wǎng)絡方程可表示為

式中：uw為所有直驅(qū)風電機組的并網(wǎng)點電壓構成的列向量；yw輸出電流構成的列向量；Zk為集電線路的阻抗。

2 基于ADI的平衡截斷方法

線性化后的狀態(tài)空間表示如下：

式中：A∈Rn×n為狀態(tài)矩陣；b∈Rn×t為輸入矩陣；cT∈Rm×n為輸出矩陣；x(t)∈Rn為系統(tǒng)的狀態(tài)向量；u(t)∈Rt為系統(tǒng)的輸入向量；y(t)∈Rm為系統(tǒng)的輸出向量。

原始系統(tǒng)(22)根據(jù)變換矩陣V和W可以得到其降階系統(tǒng)，一般表示為：

式(24)的傳遞函數(shù)可表示為

得到其降階系統(tǒng)后，式(25)的傳遞函數(shù)可表示為

一般情況下，傳遞函數(shù)Hk(s)是傳遞函數(shù)H(s)的近似，降階系統(tǒng)與原系統(tǒng)在次同步頻段內(nèi)保持相同的性質(zhì)，并且降階系統(tǒng)的階數(shù)k要遠遠小于原系統(tǒng)階數(shù)n，降階算法過程耗時較短，計算速度較快。

平衡截斷方法[19-22]最早由Moore提出，該方法首先通過計算系統(tǒng)的可控性和可觀性Gramian 矩陣，得到系統(tǒng)的Hankel 奇異值(Hankel single values，HSV)，并以此為標準找到對系統(tǒng)輸入、輸出影響較小的狀態(tài)；然后構造平衡變換矩陣，對系統(tǒng)進行平衡變換；最后將對系統(tǒng)輸入、輸出影響較小的狀態(tài)從平衡變換后的系統(tǒng)中截斷，從而達到模型降階的效果。平衡截斷方法實現(xiàn)降階過程的步驟如下：

1）通過解李亞普諾夫方程(28)得到可控的格萊姆矩陣P與可觀的格萊姆矩陣Q。

2）根據(jù)P、Q以及式(29)，計算全部的奇異值σi(i=1，2，…，n)。

當存在一個正整數(shù)q，使得σq遠大于σq+1時，q就是理想降階階數(shù)。

3）應用式(30)進行Cholesky 分解，根據(jù)得到的Lc、Lo以及式(31)，進行SVD分解。

4）計算平衡矩陣T，并對狀態(tài)矩陣A、輸入矩陣b和輸出矩陣cT進行平衡變換，得到。對平衡變換后的系統(tǒng)狀態(tài)矩陣進行截斷，得到降階系統(tǒng)狀態(tài)空間。

平衡截斷方法的優(yōu)點在于有明確的指標確定階數(shù)，同時降階系統(tǒng)可以保持原系統(tǒng)的可控性、可觀性以及穩(wěn)定性。但是其缺點是求解李亞普諾夫方程過程非常緩慢，這在處理大型風電場系統(tǒng)時不可行，利用低階矩陣近似求解P與Q，可以解決這個問題，這種方法稱為ADI。

ADI 迭代最初用于求解一個線性方程組，為了符合ADI迭代規(guī)律，這里將式(28)改為

因此，這個迭代方程表示如下：

式中αi為移位參數(shù)，將方程(33)合并成一個方程：

此外，在平衡截斷方法中，目標是計算P的低秩因子分解，假定Lc，0=0且Lc，iLc，iT=Pi，式(34)可以直接寫成：

同理，按照以上步驟也可計算Q的低秩因子分解。

3 仿真算例分析

為了驗證基于ADI的平衡截斷方法的實用性，本次仿真選取并聯(lián)式集電線路的100 臺金風1.5 MW 直驅(qū)風電機組來組成算例風電場，其中100 臺直驅(qū)風電機組參數(shù)及控制策略完全相同。風速在4～8 m/s，以0.01 m/s 為步長隨機分布；集電線路長度在400～800 m，以1 m 為步長隨機分布。在式(22)的基礎上，以場內(nèi)第1臺直驅(qū)風機的直流電壓擾動Δudcref1作為輸入，將風電場總輸出功率增量ΔPΣ作為輸出，可得如下狀態(tài)空間：

式中：ΔX為全系統(tǒng)狀態(tài)變量增量；Ar為狀態(tài)矩陣；Br為輸入矩陣；Cr為輸出矩陣。

構建面向直驅(qū)風電場次同步振蕩分析的降階模型并應用在仿真分析中，通過對比系統(tǒng)降階前后次同步振蕩模式、Bode圖及時域仿真波形的吻合情況，驗證該降階方法的有效性，從降階模型階數(shù)、計算耗時等方面分析該降階方法的實用性。

3.1 直驅(qū)風電場次同步振蕩分析的模型降階

本文采用ADI 迭代的方法，以克服平衡截斷法求解李亞普諾夫方程過程緩慢的缺點，在進行降階之前，需要尋找恰當?shù)囊莆粎?shù)αi以保證在ADI 迭代中精度最高且速度最小。計算環(huán)境采用DESKTOP-232GOKN，Windows 10 系統(tǒng)，4 GB內(nèi)存，Intel(R) Core(TM) i5-6300HQ的CPU。通過大量算法演示，最終確定移位參數(shù)αi為39 時，ADI迭代效果最佳。100臺風機系統(tǒng)龐大，利用傳統(tǒng)方法求解李亞普諾夫方程一般需要37.83 s，本次仿真利用ADI 迭代方法求解李亞普諾夫方程則需要12.11 s，計算速度大大提高。

在Simulink 上搭建風電場模型，仿真實現(xiàn)100 臺直驅(qū)風電機組聯(lián)網(wǎng)仿真，利用Matlab 計算平臺自帶算法對風電場系統(tǒng)系數(shù)矩陣AW的Hankel奇異值進行求解，如圖4(a)所示。同時，利用ADI迭代方法對AW的Hankel奇異值進行求解，如圖4(b)所示。圖4(c)為圖4(b)的局部放大圖。由圖4(a)、(b)可見，2 幅圖的數(shù)據(jù)分布情況大致相同，這說明用ADI 迭代方法求解李亞普諾夫方程與傳統(tǒng)求解方法所得到的結果基本一致。同時，通過觀察圖4 可知，Hankel 奇異值在第3、6、11 階時有了顯著的下降，因此可將3、6、11階作為降階的目標階數(shù)。

圖4 風電場Hankel奇異值Fig.4 Hankel singular value of wind farm

通過Hankel 奇異值得到目標階數(shù)之后，分別將降階得到的3、6、11階風電場等值降階系統(tǒng)與升壓變壓器相連并接入電網(wǎng)中，圖5 為風電場仿真接線圖，表1為直驅(qū)風電機組參數(shù)，表2為風電場模型參數(shù)。

表1 直驅(qū)風電機組參數(shù)Tab.1 Direct drive wind turbine parameters

表2 風電場模型參數(shù)Tab.2 Wind farm model parameters

圖5 風電場仿真接線圖Fig.5 Wind farm simulation wiring diagram

3.2 仿真結果及分析

6階降階系統(tǒng)的狀態(tài)空間表示如下：

圖6 為全階系統(tǒng)與風電場降階系統(tǒng)Bode 圖對比曲線。其中，6 階與11 階風電場降階系統(tǒng)在次同步頻段內(nèi)與全階系統(tǒng)基本吻合，而3 階降階系統(tǒng)與原系統(tǒng)相差較多，3、6、11 階風電場降階系統(tǒng)都能夠捕獲頻率的峰值40.445 1 Hz，因此降階目標階數(shù)選擇6階。

圖6 全階系統(tǒng)與降階系統(tǒng)Bode圖Fig.6 Bode diagram of full order system and reducedorder system

為了更好地驗證降階模型保持系統(tǒng)的主導振蕩特性，本文以全階模型和降階模型對應特征根的平均歐氏距離作為衡量降階模型準確度的指標。降階模型的主導特征根應當與全階模型的主導特征根盡可能地接近，即平均歐氏距離越小，降階模型的降階效果越好，降階的準確度越高。平均歐氏距離的表達式為

式中：λm，e和λm，i分別為降階模型和全階模型的主導特征根；n為振蕩特征根的數(shù)量。

特征值分布中次同步振蕩模式對應的幾何距離小于10-7，次同步振蕩模式得到了保持。圖7為6階與全階系統(tǒng)的特征值對比。圖8為6階與全階系統(tǒng)的時域波形對比。

圖7 6階與全階系統(tǒng)的特征值對比Fig.7 Comparison of eigenvalues of 6th order and full order systems

圖8 6階與全階系統(tǒng)時域波形對比Fig.8 Comparison of time-domain waveforms of 6th order and full order systems

從圖7可以看出，在特征值分布中，全階模型的次同步振蕩模式為-1.932 4+239.287j，降階模型的次同步振蕩模式為-1.935 6+239.287j，說明次同步頻段內(nèi)的特征根在降階系統(tǒng)中得到了保留，篩選出了對系統(tǒng)次同步特性影響較大的特征根。從圖8可以看出，t=0.5 s時在第一臺直驅(qū)風電機組處施加Δudcref=0.01 pu的階躍擾動，前后時域仿真波形基本吻合，雖然降階系統(tǒng)與全階系統(tǒng)存在一定誤差，但該誤差處于標準范圍內(nèi)，且降階系統(tǒng)與全階系統(tǒng)在次同步頻段內(nèi)具有相同的運行效果，能夠較好地等效全階系統(tǒng)在次同步頻段內(nèi)的響應特性。

表3 為利用基于ADI 的平衡截斷方法得到的降階模型與全階模型的對比。可以看出，利用基于ADI 的平衡截斷方法得到的降階模型，其計算次同步振蕩模式和時域仿真的時間均遠小于全階模型，并且將模型階數(shù)從1 200 階降到6 階；同時，利用ADI 迭代方法求解李亞普諾夫方程，有效縮短了仿真計算時間，提高了計算效率。

表3 降階模型與全階模型的比較Tab.3 Comparison between reduced order model and full order model

4 結論

利用基于ADI 的平衡截斷方法，對直驅(qū)式永磁同步風電機組的風電場進行模型降階研究，提出了面向直驅(qū)風電場次同步振蕩分析的降階模型，得到以下結論：

1）根據(jù)Hankel 奇異值以及全階系統(tǒng)與風電場降階系統(tǒng)Bode 圖的吻合情況，與全階模型相比，降階模型將階數(shù)從1 200階降為6階，同時在次同步頻段內(nèi)其響應與全階模型基本相同。

2）選擇合適的移位參數(shù)，利用ADI迭代的方法求解李亞普諾夫方程比傳統(tǒng)方法的求解時間縮短了近1/3，提高了降階速度。