提高電力系統潮流計算收斂性方法的比較分析

范翔，嚴正，趙文愷，范金燕，王虹富

（1．上海交通大學電力傳輸與功率變換控制教育部重點實驗室，上海 200240；2.上海交通大學應用數學系，上海 200240；3．中國電力科學研究院，北京 100192）

隨著輸電規模的擴大，電網將形成互聯大區域的復雜結構，系統運行工況變得多樣化。電網方式安排人員在得到合理的系統運行工況前，會遇到許多潮流難以收斂的情況，而此時傳統牛頓法計算失敗，難以為后續的穩定分析及方式調整提供有效信息，因此，研究提高潮流計算收斂性的方法具有普遍的現實意義和重要的實用價值。

造成傳統牛頓法或者PQ解耦法潮流計算難以收斂的主要原因是由于迭代過程中的雅可比矩陣奇異或接近奇異[1]。目前能夠用于提高大規模電力系統潮流計算收斂性的典型方法主要有：最優乘子法[1-5]、張量法[6-9]和自適應LM方法[10]。文獻[1]提出最優乘子法計算潮流，具有永不發散的特點，且與現有的潮流程序接口簡單、易于實現，是一種通用提高潮流收斂性的方法，在各類難收斂潮流計算中效果較好；文獻[6-7]將張量法引入潮流計算，考慮了潮流方程泰勒展開近似二階項的影響，在部分區域重負荷情況下，張量法能夠取得到較好解；文獻[10]引入自適應LM方法計算潮流方程，該方法在迭代過程中雅可比矩陣始終非奇異，通過算例仿真驗證了自適于應LM方法，減小了對初值選取的要求，擴大了收斂范圍；文獻[11]基于優化技術的潮流方法[11]，將求解潮流方程組轉化成計算優化模型，使得問題復雜化。基于大范圍收斂同倫法提高潮流計算收斂性，擴大了潮流方程組的收斂范圍，但計算量大，未引入商業應用[12-14]。

本文研究在常規潮流基礎上比較容易實現的3類方法：張量法、最優乘子法和自適應LM方法，及其數值過程。總結了各方法進行電力系統潮流計算的特點。基于數值過程，介紹了3種方法關鍵步的稀疏實現，比較了3種方法相對牛頓法單步迭代的額外計算開銷。在1個標準系統和2個實際系統中進行仿真，仿真計算結果表明，相比其他2種方法，自適應LM方法應用于實際大系統需結合稀疏技術。

1 提高潮流收斂性的3種方法

令潮流方程為

式中，狀態變量x=[V，θ]T。當潮流方程雅可比矩陣J（x）奇異時，常規牛頓法失效，需采用提高潮流收斂性的方法。

1.1 最優乘子法

所謂最優乘子法，即在每一步計算狀態變量的修正量Δx之后，不直接用Δx去修正狀態變量x，而是乘以一個標量乘子λ去修正，即

最優乘子法的核心在于乘子λ的計算，思想是在牛頓迭代方向上尋求最優步長，即

式中：λk為第k步迭代時的最優乘子；dk為第k次迭代的牛頓步。若在當前牛頓方向上搜尋不到最優步長，則λ→0，方法停止。若潮流結果收斂，則d→0，當方法收斂時λ→1。由此可見，引入最優乘子可有效避免因迭代步步幅過大而導致方法振蕩或發散。最優乘子方法的具體計算流程參見文獻[1-2]。

1.2 張量法

張量法潮流計算的實質為含二階項潮流計算。通過對潮流方程的二階展開項近似，計算出張量步，修正牛頓迭代步。

1.2.1 插值張量法

在當前迭代點xk處將式左側展開至二階，即

若潮流方程采用直角坐標形式，則張量Ak為常值三維矩陣。由文獻[12]知，Ak可用前述迭代點表示為其中，sk=xi-xk，i∈[1，k-1]且使得{sk}線性無關；ak為已知向量。將Ak帶入式（4）計算迭代步dk，得方程

（i=1，…，p）的二次方程組，即

式（6）稱為張量方程組，若僅取前一次迭代點xk-1做插值，則式（6）僅有一個二次方程和一個待求參數β，回代入式（5），得

