周陶宏，李宏仲，鄭 健

（1.上海市電力公司，上海 200122；2.上海電力學(xué)院電力與自動化工程學(xué)院，上海 200090）

0 引言

采用霍普夫（Hopf）分岔理論分析電力系統(tǒng)電壓的穩(wěn)定性，必須準(zhǔn)確求解霍普夫分岔點的位置，特別是確定最先發(fā)生穿越虛軸現(xiàn)象的共軛特征值，即關(guān)鍵特征值。

George分別利用牛頓法、冪法、反冪法和Rayleigh商迭代法來計算占主導(dǎo)地位的關(guān)鍵特征值，并對這些計算方法的魯棒性和計算效率分別進行了對比［1－2］；文獻［3］利用改進的矩陣變換法來求解大規(guī)模系統(tǒng)動態(tài)模型的關(guān)鍵特征值；文獻［4］則進一步提出了利用基于多處理器的并行算法來提高計算效率；L.Wang等人提出了充分利用增廣矩陣稀疏特性的計算方法［5］；文獻［6］和文獻［7］分別提出了一種對矩陣的特征值進行連續(xù)追蹤的方法，Xiaoyu Wen等人將該方法引入了電力系統(tǒng)的穩(wěn)定研究中，提出了一種追蹤關(guān)鍵特征值的連續(xù)法［8］。

本文將在這些計算方法的基礎(chǔ)上做進一步的分析與研究。首先，利用特征值實部、特征值關(guān)于分岔參數(shù)的一階和二階靈敏度系數(shù)，構(gòu)造一個能夠判斷關(guān)鍵特征值的指標(biāo)，然后直接利用連續(xù)性方法，對關(guān)鍵特征值進行連續(xù)追蹤，直至霍普夫分岔點。

1 特征值連續(xù)追蹤的基本步驟

對于某一不對稱矩陣J（μ），以及該矩陣的某一特征值λ（μ），假設(shè)J（μ）和λ（μ）均為某一參數(shù)μ的連續(xù)函數(shù)，令w和v分別為λ對應(yīng)的左特征向量和右特征向量，則：

等式兩邊對參數(shù)μ求導(dǎo)并整理：

由于：

得到：

寫成矩陣形式：

將μ視為變化因子，結(jié)合龍格－庫塔法就可以對特征值λ進行追蹤，并且可以得到λ關(guān)于μ的變化曲線。

對于電力系統(tǒng)的DAE模型：

經(jīng)過線性化處理后：

式中：A＝Dxf（x，y）；B＝Dyf（x，y）；C＝Dxg（x，y）；D＝Dyg（x，y）。

由于矩陣A、B、C、D可以視為參數(shù)μ的函數(shù)（這里μ為分岔控制參數(shù)，即負(fù)荷增長控制參數(shù)），因此，系統(tǒng)雅可比矩陣則可寫為：

顯然，矩陣J基本上失去了稀疏性，因此本文將給出直接利用稀疏矩陣A、B、C、D追蹤λ變化曲線的表達(dá)式。

式中：wy和vy分別是對應(yīng)特征向量的拓展部分。

對于某一特征值λ，對應(yīng)的左右特征向量分別表示為w、v：

套用式（1）的形式得到：

對兩邊求導(dǎo)：

整理并寫成矩陣形式：

如果狀態(tài)變量個數(shù)為m，代數(shù)變量個數(shù)為n，則有向量：

定義：

式中：

可以拓展為：

式中：

式中：下標(biāo)ij表示該元素為矩陣的第i行、第j列元素。

從理論上講，利用式（2）可以對感興趣的任意一個特征值從初始狀態(tài)開始進行連續(xù)跟蹤，但是為了提高精確度，將利用連續(xù)性方法對系統(tǒng)平衡解流形進行連續(xù)追蹤。當(dāng)平衡解流形比較靠近霍普夫分岔時，再轉(zhuǎn)入對關(guān)鍵特征值的連續(xù)追蹤。

2 關(guān)鍵特征值的選擇

關(guān)鍵特征值指得是隨著分岔控制參數(shù)μ的變化，實部最先變?yōu)榱愕囊粚曹椞卣髦担?/p>

