基于耦合振子模型的含風電電力系統聚類同步

李靜宇 周明 李庚銀

關鍵詞： 圖論; Kuramoto模型; Kron簡化; 聚類同步; 耦合振子; 風電

中圖分類號： TN99?34 ? ? ? ? ? ? ? ? ? ? ? ? ?文獻標識碼： A ? ? ? ? ? ? ? ? ? ? ? ? ? 文章編號： 1004?373X（2019）05?0135?06

Clustering synchronization of wind power system based on coupled oscillator model

LI Jingyu， ZHOU Ming， LI Gengyin

（State Key Laboratory of New Energy Power System， North China Electric Power University， Beijing 102206， China）

Abstract： A network coupled clustering synchronization method based on Kuramoto model is proposed to solve the serious system synchronization problem caused by the disturbance of the wind power electric system. The graph theory is used to establish the power network contraction model. On this basis， the electric system equivalent model based on Kuramoto coupled oscillator model is established， and a power network clustering algorithm based on dynamic coupling strength is proposed. This method can quickly and accurately cluster the synchronous state of oscillator in the system after system disturbance. The 39?node New England system with wind farms is verified. The simulation results show that the cluster analysis can understand the synchronization state of the power grid clearly， and avoid the further collapse of the system.

Keywords： graph theory; Kuramoto model; Kron reduction; clustering synchronization; coupled oscillator; wind power

0 ?引 ?言

隨著電力系統的不斷互聯，電網的規模變得越來越復雜，而新能源發電技術的發展，也使得越來越多的間歇性電源接入到了電力系統中。由于以風電為代表的間歇性電源具有很強的隨機性和波動性，且通過電力電子設備與電網相連，對電網缺少必要的慣性支撐。因此，對含風電的電力系統遭受嚴重暫態擾動后系統保持同步的能力提出了更高的要求[1]。目前常用的電力系統同步能力評價的方法大致可以分為：采用數值積分手段直接求解系統的狀態方程，觀察系統受擾后狀態變量的變化軌跡的時域仿真法[2]，根據動力學系統的穩定性理論來構造系統的暫態能量函數，并由暫態能量函數的性質判斷受擾系統穩定性的直接法[3]。傳統電力系統同步能力評價方法需要先確定系統受到擾動后的臨界能量或不平衡點，當對大系統求解時方程維度高，求解難度大，多機系統適應性差。通過對系統進行等值雖然可以減小計算復雜度，但系統很多重要的節點信息也隨著系統等值被忽略。近年來隨著復雜網絡同步理論以及圖論理論的發展，在最大限度保留網絡中節點和支路信息的情況下，計及系統中各節點之間的耦合關系對系統受到擾動后聚類同步過程的影響逐漸成為研究熱點[4]。

1 ?基于圖論的電力系統表征

在圖論中通過定義一個無向圖[G]，包含一個[n]維節點集[V]和一個[e]維弧集[E]來表征網絡所包含的信息。如果節點[vi]和節點[vj]由一條弧進行連接，即[e=（vi，vj）]。那么就稱作節點[vi]是節點[vj]的一個鄰居節點。如果圖[G]中任意兩對節點可以通過一條弧進行連接，那么該圖形就是強連接的。因此，從網絡拓撲的角度可以建立電力系統的電力網絡保持模型或電力網絡收縮模型。而對于任意一個電力網絡都可用其鄰接矩陣表示，鄰接矩陣反映節點與弧之間的連接關系。在進行潮流計算時可用式（1）表示：

[Ax=y] （1）

式（1）與無向圖存在如下映射關系。[V]中每一個節點[vi]與[x，y]中每一個元素[xi，yi]形成雙映射。同時無向圖[G=][（V，E）]存在映射[φ]：[E→R]，使得[ωij=φ（eij）]，則稱[ωij]為弧[eij=（Vi，Vj）∈E]的權重，并用[G=（V，E，W）]表示加權圖[5]。與圖形相關的[N]維無向圖[G]的鄰接矩陣[A=aij∈Rn×n]表示為：

[aij=1， ? ? （vi，vj）∈E0， ? ?其他] （2）

