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

王進釗，嚴干貴，劉侃

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

0 引言

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

在一個由直驅風機構成的風電場中，若對風電機組構建詳細的數學模型需要上千個微分方程，“自下而上”堆積木式的建模方法面臨“維數災”問題，仿真計算代價大，并且海量仿真結果難以分析，尤其是機群在功率振蕩中的作用難以估量，使得搭建模型的有效性難以檢驗。因此，建立一個能夠捕獲高階系統重要動態的風電場降階模型對系統研究具有重要意義。模型降階方法在滿足問題分析需要的同時，將一個較大系統轉化為一個近似的較小系統，使得降階后的系統仿真時長和方程維數大大降低[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）模態截斷方法。該方法識別系統特定的模式，并且將模式保留下來。

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

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

1 直驅永磁同步風電場模型的建立

1.1 直驅式永磁同步風電機組結構

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

圖1 直驅永磁同步風電系統示意圖Fig.1 Diagram of direct drive permanent magnet synchronous wind power system

1.2 風力機的風能捕獲模型

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

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

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

槳距角控制系統數學模型可用以下方程表示：

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

1.3 直驅式永磁同步風電機組模型

為建立永磁同步發電機數學模型，對定子繞組回路電壓電流、磁鏈正方向進行以下規定：定子繞組電壓、電流正方向為非關聯參考方向；定子負值電流產生正值磁鏈。

兩相旋轉坐標系下電壓方程為：

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

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

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

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

根據牛頓運動定律，轉子運動方程為

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

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

圖2 發電機側變流器控制框圖Fig.2 Control block diagram of generator side converter

發電機側變流器的數學模型表示如下：

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

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

圖3 電網側變流器控制框圖Fig.3 Control block diagram of grid side converter

電網側變流器的數學模型表示如下：

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

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

1.4 并網風電場線性化狀態空間模型

i臺直驅風電機組的線性化狀態空間可表示為：

式中xi、ui、yi分別為第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。

忽略集電線路的電阻，風電場網絡方程可表示為

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

2 基于ADI的平衡截斷方法

線性化后的狀態空間表示如下：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 仿真算例分析

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

式中：ΔX為全系統狀態變量增量；Ar為狀態矩陣；Br為輸入矩陣；Cr為輸出矩陣。

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

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

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

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

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

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

表1 直驅風電機組參數Tab.1 Direct drive wind turbine parameters

表2 風電場模型參數Tab.2 Wind farm model parameters

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

3.2 仿真結果及分析

6階降階系統的狀態空間表示如下：

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

圖6 全階系統與降階系統Bode圖Fig.6 Bode diagram of full order system and reducedorder system

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

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

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

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

圖8 6階與全階系統時域波形對比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，說明次同步頻段內的特征根在降階系統中得到了保留，篩選出了對系統次同步特性影響較大的特征根。從圖8可以看出，t=0.5 s時在第一臺直驅風電機組處施加Δudcref=0.01 pu的階躍擾動，前后時域仿真波形基本吻合，雖然降階系統與全階系統存在一定誤差，但該誤差處于標準范圍內，且降階系統與全階系統在次同步頻段內具有相同的運行效果，能夠較好地等效全階系統在次同步頻段內的響應特性。

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

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

4 結論

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

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

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