當(dāng)關(guān)鍵特征值的實部穿越虛軸時，意味著系統(tǒng)運行到了霍普夫分岔點，此時系統(tǒng)的平衡解不能保持穩(wěn)定，會出現(xiàn)振蕩性現(xiàn)象。而在利用連續(xù)法求解霍普夫分岔點的追蹤過程中，對于系統(tǒng)的每個平衡解都可以求解出系統(tǒng)雅可比矩陣的全部特征值中實部比較靠近虛軸的部分特征值。但在準(zhǔn)確追蹤到霍普夫分岔點之前，無法準(zhǔn)確地確定到底哪一個特征值是最終引發(fā)霍普夫分岔的關(guān)鍵特征值，因此需要在追蹤過程的每一步計算中根據(jù)一定的判斷標(biāo)準(zhǔn)，從求得的實部絕對值較小（即比較靠近虛軸）的特征值λ（μ）中再選取一小部分，做為準(zhǔn)關(guān)鍵特征值進行深入分析比較。

特征值的實部η（μ）可以做為一種判斷標(biāo)準(zhǔn)，實部越靠近虛軸說明其出現(xiàn)霍普夫分岔的可能性越大。但是隨著負(fù)荷的增長，在前一個平衡點相對離虛軸較遠(yuǎn)的特征值，可能在下一個平衡點距離虛軸就會變得比較近，也就是說隨著負(fù)荷的增長，特征值實部的變化率是各不相同的。由此，文獻［8］提出利用特征值實部對實部關(guān)于負(fù)荷增長的一階靈敏度系數(shù)的比值η（μ）／（?η（μ）／?μ）做為判斷標(biāo)準(zhǔn)。即在判斷關(guān)鍵特征值時考慮特征值實部的變化情況，將特征值實部較小而且變化幅度又比較大的特征值視為關(guān)鍵特征值。而文獻［9］指出，在平衡點接近霍普夫分岔點的時候，隨著負(fù)荷的增長，關(guān)鍵特征值的實部大小可能會發(fā)生突變，此時僅通過判斷一階靈敏度系數(shù)難以準(zhǔn)確地確定關(guān)鍵特征值。為此，本文提出利用特征值實部大小η（μ）、特征值實部關(guān)于負(fù)荷增長的一階靈敏度系數(shù)?η／?μ和二階靈敏度系數(shù)?2η／?μ2來綜合判斷關(guān)鍵特征值，同時參考速度、加速度與位移的關(guān)系表達(dá)式給出評判標(biāo)準(zhǔn)的計算表達(dá)式。一階和二階靈敏度系數(shù)的計算方法詳見文獻［8］和文獻［9］。

對于勻加速度運動，假設(shè)物體的位移為S，初始速度為v0，加速度為a，物體完成位移S所需的時間為t，則有：

得到時間：

如果將特征值實部絕對值｜η（μ）｜視為“位移”，特征值實部關(guān)于負(fù)荷增長的一階靈敏度系數(shù)?η／?μ視為“初始速度”，即隨著負(fù)荷做單位增長，特征值實部的相應(yīng)變化量。如果?η／?μ大于0，說明隨著負(fù)荷的增長，特征值實部η（μ）逐漸增大，假設(shè)系統(tǒng)未達(dá)到霍普夫分岔點之前，系統(tǒng)的平衡點為穩(wěn)定的平衡點，此時所有特征值實部應(yīng)該小于零。那么?η／?μ大于0就意味著隨著分岔控制參數(shù)μ的增大，該特征值實部是趨近于虛軸的，反之則該特征值遠(yuǎn)離虛軸。

進一步的將特征值實部關(guān)于負(fù)荷增長的二階靈敏度系數(shù)?2η／?μ2視為“加速度”，那么如果?2η／?μ2大于0就意味著隨著分岔控制參數(shù)μ的增大，?η／?μ是逐步增大的，也就是說隨著分岔控制參數(shù)μ的增大，該特征值實部靠近虛軸的趨勢在逐步加大；換言之，當(dāng)負(fù)荷增長時該特征值的實部接近虛軸的“速度”在逐步增大。由此，可以構(gòu)造判斷關(guān)鍵特征值的指標(biāo)如下：