圖的度矩陣是一個[n×n]的對角線矩陣[Δ=diag（degout（Vi））]，對角線元素表示各節點的出度，而其他元素均為零。圖[G]的定向可以為每條邊賦予一個方向，所以邊[（i，j）]是一條從節點[i]到節點[j]的弧線。定向為[σ]的圖[G]可以表示為[Gσ]。定向圖[Gσ]的關聯矩陣[B（Gσ）]是一個[n×e]矩陣，該矩陣的行和列可以分別通過圖[G]的節點及與之相連的弧進行表示。如果節點[j]通過弧進入節點[i]，那么[B（Gσ）]的[（i，j）]就等于1;如果弧是從節點[i]進入到節點[j]，則等于-1;其他情況則為0。

由于電網拓撲結構中的節點和邊是不同質的元件，采用有權無向網絡可以更好地表達電網中不同元件的獨特信息即[aij=aji=ωij]，分析過程中通常采用導納作為邊權。因此，通過將網絡的中節點分為邊界節點和內部節點，采用Kron簡化對內部節點進行消除，簡化后的網絡不僅保留了原來電力網中同步發電機終端電壓和電流之間的關系，同時節點數量得到了大幅簡化。網絡中的電流平衡方程可以寫為：

[I1I2=Q11Q12Q21Q22U1U2] （3）

式中：[I1]為外部節點的電流向量;[I2]為內部節點的電流向量;[U1]為外部節點的電壓向量;[U2]為內部節點的電壓向量;[Q11]為外部節點導納矩陣;[Q22]為內部節點導納矩陣;[Q12]，[Q21]為內外節點互導納矩陣。通過消去內部電壓節點[U2]可以得出簡化后的電流平衡方程為：

[I1-Q12Q-122I2=（Q11-Q12Q-122Q21）U1] （4）

令矩陣[Q11-Q12Q-122Q21=Qred]為[Q]中內部節點的爾補碼，因此[Qred]為相對于電力網絡圖[G]外部節點[U1]的Kron簡化矩陣。利用圖形理論基本原則，可以將電力系統建模成如圖1所示的電力系統收縮網絡。

因此將含風電的電力系統同步問題作為一個特殊的廣義同步問題，通過Kuramoto模型進行分析，就可以通過Kron簡化實現所有發電機之間的完全互連。

2 ?Kuramoto模型

2.1 ?耦合振子模型

一個包含[N]個耦合振子的網絡，其動力學方程可以通過Kuramoto模型描述為：

[θi=ωi-j=1nKijsin（θi-θj）] （5）

式中：[θi]表示振子相位;[K]表示耦合強度;[ωi]表示振子固有頻率。耦合強度是表征系統同步的重要參數，如果系統耦合強度弱，則系統中振子就會以與自身固有頻率不同步的方式運行;反之，如果耦合強度超過相位差保持穩定的臨界值，則會導致振子之間耦合，進而實現同步。對于包含[N]個振子的系統而言，可以通過定義全局序參量[r]，對系統的同步進行解釋。假定[N]個振子是一個圓上的各個點，且這些點會隨不同角速度進行移動，那么在完全同步的情況下，所有點都會同步地圍繞著這個圓進行移動，如圖2所示。作為同步自然測度的各點中心的有序參數可以被定義為：

[r（t）e（i?（t））=1Nj=1ne（iθj（t））] （6）

式中：[r（t）∈]（0，1），用于衡量這個網絡中振子的相位同步程度;[?（t）]為平均相位。[r]值反映了系統內處于同步簇中的節點數量占整個系統中節點數量的比例，[r]越接近1，表明網絡中的振子相位同步程度越高，[r=0]表示網絡中的振子相位不同步[6]。將式（6）左右兩邊同時乘以[e-iθi]后代入到式（5）中，得到通過全局序參數表示的Kuramoto模型：

[θi=ωi-Krsin（?i-θi）] （7）

由式（7）可知：隨著耦合強度[K]趨于零，各振子會根據其自由頻率進行震蕩，進而導致[r]趨于0;反之，當耦合強度[K]趨近于無窮大時，振子會逐漸與它們的平均相位同步，即[r]逐漸趨近于1。

