一種用于電力系統(tǒng)電壓穩(wěn)定分析的雅可比矩陣關(guān)鍵特征值算法

周陶宏，李宏仲，鄭 健

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

0 引言

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

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

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

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

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

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

由于：

得到：

寫成矩陣形式：

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

對(duì)于電力系統(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ù)荷增長(zhǎng)控制參數(shù)），因此，系統(tǒng)雅可比矩陣則可寫為：

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

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

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

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

對(duì)兩邊求導(dǎo)：

整理并寫成矩陣形式：

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

定義：

式中：

可以拓展為：

式中：

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

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

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

關(guān)鍵特征值指得是隨著分岔控制參數(shù)μ的變化，實(shí)部最先變?yōu)榱愕囊粚?duì)共軛特征值：

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

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

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

得到時(shí)間：

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

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

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

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

3 算法步驟

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

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

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

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

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

（4）形成計(jì)算式（3）所需的矩陣；

（5）完成一步龍格－庫塔計(jì)算；

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

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

（8）確定步長(zhǎng)hi；

（9）完成一步龍格－庫塔計(jì)算；

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

4 算例分析

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

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

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

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

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

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

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

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

5 結(jié)語

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