王醫(yī)強,張莉萍 ,陳宇晨,范 蕊
(上海工程技術(shù)大學(xué) 電子電氣學(xué)院,上海 201620)
基于C#牛頓潮流算法的分析
王醫(yī)強,張莉萍 ,陳宇晨,范 蕊
(上海工程技術(shù)大學(xué) 電子電氣學(xué)院,上海 201620)
針對Matlab與C#混合開發(fā)潮流計算程序時平臺間存在通信不便的問題,需要通過分析牛頓拉夫遜潮流算法在C#上的運行效率和可靠性,驗證使用C#是否能更好的開發(fā)潮流計算程序,從而統(tǒng)一開發(fā)平臺。選擇在Visual Studio2013平臺上,使用C#語言實現(xiàn)牛頓拉夫遜法潮流計算,并且對C#潮流計算的結(jié)果和運行時間與Matlab潮流計算的結(jié)果和運行時間進行對比分析。求解IEEE4系統(tǒng)的電網(wǎng)數(shù)據(jù)得出C#語言運算時間比Matlab要少55 ms,所求節(jié)點數(shù)據(jù)與Matlab所求出的數(shù)據(jù)相比,相差精度在10-5數(shù)量級。通過實驗得出基于C#牛頓潮流算法可靠性高,運行效率高,可用于開發(fā)潮流計算程序。
導(dǎo)納矩陣;雅可比矩陣;關(guān)聯(lián)矩陣;牛頓法
潮流計算中常用的算法是牛頓拉夫遜法,分析該算法時大多是在Matlab軟件上。隨著當(dāng)前智能電網(wǎng)的不斷發(fā)展,開發(fā)軟件功能的不斷完善,數(shù)據(jù)庫技術(shù)的突破,對潮流計算的要求也有所變化,以前的Matlab平臺無法滿足在線潮流計算,計算速度也無法與編譯型語言相比。
文獻[1]中利用Matlab完成潮流計算算法,VB完成程序的界面,通過ActiveX技術(shù)實現(xiàn)了Matlab與VB間的數(shù)據(jù)通信[1]。由于Matlab和C#均是基于C語言的高級語言的特點,文獻[2]中利用該特點實現(xiàn)C#和Matlab的混合編程[2],C#語言用于實現(xiàn)界面的開發(fā),Matlab用于完成龐大的數(shù)值運算和圖像生成。文獻[3]中分析了Basic、Fortran和Matlab等高級語言在潮流計算中的各自特點[3],針對潮流計算數(shù)學(xué)模型,在Matlab上開發(fā)出了潮流算法。文獻[4]中介紹了一種基于GIS[4]的潮流算法,并在C#語言的基礎(chǔ)上開發(fā)與該算法相關(guān)的應(yīng)用程序。文獻[5]中從潮流計算的意義、目的、方法以及實驗等方面對潮流計算進行了分析[5]。文獻[6]中運用MatPower軟件進行潮流計算,并且把MatPower潮流計算的最終結(jié)果與Matlab和SimuLink[6]潮流計算的結(jié)果做比較。文獻[7]中針對牛頓類潮流計算的初值敏感性問題提出牛頓潮流計算的收斂性定理,給出牛頓類潮流收斂的充分條件,解決了用牛頓法進行潮流計算時因初值選取不當(dāng)使潮流收斂過慢甚至無法收斂的問題[7]。文獻[8]中在傳統(tǒng)Fast-ICA算法的基礎(chǔ)上,為了加快算法的收斂速度,提高算法的運行效率利用8階收斂的牛頓迭代方法對Fast-ICA算法進行優(yōu)化[8]。文獻[9]中利用Matlab研究改進牛頓算法在配電網(wǎng)潮流計算中的應(yīng)用[9]。文獻[10]中對直角坐標(biāo)和極坐標(biāo)牛頓法潮流計算的速度進行比較,在PV節(jié)點數(shù)較少的情況下直角坐標(biāo)牛頓法的潮流計算速度快于極坐標(biāo)牛頓法[10],在PV節(jié)點數(shù)較多的情況下則相反。文獻[11]中在一階牛頓法的基礎(chǔ)上利用一階牛頓法的修正量作為預(yù)測量,用傳統(tǒng)牛頓法進行迭代計算,使用Matlab編程仿真,實驗結(jié)果表明該方法迭代次數(shù)少,收斂性好[11]。
圖1 單變量函數(shù)牛頓迭代法[12]。
圖1 單變量函數(shù)牛頓迭代法
為了更好地引用牛頓法的數(shù)學(xué)原理,首先對單變量方程和相關(guān)公式分析。設(shè)上圖中單變量方程為
f(x)=0
(1)
當(dāng)圖中k=0時設(shè)x(0)為解的近似值,與真解的誤差為Δx(0),則
f(x(0)+Δx(0))=0
(2)
將上式左邊的函數(shù)在x(0)附近展開為泰勒級數(shù),由于Δx(0)很小,則Δx(0)的二次及以上的階次可略去,從而得到以下修正方程
f(x(0)+Δx(0))=f(x(0))+f′(x(0))Δx(0)=0
(3)
解該方程可得到修正量為
(4)
根據(jù)以上修正量可得迭代一次后近似解為
(5)
接下來將單變量的牛頓拉夫遜法應(yīng)用到求解多變量非線性代數(shù)方程。設(shè)多變量非線性代數(shù)方程為
(6)
將各變量的初值和各變量的修正量代入以上方程組得
(7)
將上式中的n個多元函數(shù)在初始值附近分別展成泰勒級數(shù),并略去含有修正量的二次及以上階次的各項,可得到以下方程組
(8)
根據(jù)方程組可轉(zhuǎn)化為矩陣形式
(9)
多變量非線性方程的迭代過程如單變量非線性方程一樣,因此可將第k次迭代的修正方程寫成如下矩陣形式,便于應(yīng)用到潮流計算中
F(X(k))=-J(k)ΔX(k)
X(k+1)=X(k)+ΔX(k)
(10)
以上J為雅可比矩陣,X為多變量列向量。
通過從單變量到多變量牛頓法數(shù)學(xué)原理的分析,能夠更好的應(yīng)用牛頓法。
基于以上牛頓法的數(shù)學(xué)原理分析,將牛頓法應(yīng)用到潮流計算,將電網(wǎng)中的節(jié)點數(shù)據(jù)和支路數(shù)據(jù)引入牛頓法中得到牛頓法潮流計算數(shù)學(xué)模型。由于在電網(wǎng)節(jié)點數(shù)比較少的情況下,直角坐標(biāo)表示的節(jié)點電壓方程比極坐標(biāo)表示的節(jié)點電壓方程收斂速度快,此次分析的節(jié)點電壓數(shù)比較少,因此潮流計算中的節(jié)點電壓用直角坐標(biāo)表示。
節(jié)點電壓可表示為
Ui=ei+jfi
(11)
導(dǎo)納矩陣為
Yij=Gij+jBij
(12)
根據(jù)電力系統(tǒng)實際運行條件,電網(wǎng)節(jié)點一般分為PQ節(jié)點、PU節(jié)點、平衡節(jié)點,平衡節(jié)點電壓相位一般為零,幅值給定,因此不參與到牛頓法迭代計算中。約束條件包括電壓約束、有功功率和無功功率約束以及相位差約束。潮流計算的數(shù)學(xué)模型是根據(jù)以下式(13)和式(14)、式(15)求得的步長進行逐步線性化,從而求得滿足電網(wǎng)運行要求的近似解
(13)
通過式(13)和導(dǎo)納矩陣的實部以及虛部數(shù)據(jù)可求得PQ節(jié)點的有功功率和無功功率。將式(13)代入式(14)中可求得有功功率和無功功率的迭代步長
(14)
對于PU節(jié)點時,可通過式(15)求得有功功率迭代步長和電壓迭代步長
(15)
將以上各個變量的步長放入到矩陣中,逐步線性化的過程可以用修正方程(16)表示,其中J為雅可比矩陣,ΔW里的變量為節(jié)點功率的步長,ΔU是每次迭代時電壓的變化量
ΔW=-JΔU
(16)
雅可比矩陣中的元素為各個功率步長對節(jié)點電壓的偏導(dǎo)數(shù),對偏導(dǎo)數(shù)化簡后得到雅可比矩陣元素關(guān)于導(dǎo)納矩陣元素間的關(guān)系,通過對雅可比矩陣中的不同節(jié)點分布規(guī)律求出該矩陣中的各個元素。
通過高斯求逆法可對修正方程進行求解,求解后可得到節(jié)點電壓的步長,從而進一步求得下一次計算的初始電壓數(shù)據(jù)。
使用Visual Studio2013和Matlab2014a作為程序開發(fā)平臺,在主頻為3.30 GHz,運行內(nèi)存為4 GB的計算機上對IEEE4系統(tǒng)進行潮流計算。在Visual Studio2013上建立C#語言的工程后,根據(jù)圖2開發(fā)潮流計算程序。牛頓法潮流計算需要先求出導(dǎo)納矩陣Y。求解導(dǎo)納矩陣時,采用關(guān)聯(lián)矩陣法[13]。根據(jù)網(wǎng)絡(luò)拓?fù)鋱D對支路編號和方向標(biāo)識,對于n個節(jié)點m條支路的有向圖,可以得到一個n行m列的關(guān)聯(lián)矩陣A。通過網(wǎng)絡(luò)已給的支路參數(shù)可得支路阻抗矩陣Z,根據(jù)以下公式可求出導(dǎo)納矩陣
Y=AZ-1AT
(17)
編程時對關(guān)聯(lián)矩陣的處理需要使用非線性映射法對矩陣中的數(shù)據(jù)進行編號,根據(jù)編號可訪問到矩陣中的元素。
圖2 牛頓法潮流計算流程圖
導(dǎo)納矩陣求出后,需要對原始數(shù)據(jù)的存放進行處理,根據(jù)節(jié)點編號的順序?qū)⒐?jié)點數(shù)據(jù)分開存放到各個對應(yīng)數(shù)組。迭代計算中編程重點在于求雅可比矩陣,雅可比矩陣中存放的是對各個電壓自變量的偏導(dǎo)數(shù),偏導(dǎo)數(shù)的計算公式也各自不一。如何根據(jù)不同情況計算偏導(dǎo)數(shù),并且存放到雅可比矩陣對應(yīng)的位置上需要使用一些相關(guān)的技術(shù)進行處理。其次如何選擇適當(dāng)?shù)姆椒ㄇ蠼庑拚匠桃埠苤匾?,常用的方法包括高斯消去法和三角分解法以及稀疏矩陣求解技巧。由于算例中的?jié)點數(shù)比較少,程序中使用三角分解法求解修正方程。對于矩陣類[14]開發(fā)的要求也就增加了,需要在該類中設(shè)計出能夠進行三角分解的方法,程序中使用的是高斯法對雅可比矩陣求逆。其中復(fù)數(shù)類的設(shè)計是便于節(jié)點中的復(fù)數(shù)數(shù)據(jù)的存放,以及潮流計算后能夠使輸出的數(shù)據(jù)更加直觀,便于與Matlab計算出來的數(shù)據(jù)比較。
根據(jù)IEEE4節(jié)點系統(tǒng)對應(yīng)的標(biāo)幺制初始節(jié)點電壓和功率數(shù)據(jù)在Visual Studio2013上使用C#語言開發(fā)的牛頓法潮流計算程序求得每次迭代節(jié)點不平衡量數(shù)據(jù)如表1所示。
表1 C#節(jié)點非平衡量
通過上面求出的節(jié)點非平衡量和修正方程,使用高斯求逆法可得到節(jié)點電壓修正量,與對應(yīng)的節(jié)點數(shù)據(jù)。根據(jù)以上方法可求得節(jié)點電壓如下表2所示。
表2 C#的節(jié)點電壓數(shù)據(jù)
根據(jù)C#牛頓潮流算法同一原理,使用Matlab2014a平臺對同一節(jié)點系統(tǒng)進行運算得到以下表3的節(jié)點非平衡量和表4的每次迭代過程中的節(jié)點電壓數(shù)據(jù)。
表3 Matlab節(jié)點非平衡量
表4 Matlab節(jié)點電壓
對以上兩個不同平臺上計算出的數(shù)據(jù)進行比較分析,C#上求解出的各個變量與Matlab上求出的各個變量相比,其相差<0.000 01。由此可得出C#用于開發(fā)潮流計算程序是比較合理的,而且使用C#的優(yōu)勢在于能夠更好的實現(xiàn)人機交互界面,不需要在Matlab與C#之間進行跨平臺操作。
文中得出在兩種平臺上算出的數(shù)據(jù)誤差很小,但運行時間相差較大。如表5所示。
表5 C#與Matlab運行時間比較
對表6中運行時間進行比較,C#的計算速度明顯比Matlab快,表明C#語言用于實現(xiàn)潮流計算算法是可行的,并且由于ASP.NET技術(shù)也是以C#為核心,這對以后開發(fā)實時潮流計算軟件提供了可靠的依據(jù)。
目前小型電力系統(tǒng)的潮流計算程序都是使用C或者C++以及Matlab等開發(fā),MatPower也只提供牛頓拉夫遜法、高斯賽德爾法、快速解耦算法(XB版)、快速解耦算法(BX版)等4種潮流計算方法[15],而且MatPower常被用來研究電力系統(tǒng)負(fù)荷[16]用于實現(xiàn)在線潮流計算比較困難。同時文中算法對于節(jié)點數(shù)較多的電力系統(tǒng)進行潮流計算時,初值選取困難,無法收斂[17]。因此,在以后的研究中需要對牛頓法數(shù)學(xué)模型進行改進;另外由于C#訪問數(shù)據(jù)庫的便利性和潮流計算的高效性等特點,在后續(xù)研究中可以充分利用C#的這些特點實現(xiàn)從算法到交互界面以及調(diào)取數(shù)據(jù)庫信息平臺的統(tǒng)一,從而快速的進行潮流計算。
[1] 張菁,陳宇晨.Matlab與VB的集成在電力系統(tǒng)潮流計算中的應(yīng)用[J]. 微計算機信息,2007,23(20):60-62.
[2] 李佳澤,陳聰,張鵬. C#與Matlab混合編程在簡單潮流計算可視化窗體設(shè)計中的應(yīng)用[J]. 數(shù)字技術(shù)與應(yīng)用,2016,33(4):159-159.
[3] 張寧,江紅梅,張渭. 基于Matlab的電力系統(tǒng)潮流計算[J]. 西北農(nóng)林科技大學(xué)學(xué)報:自然科學(xué)版,2004,32(12):124-126.
[4] 張立國,趙睿明,霍利民,等. 基于GIS的配電網(wǎng)實時潮流計算軟件的開發(fā)[J].微計算機信息,2006,22(5):287-289.
[5] 杜偉偉.石河子電網(wǎng)的潮流計算及分析[J].石河子科技,2016,40(3):53-55.
[6] 郭桂靜,徐道安,師秀鳳. 基于MatPower的電力系統(tǒng)潮流計算[J].企業(yè)導(dǎo)報,2016,16(12):172-174.
[7] 孫秋野,陳會敏,楊家農(nóng),等.牛頓類潮流計算方法的收斂性分析[J].中國電機工程學(xué)報,2014,34(13):2196-2200.
[8] 陳夢,何選森.基于八階收斂牛頓迭代的Fast-ICA改進算法[J].計算機工程與應(yīng)用,2016,52(5):2-6.
[9] 劉云舒,吳瑋華,余志遠(yuǎn),等.基于改進牛頓算法的配電網(wǎng)潮流計算[J].通信電源技術(shù),2015,32(6):162-164.
[10] 邵尉哲,王宇俊,萬新儒,等.直角坐標(biāo)與極坐標(biāo)牛頓法潮流計算速度的比較[J].南昌大學(xué)學(xué)報:工科版,2016,38(1):98-102.
[11] 段俊東,薛靜杰,栗維冰.基于牛拉法的預(yù)估校正潮流計算算法[J].河南理工大學(xué)學(xué)報:自然科學(xué)版,2015,34(3):396-399.
[12] 何仰贊,溫增銀. 電力系統(tǒng)分析[M].4版.武漢:華中科技大學(xué)出版社,2016.
[13] 陳明,李銀紅,石東源,等.節(jié)點導(dǎo)納矩陣和阻抗矩陣的互感支路組整體追加方法[J].電工技術(shù)學(xué)報,2016,31(21):94-101.
[14] 周長發(fā).C#數(shù)值計算算法編程[M].北京:電子工業(yè)出版社,2007.
[15] 陳辰,李森,李領(lǐng)娟.基于MatPower的負(fù)荷中斷分析[J].電子科技,2016,29(6):11-14.
[16] 孟祥,沈澍東.基于Matlab和MatPower的孤島判別與處理[J].電子科技,2015,28(10):123-125.
[17] Milano F.Continuous Newton’s method for power flow analysis[J].IEEE Transactions on Power Systems,2009,24(1):50-57.
Analysis of Newton Power Flow Algorithm Based on C#
WANG Yiqiang,ZHANG Liping,CHEN Yuchen,F(xiàn)AN Rui
(School of Electronic and Electrical Engineering,Shanghai University of Engineering and Science,Shanghai 201620,China)
The Matlab and C# mixed development flow calculation program of communication problems between the platform and the inconvenience,it needs to run the efficiency and reliability of Newton Ralph Xun power flow algorithm in the analysis of the C#,and verify whether the use of C# can better develop the flow calculation program,then unify development platform. The result of C# language computing time is less than 55 ms Matlab when solve the grid data of IEEE4 system,the difference in accuracy of the node data by C# compared to the data obtained by Matlab in the 10-5orders of magnitude. The experimental results show that the C# Newton power flow algorithm has high reliability and high efficiency. It can be used to develop the power flow calculation program.
admittance matrix;jacobian matrix;incidence matrix;newton method
O151.21;TM74
A
1007-7820(2017)11-063-05
2017- 01- 05
國家自然科學(xué)基金(61673257)
王醫(yī)強 1991,男,碩士研究生。研究方向:配電網(wǎng)網(wǎng)損分析等。張莉萍 1962,女,碩士,碩士生導(dǎo)師。研究方向:智能控制等。陳宇晨(1958-),男,博士,碩士生導(dǎo)師。研究方向:電力系統(tǒng)分析等。
10.16180/j.cnki.issn1007-7820.2017.11.018