若式（7）不存在實數解，則采用正交變換將關于迭代步dk的最小二乘min‖F（xk+dk）‖2轉換為關于張量參數β的最小二乘，計算出dk做為迭代步[10]，該方法依然要求雅可比矩陣非奇異。方法具體流程詳見文獻[10]。

1.2.2 直接張量法

設實際迭代步為d=dN+dT，其中dN為牛頓步，dT為張量步，且參考直接張量法原始文獻[6]知，‖dN‖2?‖dT‖2。將迭代步帶入式（4），令其等于0，計算得到

1.2.3 數值特點

在良態迭代過程中，牛頓法模型準確，直接張量法的假設條件可以被滿足。當潮流方程難以收斂，按照含二階項潮流計算理論，牛頓法線性化模型的誤差較大，牛頓步不準確，需引入張量步進行修正。根據張量修正步dT的計算表達式知，在雅可比矩陣趨近奇異時dT的幅值較大從而‖d‖2也較大，也會出現與常規牛頓法相同的振蕩發散問題。且直接張量法缺乏收斂性理論證明。

兩種張量法提高收斂性的特點為張量步對牛頓步的修正使得迭代步d變為

易知張量法收斂的充分條件是‖F+C‖2為雅可比矩陣最小奇異值σmin的高階無窮小。由此可知，C所起補償作用是在雅可比矩陣條件數較大或接近奇異時減小部分潮流方程的失配量，防止迭代步幅過大而發散。

綜上，兩類張量法的特點是：①引入二階補償項對潮流偏差量進行補償，當迭代過程中雅可比矩陣接近奇異或條件數很大時，減小迭代步長；②基于插值的張量法由文獻[9]奠定了數學基礎，而直接張量法缺乏相關理論依據。

1.3 自適應LM方法

將式（1）左側在當前迭代點處一階展開為

式中，dk=xk+1-xk。

潮流方程的最小二乘模型為

當式（11）所得解x~滿足G（x~）=0時，x~即為潮流方程的解。

將式（10）代入式（11）并引入步長約束，得到最初計算自適應LM方法迭代步的模型，即

式中，參數δ按照一定方式更新[16]。由文獻[16]可得式（12）的解為

式中，μk為自適應LM方法阻尼因子，μk=αk‖F（xk）‖，其中αk為非負調節因子。自適應LM方法的基本計算流程[17]為

步驟1平啟動，迭代步數k置1，設置α1和常數m及收斂精度ε，使得α1＞m，0＜p0＜p1＜p2＜1。

步驟2計算μk，按式（13）計算dk。

步驟3計算τk

選擇是否接受dk

步驟4自適應因子的調整

步驟5用潮流收斂判據‖F（xk）‖2＜ε判別潮流收斂與否，不收斂則返回步驟2繼續迭代，至最大迭代次數方法停止輸出結果。

對比最優乘子法，當雅可比矩陣條件數很大或接近奇異時，牛頓步非下降方向，最優乘子強制為0，由于不能改變迭代方向，方法停止在潮流失配量較小的近似潮流解上。而式（13）中自適應LM方法通過改變阻尼因子μk來調整迭代步步幅和方向。因此在計算難收斂潮流工況時，自適應LM方法較最優乘子法有更好的適應性。另外，自適應LM方法在迭代過程中求逆矩陣JTJ+μI始終非奇異，因而自適應LM方法較具有更好的數值穩定性和收斂性。

綜上所述，自適應LM方法在潮流計算有兩個特點：①在計算過程中保持求逆矩陣JTJ+μI對稱正定，不會遇到奇異問題，因此自適應LM方法具有較好的數值穩定性；②自適應LM方法通過改變阻尼因子μ的值，來同時改變迭代步長和迭代方向，對難收斂潮流具有更好的適應性。同時方法可收斂至潮流方程的最小二乘解。

2 各方法的稀疏實現及計算量比較

2.1 稀疏實現

張量法和最優乘子法的求逆矩陣結構與牛頓法雅可比矩陣相同，故可采用與常規牛頓法相同的矩陣LU分解（LU decomposition/三角分解）及前代回代提高計算效率。因此，張量法和最優乘子法均可計算大規模系統潮流。