2.2 ?基于Kuramoto模型電力系統等值模型

一個多機電力系統的詳細動態模型包括有轉子擺動、轉子磁鏈、動態控制等耦合非線性模型。因此，通過建立完整上述模型來表示整個電力系統是非常困難的。對于具備非轉移電導的多機電力系統而建立的Kuramoto模型，各發電機可以通過潮流進行耦合，將同步問題視為暫態穩定問題的一個特例。因此，可以從發電機動力學方程入手，建立基于Kuramoto模型的電力系統等值模型[7]。

一個由[n]臺發電機組構成的電力系統中，第[i]臺同步發電機的轉子運動方程可表示為：

[Hiπfδi+KDiπfδi=PTi-PEi] （8）

式中：[Hi，δi，KDi，PTi，PEi]表示系統中第[i]臺同步發電機的慣性時間常數、轉子角度、阻尼系數、輸入的機械功率以及輸出的電磁功率的標幺值; [f]表示同步頻率。

令[Mi=2Hiω0，Di=KDiω0]，則式（8）變為：

[Miδi=PTi-PEi-Diδi] （9）

在電力網絡收縮模型中，可以將系統中的同步發電機視為電壓源，在發生短路故障后通過暫態電抗[x′di]與電網中相應的同步發電機節點相連。設[E′i]為受到擾動后的電勢，同步發電機發生短路故障后的電勢可寫成如下形式：

[E′i=Vi+jx′diIi=Vi+jx′di（PGi-jQGiVi）] ?（10）

同時，在電力系統網絡收縮模型中，為了分析發電機節點的動態特性，通常將負荷視為恒定阻抗。對于由[n]個發電機節點和[m]個負荷節點組成的系統，通過在每個發電機節點處增加一個發電機內節點，以此來形成增廣網絡，此時網路中節點總數[N=2n+m]。考慮到在增廣網絡中負荷不僅存在于負荷節點，還有可能在所有發電機的機端節點處[8]。因此，電力系統網絡收縮模型中的負荷導納為：

[yLi=（PLi-jQLi）Vi2] （11）

式中：[Pi，Qi]分別為節點[i]處的有功、無功負荷;[Vi]為穩態時節點[i]電壓的模。

令[Y′di=diag（yd）]代表發電機暫態導納矩陣;[Y′GL=diag（yGL）]代表發電機機端節點所帶負荷等值之后的導納矩陣;[Y′LL=diag（yLL）]代表負荷節點所形成的導納矩陣。因此增廣網絡導納矩陣可表示為：

[Y=Y′d-Y′d0-Y′dYGG+Y′d+Y′GLYGL0YLGYLL+Y′LL] （12）

令[Y11=Y′d]，[Y12=-Y′d0T]，[Y21=-Y′d0T]，[Y22=YGG+Y′d+Y′GLYGLYLGYLL+Y′LL]。

因此得到修改后網絡的節點電壓方程如下：

[In0=Y11（n×n）Y21（m×n） ? ?Y12（n×m）Y22（m×m）VnVm] （13）

使用Kron簡化的方法消去輸出電流為0的節點來對導納矩陣進行化簡，得到僅由[n]個外部發電機內節點構成網絡，其內部則是全聯通的。簡化后的導納陣為：

[Yred=Y1-Y12Y-122Y21] （14）

用[Yij，red]表示簡化后導納陣[Yred]中第[i]行、第[j]列的元素，則有：

[Yij，red=Gij，red+jBYij，red=Yij，red∠φij] （15）

式中[φij]表示簡化后網絡中的導納陣[Yred]中第[i]行，第[j]列的元素所對應阻抗的相角。因此，簡化后系統中任意一臺同步發電機組輸出的功率可表示為：

[PEi=E′i2Gii，red+E′ij∈I，j≠iE′jYii，redsinδij+φij] （16）

式中：[δij]為簡化后系統中節點[i]與節點[j]間短路故障后電位的相角差;[E′i2Gii，red]為簡化后網絡中節點[i]處的負荷。

在這里假設簡化后的網絡是無損的，即[Gii，red=]0，[i≠j]，則式（16）變為：

