本發(fā)明涉及一種電力系統(tǒng)牛頓法潮流計算方法，特別是一種適合研究目的使用的極坐標牛頓法潮流計算方法。
背景技術：
：電力系統(tǒng)潮流計算是研究電力系統(tǒng)穩(wěn)態(tài)運行的一項基本計算，它根據(jù)給定的運行條件和網(wǎng)絡結構確定整個網(wǎng)絡的運行狀態(tài)。潮流計算也是電力系統(tǒng)其他分析的基礎，如安全分析、暫態(tài)穩(wěn)定分析等都要用到潮流計算。極坐標牛頓法潮流計算方法是一種最常用的潮流計算方法，科研人員經(jīng)常以極坐標牛頓法潮流計算為基礎進行進一步地研究。實用的商業(yè)軟件采用c語言等高級編程語言編寫且使用稀疏矩陣技術和節(jié)點優(yōu)化編號等高級技術。這些技術雖然能大幅度提高潮流計算的速度、降低內(nèi)存占用量，但編程非常麻煩且難以修改和維護，不易增加新的功能，因而不適合科研人員用于研究目的使用。matlab軟件以矩陣為最基本的數(shù)據(jù)單位，可以方便地處理各種矩陣和向量運算，也可以很方便自然地處理復數(shù)類型，其指令表達式與數(shù)學中常用的形式很接近，還有大量常見實用的函數(shù)，給編程帶來很大便利。matlab軟件簡單易用、代碼短小易操作，易于編程和調(diào)試，計算功能強大，同時還具有非常強大的可視化圖形處理和交互式功能，為科學研究以及工程應用提供了一種高效的編程工具，目前已經(jīng)成為許多科學領域的基本工具和首選平臺，在各種科學和工程計算領域得到了廣泛的應用。為了適應越來越多的科研人員需要在matlab平臺上以極坐標牛頓法潮流計算為基礎進行進一步地研究的需求，迫切需要一種基于matlab軟件的易于編程、修改和調(diào)試的極坐標牛頓法潮流計算方法。根據(jù)電力系統(tǒng)節(jié)點的特點，潮流計算把電力系統(tǒng)節(jié)點分成3類：節(jié)點有功功率和無功功率已知、節(jié)點電壓幅值和電壓相角未知的節(jié)點稱為pq節(jié)點；節(jié)點有功功率和電壓幅值已知、節(jié)點無功功率和電壓相角未知的節(jié)點稱為pv節(jié)點；節(jié)點電壓幅值和電壓相角已知，節(jié)點有功功率和無功功率未知的節(jié)點稱為平衡節(jié)點。牛頓法潮流計算分為兩類：牛頓法潮流計算中節(jié)點電壓采用極坐標表示時的計算方法，稱為極坐標牛頓法潮流計算方法；牛頓法潮流計算中節(jié)點電壓采用直角坐標表示時的計算方法，稱為直角坐標牛頓法潮流計算方法。極坐標牛頓法潮流計算主要方程如下：節(jié)點導納矩陣為：式中，yik為節(jié)點導納矩陣元素，當下標i≠k時，為節(jié)點i和節(jié)點k的互導納，當下標i＝k時，為節(jié)點i的自導納；n為節(jié)點數(shù)。節(jié)點功率方程為：式中，pi、qi分別為節(jié)點i的節(jié)點有功功率和無功功率；ui、uk分別為節(jié)點i和節(jié)點k的節(jié)點電壓幅值；θik＝θi-θk，θi和θk分別為節(jié)點i和節(jié)點k的節(jié)點電壓相角；gik、bik分別為節(jié)點導納矩陣元素yik的實部和虛部。功率不平衡量方程為：式中，δpi、δqi分別為節(jié)點i的節(jié)點有功功率不平衡量和無功功率不平衡量；pis、qis分別為節(jié)點i給定的節(jié)點注入有功功率和注入無功功率；m為pq節(jié)點數(shù)。修正方程組為：式中，j為雅可比矩陣，h、n、m、l分別為雅可比矩陣的分塊子矩陣；δu/u為節(jié)點電壓幅值修正量除以電壓幅值后的列向量。雅可比矩陣各元素計算公式為：當j≠i時當j＝i時或用下列公式計算：式中，pi、qi分別為節(jié)點i的有功功率和無功功率，按式(2)計算。如圖1-2所示，現(xiàn)有極坐標牛頓法潮流計算方法，主要包括以下步驟：a、原始數(shù)據(jù)輸入和電壓初始化；原始數(shù)據(jù)包括線路和變壓器支路數(shù)據(jù)、節(jié)點注入有功功率和無功功率、節(jié)點電壓幅值、節(jié)點無功補償數(shù)據(jù)，以及收斂精度、最大迭代次數(shù)。電壓初始化采用平啟動，即pv節(jié)點和平衡節(jié)點的節(jié)點電壓幅值取給定值，pq節(jié)點的節(jié)點電壓幅值取1.0；所有節(jié)點電壓的相角都取0.0。這里相角單位為弧度，其他量采用標幺值。b、形成節(jié)點導納矩陣；根據(jù)輸入的線路和變壓器支路數(shù)據(jù)形成如式(1)所示的節(jié)點導納矩陣。c、形成雅可比矩陣；按式(5)～式(16)計算雅可比矩陣的各元素。d、計算節(jié)點功率及功率不平衡量；按式(2)計算節(jié)點功率，按式(3)計算節(jié)點功率不平衡量。e、解方程及修正節(jié)點電壓幅值u和相角θ；解修正方程組(4)，求出電壓幅值修正量列向量δu及電壓相角修正量列向量δθ。電壓修正公式為：式中，上標(t)表示第t次迭代的值；δui和δθi分別為節(jié)點i的節(jié)點電壓幅值修正量和節(jié)點電壓相角修正量。f、判斷功率最大不平衡量|δp|max和|δq|max是否都小于收斂精度ε；如果都小于收斂精度ε，進行步驟g，否則返回步驟c進行下一次迭代；g、計算平衡節(jié)點的有功功率和無功功率及pv節(jié)點的無功功率，計算各支路有功功率和無功功率，結束。直接采用上述原理實現(xiàn)的極坐標牛頓法潮流計算軟件計算速度較慢，商業(yè)使用的極坐標牛頓法潮流計算軟件采用稀疏矩陣技術和節(jié)點優(yōu)化編號技術，比較復雜，不適合科研人員以此為基礎進一步進行科學研究。因此，中國專利zl201610863886.1提出了一種基于matlab的極坐標牛頓法潮流計算方法，可以充分利用matlab特有的擅長矩陣運算和復數(shù)運算的特點，設計出了簡潔又有較快計算速度的潮流計算方法，為以極坐標牛頓法潮流計算為基礎進行進一步研究的科研人員提供一個易于修改和維護的牛頓法潮流計算方法，其特點如下：1、在matlab平臺實現(xiàn)，便于科研人員使用matlab提供的各種工具和函數(shù)對計算結果進行測試和分析；2、采用矩陣運算和復數(shù)運算，減少了程序代碼，簡化了編程，使得程序更加清晰，便于科研人員修改程序、對程序進行調(diào)試和改進、添加新功能；3、采用矩陣運算并直接調(diào)用matlab的方程求解算法，大大提高了計算速度。中國專利zl201610863886.1提出的一種基于matlab的極坐標牛頓法潮流計算方法，為從事電力系統(tǒng)研究的科研人員提供了一種基于matlab平臺的易于修改和維護且計算較為快速的極坐標牛頓法潮流計算方法。該方法采用matlab實現(xiàn)，充分利用matlab擅長矩陣運算和復數(shù)運算的特點，但未使用稀疏矩陣技術，計算速度相對較慢，仍有待進一步提高計算速度。技術實現(xiàn)要素：為解決現(xiàn)有技術存在的上述問題，本發(fā)明要提出基于matlab稀疏矩陣的極坐標牛頓法潮流計算方法，充分利用matlab特有的擅長矩陣運算和復數(shù)運算的特點，并且運用matlab提供的稀疏矩陣技術，實現(xiàn)提高潮流計算的計算速度的目的。為了實現(xiàn)上述目的，本發(fā)明的技術方案如下：基于matlab稀疏矩陣的極坐標牛頓法潮流計算方法，采用矩陣運算和復數(shù)運算并使用matlab提供的稀疏矩陣技術，包括以下步驟：a、原始數(shù)據(jù)輸入和電壓初始化；b、記錄相關節(jié)點類型的節(jié)點號；極坐標牛頓法潮流計算的修正方程組方程個數(shù)及變量個數(shù)與電力系統(tǒng)的節(jié)點類型有關，δp方程組中沒有平衡節(jié)點有功功率不平衡量對應的方程，δq方程組中僅有pq節(jié)點無功功率不平衡量對應的方程；變量則不包含平衡節(jié)點的相角變量和電壓幅值變量以及pv節(jié)點的電壓幅值變量。為了提高計算速度，形成雅可比矩陣及方程右端向量時先不考慮節(jié)點類型，形成雅可比矩陣及方程右端向量后解修正方程時，再去掉無關的行和列。為此，設置3個數(shù)組記錄有關節(jié)點類型的節(jié)點號，其中數(shù)組bt1記錄pq節(jié)點和pv節(jié)點的節(jié)點號，數(shù)組bt2記錄pq節(jié)點的節(jié)點號，數(shù)組bt記錄雅可比矩陣及方程右端向量需要的行列號。記錄相關節(jié)點類型的節(jié)點號使用matlab的find函數(shù)實現(xiàn)：bt1＝find(bus_type～＝vθ)(19)bt2＝find(bus_type＝＝pq)(20)式中，bus_type為節(jié)點類型列向量；～＝為不等于關系運算符；＝＝為等于關系運算符；vθ為平衡節(jié)點類型；pq為pq節(jié)點類型。形成數(shù)組bt1和bt2后，把數(shù)組bt2的所有元素都加上節(jié)點數(shù)n后，添加到數(shù)組bt1后形成數(shù)組bt，用來記錄雅可比矩陣及方程右端向量需要的行列號。bt＝[bt1bt2+n](21)c、形成節(jié)點導納矩陣，并轉(zhuǎn)化為稀疏矩陣y；d、形成雅可比矩陣及計算節(jié)點功率；采用matlab編程，推導出基于稀疏矩陣技術的雅可比矩陣的矩陣運算和復數(shù)運算計算方法。雅可比矩陣元素與節(jié)點類型有關，常規(guī)形成雅可比矩陣時要判斷節(jié)點類型，根據(jù)節(jié)點類型決定哪些節(jié)點需要形成雅可比矩陣元素。這樣處理對于通過循環(huán)來實現(xiàn)的算法，容易處理，但不適合通過矩陣整體運算形成雅可比矩陣的方法。因此，本發(fā)明形成雅可比矩陣時，不判斷節(jié)點類型，對所有節(jié)點都形成雅可比矩陣元素，之后再去掉不需要的行和列。采用矩陣運算計算雅可比矩陣元素和節(jié)點功率的公式如下：對式(5)～式(8)的分析，得：mij＝-nij(22)lij＝hij(23)因此，先求hij和nij，求出hij和nij后，自然就能得到mij和lij。式(5)和式(6)兩式中都包含uiuj和θij項，說明雅可比矩陣各元素中包含乘以的共軛的項，即包含項；觀察式(5)和式(6)括號內(nèi)的表達式，知雅可比矩陣各元素中還包含yij的共軛項。因此，hij、nij、mij、lij由下式生成：式中，上標(^)表示復數(shù)的共軛；為節(jié)點i的電壓相量；為節(jié)點j的電壓相量的共軛。由式(24)元素組成的矩陣為：式中，j0為雅可比初始計算矩陣。式(25)的矩陣元素為電壓相量、導納矩陣元素的共軛和電壓相量的共軛的乘積，矩陣由下式得到：式中，為節(jié)點電壓列向量形成的稀疏對角矩陣；y為稀疏導納矩陣，上標(^)表示復數(shù)的共軛；為節(jié)點電壓共軛值列向量形成的稀疏對角矩陣。式(25)中，雅可比矩陣每行元素相加就是該行對應節(jié)點的復功率，矩陣每行元素相加，通過矩陣乘以全1列向量得到：式中，為節(jié)點復功率列向量；1n×1是全部元素都為1的n階列向量。由j0得到初始雅可比矩陣分塊子矩陣為：h0＝-im(j0)(28)n0＝-re(j0)(29)m0＝re(j0)(30)l0＝-im(j0)(31)式中，h0、n0、m0、l0為初始雅可比矩陣的分塊子矩陣；re表示取矩陣元素的實部；im表示取矩陣元素的虛部。由式(28)～式(31)得到的初始雅可比矩陣分塊子矩陣的非對角元素已經(jīng)是雅可比矩陣元素了，將對角元素修正如下：由式(13)～式(16)，并考慮到sinθii＝0和cosθii＝1，得：hii＝-uiui(giisinθii-biicosθii)+qi(32)nii＝-uiui(giicosθii+biisinθii)-pi(33)mii＝uiui(giicosθii+biisinθii)-pi(34)lii＝-uiui(giisinθii-biicosθii)-qi(35)式(32)～式(35)中等式右側的第1項就是h0、n0、m0、l0的對角元素，因此只需對得到的h0、n0、m0、l0用pi和qi修正，得到雅可比矩陣分塊子矩陣對角元素。式(32)～式(35)寫成矩陣運算形式為：式中，為節(jié)點復功率形成的稀疏對角矩陣。式(28)～式(31)代入到式(36)～式(39)，得：形成雅可比矩陣及計算節(jié)點功率，包括以下步驟：d1、計算雅可比初始計算矩陣j0；d2、計算節(jié)點復功率d3、由j0和計算雅可比矩陣分塊子矩陣h、n、m和l；d4、由雅可比矩陣分塊子矩陣h、n、m和l形成雅可比矩陣；e、計算節(jié)點功率不平衡量；式(3)計算節(jié)點功率不平衡量公式寫成矩陣運算的形成為：式中，δp、δq分別為節(jié)點有功功率不平衡量列向量和無功功率不平衡量列向量；ps、qs分別為節(jié)點給定的注入有功功率列向量和注入無功功率列向量。計算最大有功功率不平衡量δpmax和最大無功功率不平衡量δqmax。f、解方程及修正電壓幅值u和相角θ；直接調(diào)用matlab軟件的解線性方程組算法解修正方程組(4)，求出電壓幅值修正量列向量δu及電壓相角修正量列向量δθ。調(diào)用matlab軟件的解線性方程組算法時，用數(shù)組bt去掉雅可比矩陣中不需要的行和列及不平衡量中不需要的行。對電壓進行修正的式(17)和式(18)改寫成矩陣形式為：u(t+1)＝u(t)-δu(t)(46)θ(t+1)＝θ(t)-δθ(t)(47)式中，上標(t)表示第t次迭代的值；δu和δθ分別為節(jié)點電壓幅值修正量列向量和節(jié)點電壓相角修正量列向量；u為節(jié)點電壓幅值列向量；θ為節(jié)點電壓相角列向量。計算節(jié)點電壓幅值和電壓相角后，按下式計算節(jié)點電壓相量g、判斷功率最大不平衡量|δp|max和|δq|max是否都小于收斂精度ε；如果都小于收斂精度ε，轉(zhuǎn)步驟h，否則返回步驟d進行下一次迭代。h、計算平衡節(jié)點的有功功率和無功功率及pv節(jié)點的無功功率，計算各支路有功功率和無功功率，結束。與現(xiàn)有技術相比，本發(fā)明具有以下有益效果：1、本發(fā)明提出的方法在matlab平臺實現(xiàn)，便于科研人員使用matlab提供的各種工具和函數(shù)對計算結果進行測試和分析。2、本發(fā)明提出的方法主要部分都采用矩陣運算和復數(shù)運算，減少了程序代碼，簡化了編程，使得程序更加清晰，便于科研人員修改程序、對程序進行調(diào)試和改進、添加新功能；使用矩陣運算也大大提高了計算速度。3、本發(fā)明采用matlab的稀疏矩陣技術，較大幅度地提高了計算速度，同時matlab的稀疏矩陣使用非常方便，可以像全矩陣一樣用行列號直接使用稀疏矩陣的元素，也不需要設計稀疏存儲結構。4、發(fā)明根據(jù)稀疏矩陣運算特點，修改了雅可比矩陣和復功率的計算公式，使得程序簡化，速度提高。實踐證明，本發(fā)明的方法既方便了科研人員對程序進行編寫、修改和調(diào)試，同時計算速度也基本接近了在c語言平臺上實現(xiàn)的速度，為科研人員的科研工作提供了一個優(yōu)秀的分析工具。附圖說明本發(fā)明共有附圖4張。其中：圖1是現(xiàn)有極坐標牛頓法潮流計算的流程圖。圖2是現(xiàn)有極坐標牛頓法形成雅可比矩陣的流程圖。圖3是本發(fā)明極坐標牛頓法潮流計算的流程圖。圖4是本發(fā)明形成雅可比矩陣和計算節(jié)點功率的流程圖。具體實施方式下面結合附圖對本發(fā)明進行進一步地說明，按照圖3-4所示流程對一個修改后的445節(jié)點實際系統(tǒng)算例進行了計算。445節(jié)點實際大型電力系統(tǒng)有445個節(jié)點，544條支路，含有大量的小阻抗支路。為了對各種方法進行比較，把這些小阻抗支路改為正常阻抗支路以滿足各種方法的要求。采用本發(fā)明和已有專利方法對445節(jié)點實際系統(tǒng)算例進行了計算，計算時收斂精度為0.00001。兩種潮流計算方法分別為：方法1：中國專利zl201610863886.1方法，未采用稀疏矩陣技術；方法2：本發(fā)明方法。兩種方法的潮流計算的計算時間見表1，計算時間不包括數(shù)據(jù)讀入和輸出及支路功率計算的時間。表1兩種極坐標牛頓法潮流計算計算時間比較潮流計算方法計算時間(s)方法10.519方法20.060從表1可見，中國專利zl201610863886.1方法沒有采用稀疏矩陣技術，計算時間較長；本發(fā)明采用稀疏矩陣技術并對某些計算公式根據(jù)稀疏矩陣的特點進行改進可以明顯提高計算速度，計算時間僅為現(xiàn)有專利方法的1/9；同時采用matlab的稀疏矩陣技術比較簡單方便。本發(fā)明可以在任何版本的matlab編程語言實現(xiàn)，但建議使用較新版本的matlab語言。本發(fā)明不局限于本實施例，任何在本發(fā)明披露的技術范圍內(nèi)的等同構思或者改變，均列為本發(fā)明的保護范圍。當前第1頁12