由式知計算自適應LM方法迭代步dk的求逆矩陣JTJ+μI結構不同于牛頓法雅可比矩陣J，且較J矩陣稀疏度有所下降，但其仍為高度稀疏的矩陣，因此不能直接采用牛頓法的稀疏處理。

JTJ+μI在結構和數值上均對稱，且為正定矩陣，可以采用稀疏Cholesky分解[18]，來求解稀疏線性方程組。由于矩陣JTJ+μI與系統導納矩陣稀疏結構不同，所以不能沿用基于導納矩陣的排序結果。此時采用近似最小度AMD（approx-imate minimum degree）排序算法[18]對JTJ+μI矩陣排序，減少Cholesky分解時的注入元個數。由于整個迭代過程中矩陣JTJ+μI的結構保持不變，故只需在首次迭代時對矩陣JTJ+μI進行排序即可。

2.2 計算量比較

基于插值張量法的每步迭代，當p=1時，求解式（6）時需求解一次大規模稀疏線性方程組，在求解式（7）時還需求解一次大規模線性方程組。由于單步迭代中求逆矩陣因子表固定，所以p=1時插值法的單步計算量約為普通牛頓法單步計算量的2倍。

對于直接張量法的每步迭代，計算dN和dT需各求解一次大規模稀疏線性方程組，因而單步計算量約為普通牛頓法的2倍。

最優乘子法僅在計算牛頓步時求解一次大規模稀疏線性方程組，因此最優乘子法單步迭代計算量約等于普通牛頓法。

自適應LM方法的單步迭代中，需要由式（13）來求解迭代步，而式（13）中的矩陣JTJ與狀態估計中的系數矩陣有類似的結構，因此可以認為自適應LM方法單步計算量比普通牛頓法略有增加，本文的算例仿真部分驗證了該結論。

綜合第1節的數值機理分析，3種方法的單步計算量和收斂特點如表1所示。

表1 最優乘子法、張量法和自適應LM方法單步計算量及收斂特點匯總Tab.1 Computational costs per step and convergent characteristic of optimal multiplier method，tensor method and self-adapted LM method

3 算例仿真

仿真平臺為主頻3.0GHzCore2 CPU，4G內存的PC機，操作系統為Win7，采用Matlab2013a和C++的混合編程方式，C++實現與BPA的數據接口，Matlab實現各方法以及測試各算例。算例中的具體信息如表2所示。

表2 測試系統數據Tab.2 Data of test system s

比較幾種方法的單步迭代時間，其中自適應LM方法調用Matlab函數chol和amd實現稀疏Cholesky分解及AMD排序。

算例1采用IEEE39節點系統構造重負荷難收斂潮流。

算例2采用來源于mat power4.1[19]波蘭3375wp測試系統，采用mat power自帶潮流求解器因迭代初始點雅可比矩陣奇異無法收斂。

算例3源于華東電網2004年夏高的BPA數據，加重系統部分區域負荷，至BPA潮流程序不收斂。模擬局部重負荷難收斂工況。

對以上3個算例分別采用張量法、最優乘子法和自適應LM方法計算各系統的潮流，驗證并比較3種方法的收斂特性。

算例4采用國網21 479節點數據，用于比較文獻[10]的LM方法和本文的自適應LM方法，在大規模系統潮流計算中的適應性。

對于重負荷系統，有可能出現PV-PQ節點的頻繁轉換，使得潮流計算的收斂域變窄[20]，在某種程度上惡化了潮流程序的收斂性。由于本文旨在對潮流計算的3種主要方法的收斂性和計算量進行分析和比較，因此，算例中不考慮系統中PV-PQ節點轉換的情況。

在比較各方法的收斂特性之前，先通過華東04系統來比較各方法的單步計算量，結果如表3所示。

表3 單步迭代計算時間比較Tab.3 Comparisonin of computing time per iteration between each method s

由表3可知，各方法的單步計算時間與常規牛頓法的單步時間之比與表1中給出的結論基本一致。同時，也驗證了用稀疏Cholesky分解和AMD排序能夠顯著提高自適應LM迭代步的計算效率。