[PEi=E′i2Gii，red+E′ij∈I，j≠iE′jYij，redsinδi-δj] （17）

將式（17）代入同步發電機轉子運動方程中，因此對于收縮網絡有：

[Miδi+Diδi=PTi-E′i2Gii，red-E′i?j∈I，j≠iE′jYij，redsinδi-δj] （18）

令[Pi=PTi-E′j2Gii，red]表示簡化后網絡中節點[i]的注入功率，[Pmax=E′iE′jYii，red]表示簡化后的網絡中節點[i，j]間能夠傳輸功率的最大值，因此可以將式（18）寫為：

[Miδi+Diδi=Pi-j∈I，j≠iPmaxijsinδi-δj] （19）

以風電為代表的清潔能源在接入電網時，通過電力電子設備與電網進行功率交換，與傳統同步機相比缺少對系統慣性的支撐[9]，因此當考慮接入點為風電時的振子模型采用Kuramoto一階模型表示：

[Diδi=Pi-j∈I，j≠iPmaxijsinδi-δj] （20）

式中[δ]為電力電子設備的相角。

3 ?電網的動態耦合強度聚類模型及算法

通過2.2節分析可知，在類Kuramoto電力網絡模型中，各節點之間的耦合強度可由兩節點間傳輸的最大功率來表示。當系統發生故障后，潮流分布發生轉移，節點之間的耦合強度也隨之改變，網絡中各節點的穩定狀態也發生改變，此時網絡有可能存在多種同步狀態。因此這里通過引入動態耦合強度分析各節點的不同同步狀態，更快更合理地將電網中節點分群。

令式（7）中[K=αK′i]，其中[K′i]是耦合強度，[α]是耦合權重，其表達式為：

[α=ηNj=1Ncos（θi-θj）] （21）

式中[η]稱為可塑性參數，由式（21）可知，耦合權重[α]的大小取決于振子相對時序。文獻[10]對耦合振子相互作用強度與振子之間的相互作用的相對時間的動態特性進行了深入研究，此處不再贅述。

