王進釗，嚴(yán)干貴，劉侃

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

0 引言

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

在一個由直驅(qū)風(fēng)機構(gòu)成的風(fēng)電場中，若對風(fēng)電機組構(gòu)建詳細(xì)的數(shù)學(xué)模型需要上千個微分方程，“自下而上”堆積木式的建模方法面臨“維數(shù)災(zāi)”問題，仿真計算代價大，并且海量仿真結(jié)果難以分析，尤其是機群在功率振蕩中的作用難以估量，使得搭建模型的有效性難以檢驗。因此，建立一個能夠捕獲高階系統(tǒng)重要動態(tài)的風(fēng)電場降階模型對系統(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)外開展了廣泛研究。在風(fēng)電場建模方面，主要有聚合建模和降階2 種。聚合建模又分為單機等值模型和多機等值模型，多用于穩(wěn)態(tài)潮流計算和暫態(tài)計算，較早一些用于振蕩分析的等值模型采用的是單機等值模型。由于次同步振蕩現(xiàn)象是“場-網(wǎng)”作用結(jié)果，采用單機等值模型后，對于源網(wǎng)耦合作用考慮不足，尤其當(dāng)電網(wǎng)強度較弱時，等值模型有效性值得商榷；另外，等值理論依據(jù)不足，如何證明風(fēng)電場等值前后主導(dǎo)振蕩特性一致也有待分析。在風(fēng)電場的降階模型中，最初采用基于聚合的方法，該方法可進一步分為單機聚合和多機聚合[13-15]。此外，文獻(xiàn)[16]應(yīng)用模態(tài)分析法，通過分析并保留特定需要的模式，模擬風(fēng)電機組輸出功率的降階。文獻(xiàn)[17]應(yīng)用奇異值攝動的方法對復(fù)雜電力系統(tǒng)進行模型降階研究。文獻(xiàn)[18]應(yīng)用Krylov 子空間類方法對系統(tǒng)進行模型降階研究。但是目前對于面向直驅(qū)風(fēng)電場次同步振蕩分析的模型降階研究鮮見報道，值得更加深入的研究。

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

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

1.1 直驅(qū)式永磁同步風(fēng)電機組結(jié)構(gòu)

直驅(qū)式永磁同步風(fēng)電機組結(jié)構(gòu)如圖1 所示，包括風(fēng)輪機、永磁同步發(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ū)永磁同步風(fēng)電系統(tǒng)示意圖Fig.1 Diagram of direct drive permanent magnet synchronous wind power system

1.2 風(fēng)力機的風(fēng)能捕獲模型

風(fēng)力機通過葉片捕獲風(fēng)能，葉片固定在風(fēng)力機輪轂處，假設(shè)輪轂處風(fēng)速為v0，則風(fēng)力機所捕獲的風(fēng)功率為

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

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

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

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

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

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

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

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

將旋轉(zhuǎn)坐標(biāo)系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ēng)電機組拓?fù)浣Y(jié)構(gòu)中，發(fā)電機側(cè)變流器控制結(jié)構(gòu)如圖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ù)學(xué)模型表示如下：

式中：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ū)式永磁同步風(fēng)電機組拓?fù)浣Y(jié)構(gòu)中，電網(wǎng)側(cè)變流器控制結(jié)構(gòu)如圖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ù)學(xué)模型表示如下：

式中：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ū)式永磁同步風(fēng)電機組拓?fù)浣Y(jié)構(gòu)中，直流側(cè)電容C的數(shù)學(xué)模型表示如下：

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

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

式中xi、ui、yi分別為第i臺直驅(qū)風(fēng)電機組的狀態(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。

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

式中：uw為所有直驅(qū)風(fēng)電機組的并網(wǎng)點電壓構(gòu)成的列向量；yw輸出電流構(gòu)成的列向量；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要遠(yuǎn)遠(yuǎn)小于原系統(tǒng)階數(shù)n，降階算法過程耗時較短，計算速度較快。

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

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

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

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

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

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

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

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

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

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

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

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

3 仿真算例分析

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

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

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

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

本文采用ADI 迭代的方法，以克服平衡截斷法求解李亞普諾夫方程過程緩慢的缺點，在進行降階之前，需要尋找恰當(dāng)?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臺風(fēng)機系統(tǒng)龐大，利用傳統(tǒng)方法求解李亞普諾夫方程一般需要37.83 s，本次仿真利用ADI 迭代方法求解李亞普諾夫方程則需要12.11 s，計算速度大大提高。

在Simulink 上搭建風(fēng)電場模型，仿真實現(xiàn)100 臺直驅(qū)風(fēng)電機組聯(lián)網(wǎng)仿真，利用Matlab 計算平臺自帶算法對風(fēng)電場系統(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)求解方法所得到的結(jié)果基本一致。同時，通過觀察圖4 可知，Hankel 奇異值在第3、6、11 階時有了顯著的下降，因此可將3、6、11階作為降階的目標(biāo)階數(shù)。

圖4 風(fēng)電場Hankel奇異值Fig.4 Hankel singular value of wind farm

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

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

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

圖5 風(fēng)電場仿真接線圖Fig.5 Wind farm simulation wiring diagram

3.2 仿真結(jié)果及分析

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

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

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

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

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

特征值分布中次同步振蕩模式對應(yīng)的幾何距離小于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ū)風(fēng)電機組處施加Δudcref=0.01 pu的階躍擾動，前后時域仿真波形基本吻合，雖然降階系統(tǒng)與全階系統(tǒng)存在一定誤差，但該誤差處于標(biāo)準(zhǔn)范圍內(nèi)，且降階系統(tǒng)與全階系統(tǒng)在次同步頻段內(nèi)具有相同的運行效果，能夠較好地等效全階系統(tǒng)在次同步頻段內(nèi)的響應(yīng)特性。

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

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

4 結(jié)論

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

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

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