3.1 算例1

采用IEEE39節點系統構造重負荷難收斂潮流。良態情況下各方法所得結果與BPA所得潮流結果比較，電壓相角（rad）和幅值（p.u.）的絕對偏差為10-7，證明所寫程序無誤。

將非發電機節點的有功負荷擴大1.3倍，無功負荷擴大1.4倍，BPA牛頓法計算該潮流不收斂。采用張量法、自適應LM方法和最優乘子法計算潮流方程組迭代曲線如圖1和圖2所示。

自適應LM方法參數取值為：α1=0.1，m=10-5，p0=10-4，p1=0.75，p2=0.95。各方法的收斂準則均為潮流偏差‖F‖2≤10-5。

圖1 兩類張量法迭代曲線Fig.1 Iteration curves of two tensor methods

從圖1可見，直接張量法在潮流偏差較大時陷入振蕩，插值張量法補償項取得的結果更好，但在得到10-2數量級的近似潮流解時，插值張量法陷入數值振蕩。兩種方法近似二階補償項C的補償效果插值張量法優于直接張量法。

從第15次迭代可以看出，插值張量法二階補償項Cp的修正使得潮流偏差量減小為最小奇異值的高階無窮小。而直接張量法二階補償項Cd的補償性能不佳，導致步幅仍呈現增大趨向。第21次迭代時，補償項不能補償潮流偏差至合適值，此時最小奇異值為：0.03，因而產生振蕩。張量法雖在迭代中得到最小潮流偏差量近似解，但因補償項不能跟上雅可比矩陣奇異值減小的速度最終振蕩，難以通過現有的方法將迭代停止在合適的潮流解上。

圖2 最優乘子法和自適應LM方法迭代曲線Fig.2 Iteration curve of self-adaptive LM method and optimal multiplier method

從圖2可以看出，初始迭代段自適應LM方法由于潮流偏差較大阻尼因子很大，因而迭代步-（JTJ+μI）-1JTF呈現出最速下降步-1/μJTF折線下降的特點。迭代至第5步時兩種方法具有十分接近的潮流偏差量。第5、6步迭代時最優乘子為6.1×10-4和3.3×10-5；最優乘子法第5步迭代時雅可比矩陣的最小特征值為3.5×10-3，第6步迭代時雅可比矩陣的最小特征值為8.2×10-4；由此可見，雅可比矩陣趨近奇異，最優乘子趨近0，方法停滯。自適應LM方法在第5次迭代時求逆矩陣最小特征值為0.09，第6次迭代時求逆矩陣最小特征值為6.0×10-3，由此可驗證自適應LM方法阻尼因子繞開了雅可比矩陣在第5次迭代接近奇異的區域，獲得了潮流偏差量更小的解。

3.2 算例2

波蘭3375wp潮流數據取自波蘭2007年冬季晚高峰潮流斷面。應用matpower自帶潮流求解器計算該系統潮流發散。由于平啟動迭代時雅可比矩陣趨近奇異，最優乘子法停滯，張量法均失效。系統的雅可比矩陣趨近奇異，最小特征值為0，最小奇異值為0。

由于首次迭代雅可比矩陣趨近奇異，最優乘子法停滯，張量法失效，只有自適應LM方法順利計算出了該系統的潮流。自適應LM方法計算該系統潮流耗時1 s，迭代曲線如圖3所示。

圖3 自適應LM方法迭代功率失配量曲線Fig.3 Variation of power mismatches via self-adaptive LM method

從圖3可以看出，經過10次迭代后，自適應LM方法平穩收斂至‖F‖2=7.2×10-8。

采用最優乘子法計算該系統時，首次迭代遇到雅可比矩陣趨近奇異，導致最優步長趨近0，方法停滯。

采用張量法計算時，由于雅可比矩陣趨近奇異，牛頓步無法計算，因此張量法失效。

采用自適應LM方法計算得到的電壓如圖4所示。

圖4 電壓幅值柱形圖Fig.4 Voltage amplitude histogram