圖3表示序參量與振子同步模式之間的關系，從圖中可以看出，當系統發生擾動時如果沒有足夠大的耦合強度[K]，則系統將出現從圖3a）～3c）逐漸失去同步過程的情況。當[K

電網的動態耦合強度的聚類算法流程和步驟如圖4所示。

1） 設置系統故障，記錄故障前后各發電機轉子角曲線。

2） 計算所有發電機轉子角在時間序列上的標準差。

3） 計算不同聚類的數量下標準差和步長的最小值和最大值。

4） 根據步長，按轉子角同步狀態將發電機聚類。

4 ?算例分析

IEEE 39節點系統的網絡拓撲結構如圖5所示，該系統包含10個發電機節點、19個負荷節點以及46條輸電線路，其中3號和6號發電機節點由相應容量的風電場替代。

設置故障場景1參數為：仿真時間為5 s，在0.4 s時節點14發生三相短路故障，經過0.02 s故障清除。網絡全局序參量[R]隨時間變化的曲線如圖6a）所示。

在整個故障發生前系統全局序參量[R]等于1，系統同步狀態近似于圖3a）。系統發生故障后，系統內振子同步狀態發生了改變，由于故障持續時間短，因此沒有出現圖3c）中失去同步的情況，并最終恢復同步，此時序參量[R]接近于1，系統內振子同步狀態如圖6b）所示。圖7為系統內各發電機轉子角隨時間變化曲線，從圖中可知故障清除后系統內4，5，7號發電機組轉子角發生了小幅震蕩，但最終恢復到初始值。

設置故障場景2參數為：仿真時間為5 s，在0.4 s時節點14發生三相短路故障，風電場從系統中切除，0.6 s后故障清除。故障清除后將系統發電機按照同步狀態聚類成三個機群。其中編號為2，4，5的發電機組為機群1，編號為8，9，10的發電機組為機群2。7號發電機為機群3。網絡中振子聚類后的各機群的序參量[R1，R2，R3]隨時間變化的曲線如圖8a）所示。聚類后各機群本身是同步的，因此[R1，R2]的數值接近于1，由于機群3中只有7號發電機一臺機組，因此[R3]的值等于1。系統內振子同步狀態如圖8b）所示。

圖9為系統內各發電機轉子角隨時間變化曲線，從圖中可知故障清除后，系統內振子失去同步，其中2，4，5號發電機組轉子角波動趨勢相近，8，9，10號發電機轉子角波動趨勢相近，7號發電機由于距離故障節點電氣距離較近，其轉子角波動變化較大，進而驗證了聚類算法的正確性。

5 ?結 ?論

本文針對含有風電電力系統的同步問題提出基于Kuramoto模型的網絡自適應耦合聚類同步方法，通過對聚類后各機群序參數的分析可以更清晰地了解網絡的同步狀態，避免完全崩潰。

參考文獻

[1] 姜惠蘭，吳玉璋，周照清，等.含雙饋風力發電場的多機系統暫態功角穩定性分析方法[J].中國電機工程學報，2018，38（4）：999?1005.

JIANG Huilan， WU Yuzhang， ZHOU Zhaoqing， et al. A method to analyze the transient angle stability of multi?machine system with DFIG?based wind farm [J]. Proceedings of the CSEE， 2018， 38（4）： 999?1005.

[2] 蘇福，楊松浩，王懷遠，等.電力系統暫態穩定時域仿真快速終止算法研究[J].中國電機工程學報，2017，37（15）：4372?4378.

SU Fu， YANG Songhao， WANG Huaiyuan， et al. Study on fast termination algorithm of time?domain simulation for power system transient stability [J]. Proceedings of the CSEE， 2017， 37（15）： 4372?4378.

[3] 陳厚合，王長江，姜濤，等.基于投影能量函數和Pin?SVM的電力系統暫態穩定評估[J].電工技術學報，2017，32（11）：67?76.

CHEN Houhe， WANG Changjing， JIANG Tao， et al. Transient stability assessment in bulk power grid using projection energy function and support vector machine with pinball loss [J]. Transactions of China electrotechnical society， 2017， 32（11）： 67?76.

[4] 王麗娟，蔡曉東，楊超，等.基于優化結構洞的無向加權網絡關鍵節點發現方法[J].現代電子技術，2017，40（6）：35?39.

WANG Lijuan， CAI Xiaodong， YANG Chao， et al. Method to identify key nodes in undirected weighted networks on basis of optimized ″structural hole″ [J]. Modern electronics technique， 2017， 40（6）： 35?39.

[5] DORFLER F， SIMPSON?PORCO J W， BULLO F， et al. Networks and algebraic graph theory： models， properties， and applications [J]. Proceedings of the IEEE， 2018， 106（5）： 977?1005.

[6] 陸君安，劉慧，陳娟.復雜動態網絡的同步[M].北京：高等教育出版社，2016：50?79.

LU Junan， LIU Hui， CHEN Juan. Synchronization in complex dynamical networks [M]. Beijing： Higher Education Press， 2016： 50?79.

[7] LI B， WONG K Y M. Optimizing synchronization stability of the Kuramoto model in complex networks and power grids [J]. Physical review E， 2017， 95（1）： 1?7.

[8] 吳成業，劉光曄，孫瑞，等.應用潮流方程簡化網絡計算的電力系統小干擾穩定分析[J].電力系統及其自動化學報，2018（6）：10?15.

WU Chengye， LIU Guangye， SUN Rui， et al. Small signal stability analysis of power system with simplified network compu?ting using power flow equations [J]. Proceedings of the CSU?EPSA， 2018（6）： 10?15.

[9] 黃李，鄒艷麗，王意，等.分布式電站的3種入網方式比較研究[J].廣西師范大學學報（自然科學版），2017，35（3）：30?36.

HUANG Li， ZOU Yanli， WANG Yi， et al. A comparative study on three types of distributed power plant connection stra?tegies [J]. Journal of Guangxi Normal University （natural science edition）， 2017， 35（3）： 30?36.

[10] CHANDRASEKAR V， SHEEBA J H， SUBASH B， et al. Adaptive coupling induced multi?stable states in complex networks [J]. Physica D： nonlinear phenomena， 2014， 267： 36?48.