很顯然，當(dāng)特征值實部η（μ）絕對值較小、特征值實部關(guān)于負(fù)荷增長的一階靈敏度系數(shù)?η／?μ和二階靈敏度系數(shù)?2η／?μ2都比較大的情況下，Icrt的值比較小。這就意味著該特征值由于本身就比較接近虛軸，而且隨著負(fù)荷的增長實部絕對值減小的趨勢比較快，因此該特征值最有可能發(fā)生穿越虛軸的情況。利用該指標(biāo)可以判斷哪一個特征值是關(guān)鍵特征值，在每一步計算中所得的特征值中，Icrt值最小的即視為關(guān)鍵特征值（對應(yīng)的指標(biāo)即Icrtmin）。由于Icrt的計算過程中用到了特征值實部η（μ）、一階靈敏度系數(shù)?η／?μ和二階靈敏度系數(shù)?2η／?μ2的大小，從而充分的計及了特征值當(dāng)前狀態(tài)和特征值隨著負(fù)荷增長的變化趨勢，因此其指示性較好，能夠?qū)μ卣髦当平撦S的態(tài)勢給出定量的分析。

為了提高計算的可靠性，根據(jù)追蹤平衡解流形的最后一步計算結(jié)果，對特征值按照Icrt值由小到大的順序進行排序，選擇排在前面的幾個特征值分別進行連續(xù)追蹤。

3 算法步驟

為了求解結(jié)果必須先確定追蹤步長。在計算過程中假設(shè)初始步長為h0，首先計算初始狀態(tài)下系統(tǒng)的雅可比矩陣的關(guān)鍵特征值及一階靈敏度系數(shù)?η／?μ和二階靈敏度系數(shù)?2η／?μ2，并計算Icrt，min，0，即關(guān)鍵特征值對應(yīng)的指標(biāo)Icrt值；然后執(zhí)行龍格庫塔法完成一步計算；在第二步計算過程中，首先求出關(guān)鍵特征值一階和二階靈敏度系數(shù)，最后計算出Icrt，min，1。步長h1采用下式確定：

追蹤關(guān)鍵特征值的混合方法的計算步驟：

（1）利用常規(guī)連續(xù)性方法追蹤系統(tǒng)的平衡解流形，直至平衡解流形比較靠近霍普夫分岔時（Icrt，min的值小于某一預(yù)先設(shè)定的值）；

（2）根據(jù)研究的需要利用Icrt的大小確定進行連續(xù)跟蹤的特征值；

（3）設(shè)定一個極小值ε用于判定是否出現(xiàn)特征值實部穿越虛軸的現(xiàn)象；

（4）形成計算式（3）所需的矩陣；

（5）完成一步龍格－庫塔計算；

（6）判斷所追蹤特征值的實部是否小于ε，如果小于則停止計算，否則轉(zhuǎn)到步驟（7）；

（7）求出?η／?μ和?2η／?μ2，并進一步求出I crt，min，1；

（8）確定步長hi；

（9）完成一步龍格－庫塔計算；

（10）修正計算所需的矩陣，轉(zhuǎn)到步驟（6）。

4 算例分析

以New-England 39節(jié)點系統(tǒng)為例，驗證所提方法的有效性。先從初始狀態(tài)開始利用常規(guī)連續(xù)性算法追蹤平衡解流形，當(dāng)Icrt，min小于0.15時，轉(zhuǎn)入對關(guān)鍵特征值的連續(xù)追蹤。表1給出了追蹤特征值的初始狀態(tài)。選表中Icrt值最小的一對共軛特征值進行連續(xù)追蹤，h0選0.018。

表1 部分特征值靈敏度系數(shù)

圖1給出了所追蹤的共軛特征值實部的變化，隨著負(fù)荷控制參數(shù)變化，特征值λ1，2在總負(fù)荷增量為791 MW左右時出現(xiàn)穿越虛軸的情況，從而出現(xiàn)霍普夫分岔，此時系統(tǒng)將可能出現(xiàn)振蕩性失穩(wěn)現(xiàn)象。