從圖4可以看出，電壓幅值在允許范圍內。計算所得，最小相角為-37.29°，最大相角為3.21°。該潮流解未呈現出病態情況，表明自適應LM方法計算求得真實潮流解。

3.3 算例3

該算例采用2004年夏季華東電網實際潮流斷面數據。原始數據為良態潮流數據，將SH區域的有功負荷擴大21.98%，無功負荷擴大60%；JS區域的有功負荷擴大7.95%，無功負荷擴大6.82%，以此模擬部分區域重負荷難收斂工況。此時BPA潮流程序計算發散，自適應LM方法、最優乘子法所得迭代曲線如圖5所示。

圖5 自適應LM方法及最優乘子法迭代曲線Fig.5 Iterative curves of self-adaptive LM method and optimal multiplier method

采用插值張量法計算該系統潮流，獲得初始迭代點后，首次迭代雅可比矩陣奇異，由于初值非足夠接近潮流解故不能獲得牛頓步，即不能通過牛頓步求解張量修正步。采用直接張量法亦存在此問題。

取最優乘子法停滯時的雅可比矩陣Jop_end和功率偏差向量Fop_end，計算得到Jop_end的正交基Qop_end后發現同樣，取自適應LM方法終止時的雅可比矩陣JLM_end和功率偏差向量FLM_end，計算得到FLM_end的正交基QLM_end后發現，因此，自適應LM方法得到了該系統潮流的精確最小二乘解。

3.4 算例4

本算例旨在比較文獻[10]的方法與本文的自適應LM方法的收斂性。由于迭代過程中雅可比矩陣出現奇異，最優乘子法和兩類張量法均失效。

[10]提出的LM方法，其阻尼因子為μ=0.001‖F‖2與功率偏差量呈線性關系。即潮流偏差大時阻尼因子大，潮流偏差小時阻尼因子小，具備提高潮流收斂性的LM方法阻尼因子的特點。但是該阻尼因子在當前迭代步不能調節，即不能獲得當前迭代點的最優下降方向，從而增加迭代次數，計算速度慢。

本文的阻尼因子采用

式中，α為自適應因子，通過類似信賴域的方法進行調節，阻尼因子與功率偏差量呈非線性關系，能夠獲得當前迭代點合適的迭代方向和迭代步幅。通過國網21479節點算例可驗證（見圖6）。

圖6 2種方法計算國網21479系統迭代曲線比較Fig.6 Comparison of two methods in terms of iterations for state grid-21479 system

從圖6可以看出，采用自適應因子調節的阻尼因子更能適應大規模系統潮流計算。

4 結論

本文分析了張量法、最優乘子法、和自適應LM方法的數值特點，以此為基礎比較了3種方法的收斂特點及單步計算量。并對1個標準系統以及2個實際系統進行仿真計算。

（1）在重負荷情況下，插值張量法利用二階項能夠較好地補償功率偏差，對難收斂潮流具有較好地適應性；然而當雅可比矩陣趨近奇異時，易振蕩；直接張量法補償項效果不佳。

（2）當系統潮流呈現重負荷而難收斂時，最優乘子法最優乘子迅速減小至0，收斂到近似最小二乘，近似解處偏差量仍較大；當初始點難以確定甚至雅可比矩陣奇異時，最優乘子法失效。

（3）采用自適應LM方法提高潮流計算收斂性較最優乘子法和張量法有更好的魯棒性和數值穩定性。由于阻尼因子的存在，自適應LM方法迭代過程中不會發生因雅可比矩陣奇異而導致的發散問題。自適應LM方法較小依賴初值，對初值難設定、重負荷等工況具有很好的適應性。

參考文獻：

[1]Iwamoto S，Tamura Y. A load flow calculation method for ill-conditioned power systems[J]. IEEE Trans on Power Apparatus and Systems，1981，100（4）：1736-1743.

[2]王憲榮，包麗明（Wang Xianrong，Bao Liming）.極坐標系準最優乘子病態潮流解法研究（The study of ill-conditioned load flow using quasi-optimal factor in polar coordinates）[J].中國電機工程學報（Proceedings of the CSEE），1994，14（1）：40-45.