從理論上講，利用本文方法可以對任何感興趣的特征值進行追蹤，但是只有通過追蹤關(guān)鍵特征值才能確定霍普夫分岔點的位置。為了進一步說明這一問題，同樣對特征值λ3，4和特征值λ5，6進行追蹤，圖2和圖3給出了這兩對共軛特征值實部的變化曲線。

圖1 特征值λ1，2的η（μ）變化曲線圖

圖2 特征值λ3，4的η（μ）變化曲線圖

圖3 特征值λ5，6的η（μ）變化曲線圖

結(jié)果表明，當(dāng)這兩對共軛特征值實部為零時，系統(tǒng)的負(fù)荷增量大于791 MW。即這兩對共軛特征值不是最先穿越虛軸的關(guān)鍵特征值。同時也說明了指標(biāo)Icrt具有較高的指示性，能夠比較準(zhǔn)確的確定出關(guān)鍵特征值。通過追蹤該特征值得到的霍普夫分岔點是準(zhǔn)確的。

5 結(jié)語

提出了一種追蹤關(guān)鍵特征值的混合方法。該方法將追蹤平衡解流形的連續(xù)方法同關(guān)鍵特征值連續(xù)追蹤方法結(jié)合起來，先從初始狀態(tài)開始追蹤平衡解流形，當(dāng)平衡解靠近霍普夫分岔點時，利用判斷關(guān)鍵特征值的指標(biāo)Icrt來確定最有可能發(fā)生穿越虛軸現(xiàn)象的關(guān)鍵特征值并對其進行連續(xù)跟蹤，由于Icrt的計算過程中計及了特征值當(dāng)前狀態(tài)和特征值隨著負(fù)荷增長的變化趨勢，能夠較為準(zhǔn)確的判定關(guān)鍵特征值，并可利用Icrt來輔助確定跟蹤步長的大小。利用該方法對New England 39節(jié)點系統(tǒng)進行計算，得到了引起霍普夫分岔現(xiàn)象的關(guān)鍵特征值的具體變化曲線。

［1］ Semlyen A，Wang L.Sequential Computation of the Complete Eigensystem for the Study Zone in Small Signal Stability Analysis of Large Power Systems［J］.IEEE Trans.on Power System，1988，3（2）：715-725.

［2］ George Angelidis，Adam Semlyen.Improved Methodologies for the Calculation of Critical Eigenvalues in Small Signal Stability Analysis［J］.IEEE Trans.on Power Systems，1996，11（3）：1209-1217.

［3］ Leonardo T G Lima，Licio H Bezerra，Carlos Tomei，et al.New Method for Fast Small-Signal Stability Assessment of Large Scale Power Systems［J］.IEEE Transaction on Power System，1995，10（4）：1979-1985.

［4］ Jorge M Campagnolo，Nelson Martins，Jose L R Pereira，et al.Fast Small-signal Stability Assessment using Parallel Processing［J］.IEEE Transaction on Power System，1993，9（2）：949-956.

［5］ L Wang，A Semlyen.Application of Sparse Eigenvalue Techniques to the Small Signal Stability Analysis of Large Power Systems［J］.IEEE Trans.on Power Systems，1990，5（2）：635-642.

［6］ Robert Kalaba，Karl Spingarn，Leigh Tesfatsion.Variational Equations for the Eigenvalues and Eigenvectors of Nonsymmetric Matrices［J］.Journal of Optimaization Theory and Applications，1981，33（1）：1-8.

［7］ Ji W，Venkatasubramanian V.Dynamics of a Minimal Power System：Invariant Tori and Quasiperiodic Motions［J］.IEEE Trans.on Circuits and Systems I：Fundamental Theory and Applications，1995，42（12）：981-1000.

［8］ Xiaoyu Wen，Ajjarapu V.Application of a Novel Eigenvalue Trajectory Tracing Method to Identify Both Oscillatory Stability Margin and Damping Margin［J］.IEEE Trans.on Power System，2006，21（2）：817-824

［9］ Hae-Kon Nam，Yong-Ku Kim，Kwan-Shik Shim，et al.A New Eigen-sensitivity Theory of Augmented Matrix and its Applications to Power System Stability Analysis［J］.IEEE Trans.on Power Systems，2000，15（1）：363-369