[3]Overbye T J.Computation of a practical method to restore power flow solvability[J].IEEE Trans on Power Systems，1995，10（1）：280-287.

[4]石飛，趙晉泉，王毅（ShiFei，Zhao Jinquan，Wang Yi）.計及發電機無功約束的最優乘子潮流計算方法比較（A comparison of Newton -Raphson optimal multipliers load flow methods considering generator reactive power limitation）[J].電力系統保護與控制（Power System Protection and Control），2009，37（2）：6-10.

[5]胡澤春，嚴正（Hu Zechun，Yan Zheng）.帶最優乘子牛頓法在交直流系統潮流計算中的應用（Application of Newton load flow methods with optimal multiplier for AC/DC power system）[J].電力系統自動化（Automation of Electric Power Systems），2009，33（9）：26-31.

[6]Salgado R S，Zeitune A F.Power flow solutions through tensor methods[J].IET Generation，Transmission&Distribution，2009，3（5）：413-424.

[7]林濟鏗，吳鵬，袁龍，等（Lin Jikeng，Wu Peng，Yuan Long，et al）.基于張量法的電力系統潮流計算（Power flow computation based on tensor methods）[J].中國電機工程學報（Proceedings of the CSEE），2011，31（34）：113-119.

[8]Bouaricha A，Schnabel R B.Tensor methods for large sparse systems of nonlinear equations[J].Mathematical Programming，1998，82（3）：377-400.

[9]Schnabel RB，Frank PD.Tensor methods for nonlinear equations[J].SIAM Journal on Numerical Anal-ysis，1984，21（5）：815-843.

[10]Lagace PJ.Power flow methods for improving convergence[C]//38th Annual Conference of IEEE Industrial Electronics.Montreal，Canada，2012：1387-1392.

[11]覃智君，侯云鶴，吳復立（Qin Zhijun，Hou Yunhe，Wu Fulix）.大規模交直流系統潮流計算的實用化模型（Practi-cal model for large-scale AC-DC system power flow calcu-lation）[J].中國電機工程學報（Proceedings of the CSEE），2011，31（10）：95-101.

[12]周佃民，廖培金（Zhou Dianmin，Liao Peijin）.電力系統病態潮流的同倫方法求解（Homotopy method for illconditioned power system load flow calculation）[J].電力系統及其自動化學報（Proceedings of the CSU-EPSA），1999，11（5-6）：67-71.

[13]陳禮義，戴宏偉，張琦鵬（Chen Liyi，DaiHongwei，Zhang Qipeng）.一種大范圍收斂的電力系統潮流算法——同倫延拓法（A large range convergent method for power system load calculation—Homotopy calculation method）[J].電力系統及其自動化學報（Proceedings of the CSUEPSA），1993，5（1）：67-74，30.

[14]孫小舟（Sun Xiaozhou）.使用“BPA電力系統分析程序”的體會（Experiences on using‘BPA power system analysis program’）[J].電力建設（Electric Power Construction），1998，19（11）：34-35.

[15]周國標，宋寶瑞，謝建利.數值計算[M].北京：高等教育出版社，2008.

[16]Dennis JE Jr.Numerical Methods for Unconstrained Optimization and Nonlinear Equations[M].北京：科學出版社，2009.

[17]Fan Jinyan，Yuan Yaxiang. On the quadratic convergence of the Levenberg-Marquardt method without nonsingularity assumption[J]. Computing，2005，74（1）：23-39.

[18]Davis T A.DirectMethods for Sparse Linear Systems[M].Philadelphia：SIAM，2006.

[19]Zimmerman R D，Murillo-Sanchez C E，Thomas R J.MATPOWER：Steady-state operations，planning，and analysis tools for power systems research and education[J].IEEE Transon Power Systems，2011，26（1）：12-19.

[20]趙晉泉，江曉東，張伯明（Zhao Jinquan，Jiang Xiaodong，Zhang Boming）.潮流計算中PV-PQ節點轉換邏輯的研究（Study on PV-PQ bus type switching logic in power flow computation）[J].中國電機工程學報（Proceedings of the CSEE），2005，25（1）：54-59.