本發(fā)明涉及一種電力系統(tǒng)快速分解法潮流計算方法，特別是一種適合研究目的使用的快速分解法潮流計算方法。
背景技術：
：電力系統(tǒng)潮流計算是研究電力系統(tǒng)穩(wěn)態(tài)運行的一項基本計算，它根據(jù)給定的運行條件和網(wǎng)絡結構確定整個網(wǎng)絡的運行狀態(tài)。潮流計算也是電力系統(tǒng)其他分析的基礎，如安全分析、暫態(tài)穩(wěn)定分析等都要用到潮流計算。由于具有收斂可靠、計算速度快及內存需求少的優(yōu)點，快速分解法成為當前潮流計算的主流方法之一，科研人員經(jīng)常以快速分解法潮流計算為基礎進行進一步地研究。實用的商業(yè)軟件采用c語言等高級編程語言編寫且使用稀疏矩陣技術和節(jié)點優(yōu)化編號等高級技術。這些技術雖然能大幅度提高潮流計算的速度、降低內存占用量，但編程非常麻煩且難以修改和維護，不易增加新的功能，因而不適合科研人員用于研究目的使用。matlab軟件以矩陣為最基本的數(shù)據(jù)單位，可以方便地處理各種矩陣和向量運算，也可以很方便自然地處理復數(shù)類型，其指令表達式與數(shù)學中常用的形式很接近，還有大量常見且實用的函數(shù)，給編程帶來很大便利。matlab軟件簡單易用、代碼短小易操作，易于編程和調試，計算功能強大，同時還具有非常強大的可視化圖形處理和交互式功能，為科學研究以及工程應用提供了一種高效的編程工具，目前已經(jīng)成為許多科學領域的基本工具和首選平臺，在各種科學和工程計算領域得到了廣泛的應用。為了適應越來越多的科研人員需要在matlab平臺上以快速分解法潮流計算為基礎進行進一步地研究的需求，迫切需要一種基于matlab軟件的易于編程、修改和調試的快速分解法潮流計算方法。如圖1-2所示，現(xiàn)有快速分解法潮流計算方法，主要包括以下步驟：a、輸入原始數(shù)據(jù)和初始化電壓；根據(jù)電力系統(tǒng)節(jié)點的特點，潮流計算把電力系統(tǒng)節(jié)點分成3類：節(jié)點有功功率和無功功率已知、節(jié)點電壓幅值和電壓相角未知的節(jié)點稱為pq節(jié)點；節(jié)點有功功率和電壓幅值已知、節(jié)點無功功率和電壓相角未知的節(jié)點稱為pv節(jié)點；節(jié)點電壓幅值和電壓相角已知，節(jié)點有功功率和無功功率未知的節(jié)點稱為平衡節(jié)點。電壓初始化采用平啟動，即pv節(jié)點和平衡節(jié)點的電壓幅值取給定值，pq節(jié)點的電壓幅值取1.0；所有電壓的相角都取0.0。這里相角單位為弧度，其他量單位采用標幺值。b、形成節(jié)點導納矩陣；c、形成修正方程的系數(shù)矩陣b′和b″并進行因子表分解；潮流計算的基本方程是非線性方程組，通常采用逐次線性化方法迭代求解。線性化得到的方程稱為修正方程，用來求電壓幅值和相角的修正量?？焖俜纸夥ㄐ拚匠淌窃跇O坐標牛頓法潮流計算修正方程基礎上解耦并改進得到的?？焖俜纸夥ㄐ拚匠虨椋篵′δθ＝δp/v(1)b″δv＝δq/v(2)式中，δp/v和δq/v分別為有功功率和無功功率不平衡量除以電壓幅值后的列向量；δv和δθ分別為電壓幅值和電壓相角修正量列向量；b′為導納矩陣的虛部，但計算時不計及支路電阻、對地導納和非標準變比，導納矩陣中包含pq節(jié)點和pv節(jié)點相關的行和列；b″為導納矩陣的虛部，僅包括與pq節(jié)點有關的行和列。d、設置迭代計數(shù)t＝0，設置收斂標志kp＝0，kq＝0；e、計算有功功率不平衡量δp；pq節(jié)點和pv節(jié)點的有功功率不平衡量為：式中，pis為節(jié)點i的注入有功功率；vi為節(jié)點i的電壓幅值；θij＝θi-θj，θi、θj分別為節(jié)點i和節(jié)點j的電壓相角；gij和bij分別為導納矩陣元素的電導部分和電納部分；n為節(jié)點數(shù)。求各節(jié)點中有功功率不平衡量絕對值最大的值，稱為有功功率最大不平衡量，記為δpmax。f、判斷有功功率最大不平衡量絕對值|δpmax|是否小于收斂精度ε；如果小于收斂精度ε，令kp＝1，轉到步驟g；否則，解修正方程b'δθ＝δp/v，修正電壓相角，令kp＝0，轉到步驟h；求解修正方程b′δθ＝δp/v，得到δθ，按下式修正電壓相角：θ(t+1)＝θ(t)-δθ(t)(4)式中，上標(t)表示第t次迭代；θ為節(jié)點電壓相角列向量。g、判斷kq是否等于1；如果kq＝1，轉到步驟l；h、計算無功功率不平衡量δq；pq節(jié)點的無功功率不平衡量為：式中，qis為節(jié)點i的注入無功功率；m為pq節(jié)點數(shù)。求各節(jié)點中無功功率不平衡量絕對值最大的值，稱為無功功率最大不平衡量，記為δqmax。i、判斷無功功率最大不平衡量絕對值|δqmax|是否小于收斂精度ε；如果小于收斂精度ε，令kq＝1，轉到步驟j；否則，解修正方程b"δv＝δq/v，修正電壓幅值，令kq＝0，轉到步驟k；求解修正方程b″δv＝δq/v，得到δv，按下式修正電壓幅值列向量v：v(t+1)＝v(t)-δv(t)(6)j、判斷kp是否等于1；如果kp＝1，轉到步驟l；k、令t＝t+1，返回步驟e進行下一次迭代；l、計算平衡節(jié)點功率及pv節(jié)點的無功功率，計算支路功率，結束。步驟e和步驟f為p～θ迭代，即通過δp求δθ進而修正θ；步驟h和步驟i為q～v迭代，即通過δq求δv進而修正v。主流快速分解法都是按上述步驟設計方法，即先進行p～θ迭代，后進行q～v迭代。也有文獻采用先進行q～v迭代，后進行p～θ迭代的方法。直接采用上述原理實現(xiàn)的快速分解法潮流計算軟件計算速度較慢，商業(yè)使用的快速分解法潮流計算軟件采用稀疏矩陣技術和節(jié)點優(yōu)化編號技術，比較復雜，不適合科研人員以此為基礎進一步進行科學研究。因此，中國專利zl201710056328.9提出了一種基于matlab的快速分解法潮流計算方法，可以充分利用matlab特有的擅長矩陣運算和復數(shù)運算的特點，設計出了簡潔又有較快計算速度的潮流計算方法，其特點如下：(1)在matlab平臺實現(xiàn)，便于科研人員使用matlab提供的各種工具和函數(shù)對計算結果進行測試和分析；(2)采用矩陣運算和復數(shù)運算，減少了程序代碼，簡化了編程，使得程序更加清晰，便于科研人員修改程序、對程序進行調試和改進、添加新功能；(3)采用矩陣運算并直接調用matlab的三角分解法方程求解算法，大大提高了計算速度。中國專利zl201710056328.9所提出方法以快速分解法潮流計算為基礎，為進一步研究的科研人員提供了一個易于修改和維護的快速分解法潮流計算方法。該方法采用matlab實現(xiàn)，充分利用matlab擅長矩陣運算和復數(shù)運算的特點，但未使用稀疏矩陣技術，計算速度相對較慢，仍有待進一步提高計算速度。技術實現(xiàn)要素：為解決現(xiàn)有技術存在的上述問題，本發(fā)明要提出一種基于matlab稀疏矩陣技術的快速分解法潮流計算方法，在充分利用matlab特有的擅長矩陣運算和復數(shù)運算特點的基礎上，采用matlab的稀疏矩陣技術，設計出具有較快計算速度的潮流計算方法。為了實現(xiàn)上述目的，本發(fā)明的技術方案如下：基于matlab稀疏矩陣的快速分解法潮流計算方法，采用matlab的稀疏矩陣技術，包括以下步驟：a、輸入原始數(shù)據(jù)和初始化電壓；電壓初始化采用平啟動，形成節(jié)點電壓相量列向量同時形成節(jié)點電壓幅值列向量初值v(0)和節(jié)點電壓相角單位相量列向量初值b、記錄相關節(jié)點類型的節(jié)點號；快速分解法修正方程組的方程個數(shù)及變量個數(shù)與電力系統(tǒng)的節(jié)點類型有關，p～θ迭代方程組中沒有平衡節(jié)點有功功率不平衡量對應的方程和平衡節(jié)點相角變量；q～v迭代方程組中僅有pq節(jié)點無功功率不平衡量對應的方程和pq節(jié)點電壓幅值變量。為了提高計算速度，形成方程組系數(shù)矩陣及方程右端向量時先不考慮節(jié)點類型，形成系數(shù)矩陣及方程右端向量后，再去掉無關的行和列。為此，設置兩個數(shù)組記錄有關節(jié)點類型的節(jié)點號，其中數(shù)組bt1記錄pq節(jié)點和pv節(jié)點的節(jié)點號，數(shù)組bt2記錄pq節(jié)點的節(jié)點號。記錄相關節(jié)點類型的節(jié)點號使用matlab的find函數(shù)實現(xiàn)：bt1＝find(bus_type～＝vθ)(7)bt2＝find(bus_type＝＝pq)(8)式中，bus_type為節(jié)點類型列向量；～＝為不等于關系運算符；＝＝為等于關系運算符；vθ為平衡節(jié)點類型；pq為pq節(jié)點類型。c、形成節(jié)點導納矩陣，并轉化為稀疏的導納矩陣y，按bt2提取導納矩陣y的相應各行，形成僅包含pq節(jié)點對應行的稀疏導納矩陣子矩陣ypq；形成節(jié)點導納矩陣的步驟如下：c1、預定義導納矩陣y的維數(shù)為n×n；c2、根據(jù)線路參數(shù)和變壓器支路參數(shù)形成導納矩陣y的元素；c3、根據(jù)無功補償設備參數(shù)修正導納矩陣y部分對角元素；c4、把導納矩陣y轉化為稀疏矩陣；d、形成修正方程的稀疏系數(shù)矩陣b′和b″并進行因子表分解；為了提高計算速度和簡化程序，形成系數(shù)矩陣b′時不考慮節(jié)點類型，形成n階方陣，然后再按數(shù)組bt1記錄的節(jié)點號提取矩陣元素，去掉平衡節(jié)點對應的行和列，形成新的系數(shù)矩陣b′；按數(shù)組bt2記錄的節(jié)點號提取矩陣ypq中pq節(jié)點對應的列，取矩陣元素的虛部，形成系數(shù)矩陣b″。直接調用matlab軟件的lu分解法對系數(shù)矩陣b′進行三角分解形成下三角矩陣l1和上三角矩陣u1；對系數(shù)矩陣b″進行三角分解形成下三角矩陣l2和上三角矩陣u2。分解后得到的矩陣l1、u1、l2和u2都不包含無關的行和列，在迭代過程解方程時不用再提取矩陣元素。e、形成節(jié)點注入有功功率和無功功率向量；潮流計算迭代過程中，計算節(jié)點有功功率不平衡量向量和節(jié)點無功功率不平衡量向量時，要用到節(jié)點注入有功功率列向量ps和節(jié)點注入無功功率列向量qs，為了提高計算速度，先形成節(jié)點注入有功功率向量和節(jié)點注入無功功率向量。節(jié)點注入有功功率列向量為：ps＝pg-pl(9)式中，ps為節(jié)點注入有功功率列向量；pg為節(jié)點發(fā)電有功功率列向量；pl為節(jié)點負荷有功功率列向量。節(jié)點注入無功功率列向量為：qs＝qg-ql(10)式中，qs為節(jié)點注入無功功率列向量；qg為節(jié)點發(fā)電無功功率列向量；ql為節(jié)點負荷無功功率列向量。形成向量ps和qs時不考慮節(jié)點類型，然后再按數(shù)組bt1和bt2記錄的節(jié)點號提取向量元素，去掉多余的元素。按數(shù)組bt1記錄的節(jié)點號提取向量ps需要的元素，去掉平衡節(jié)點對應的元素，形成新的向量ps；按數(shù)組bt2記錄的節(jié)點號提取向量qs需要的元素，去掉平衡節(jié)點和pv節(jié)點對應的元素，形成新的向量qs。f、設置迭代計數(shù)t＝0，設置收斂標志kp＝0，kq＝0；g、計算有功功率不平衡量δp，并求有功功率最大不平衡量δpmax；采用matlab矩陣運算和復數(shù)運算編程，需要推導出基于矩陣運算和復數(shù)運算的功率計算方法。定義節(jié)點i的復功率公式為：式中，為節(jié)點i的復功率；pi和qi分別為節(jié)點i的有功功率和無功功率；為節(jié)點電壓相量；為節(jié)點電流相量的共軛，上標(^)表示復數(shù)的共軛。式(11)寫成向量相乘的形式為：式中，為節(jié)點復功率列向量；為節(jié)點電壓相量列向量；為節(jié)點電流相量的共軛值列向量；.*表示兩向量對應的元素相乘。式(12)中，節(jié)點電流相量為：將式(13)代入式(12)，得：式中，為節(jié)點電壓相量列向量；為節(jié)點電壓相量的共軛值列向量；y為稀疏導納矩陣；上標(^)表示復數(shù)的共軛。節(jié)點有功功率p為：式中，p為節(jié)點有功功率列向量；re表示取矩陣元素的實部。計算節(jié)點有功功率不平衡量的矩陣運算的形成為：式中，δp為節(jié)點有功功率不平衡量列向量；ps為節(jié)點注入有功功率列向量。求各節(jié)點中有功功率不平衡量絕對值最大的值，稱為有功功率最大不平衡量，記為δpmax。h、判斷有功功率最大不平衡量絕對值|δpmax|是否小于收斂精度ε；如果小于收斂精度ε，令kp＝1，轉到步驟i；否則，解式(17)所示的修正方程，然后按式(18)修正電壓相角，計算電壓相量列向量b'δθ＝δp/v(17)θ(t+1)＝θ(t)-δθ(t)(18)式中，上標(t)表示第t次迭代的值；δθ為節(jié)點電壓相角修正量列向量。使用步驟d形成的下三角矩陣l1和上三角矩陣u1直接調用matlab軟件的解線性方程組算法解修正方程組(17)。計算電壓相角后，按下式計算電壓相量由于三角函數(shù)計算工作量較大，在q～v迭代時，角度沒有更新，沒有必要重新計算角度的正弦和余弦，因此把式(19)進行簡化，用下式計算電壓相量式中，為節(jié)點i的電壓相角單位相量，在q～v迭代時不變，無須重新計算，為式(21)和式(20)寫成矩陣運算形式分別為：式中，為節(jié)點電壓相量列向量，v為節(jié)點電壓幅值列向量，為節(jié)點電壓相角單位相量列向量。令kp＝0，轉到步驟j。i、判斷kq是否等于1；如果kq＝1，轉到步驟n；j、計算無功功率不平衡量δq，并求無功功率最大不平衡量δqmax；按式(14)計算節(jié)點復功率但q～v迭代時平衡節(jié)點和pv節(jié)點的無功功率不需要計算，因此把式(14)修改為：式中，為pq節(jié)點復功率列向量；為pq節(jié)點電壓相量列向量；為節(jié)點電壓相量的共軛值列向量；ypq為僅包含pq節(jié)點對應行的稀疏導納矩陣子矩陣。按下式計算pq節(jié)點無功功率qpq：式中，qpq為pq節(jié)點無功功率列向量；im表示取矩陣元素的虛部。計算pq節(jié)點無功功率不平衡量矩陣運算的形成為：式中，δqpq為pq節(jié)點無功功率不平衡量列向量；qs為pq節(jié)點注入無功功率列向量。求各節(jié)點中無功功率不平衡量絕對值最大的值，稱為無功功率最大不平衡量，記為δqmax。k、判斷無功功率最大不平衡量絕對值|δqmax|是否小于收斂精度ε；如果小于收斂精度ε，令kq＝1，轉到步驟l；否則，解式(27)所示的修正方程，然后按式(28)修正電壓幅值，計算電壓相量列向量b"δv＝δq/v(27)v(t+1)＝v(t)-δv(t)(28)式中，上標(t)表示第t次迭代的值；δv為節(jié)點電壓幅值修正量列向量。使用步驟d形成的三角矩陣l2和u2直接調用matlab軟件的解線性方程組算法解修正方程組(27)。計算電壓幅值后，按式(23)計算電壓相量列向量令kq＝0，轉到步驟m。l、判斷kp是否等于1；如果kp＝1，轉到步驟n；m、令t＝t+1，返回步驟g進行下一次迭代；n、計算平衡節(jié)點的有功功率和無功功率及pv節(jié)點的無功功率，計算各支路有功功率和無功功率，結束。與現(xiàn)有技術相比，本發(fā)明具有以下有益效果：1、本發(fā)明提出的方法在matlab平臺實現(xiàn)，便于科研人員使用matlab提供的各種工具和函數(shù)對計算結果進行測試和分析。2、本發(fā)明提出的方法采用矩陣運算和復數(shù)運算，減少了程序代碼，簡化了編程，使得程序更加清晰，便于科研人員修改程序、對程序進行調試和改進、添加新功能；使用矩陣運算也大大提高了計算速度。3、本發(fā)明采用matlab的稀疏矩陣技術，較大幅度地提高了計算速度，同時matlab的稀疏矩陣使用非常方便，可以像全矩陣一樣用行列號直接使用稀疏矩陣的元素，也不需要設計稀疏存儲結構。4、本發(fā)明修改了導納矩陣的計算過程，進一步提高了計算速度。實踐證明，本發(fā)明的方法既方便了科研人員對程序進行編寫、修改和調試，同時計算速度也基本接近了在c語言平臺上實現(xiàn)的速度，為科研人員的科研工作提供了一個優(yōu)秀的分析工具。附圖說明本發(fā)明共有附圖4張。其中：圖1是現(xiàn)有快速分解法潮流計算的流程圖。圖2是現(xiàn)有快速分解法潮流計算形成導納矩陣的流程圖。圖3是本發(fā)明快速分解法潮流計算的流程圖。圖4是本發(fā)明快速分解法潮流計算形成導納矩陣的流程圖。具體實施方式下面結合附圖及實施例對本發(fā)明進行進一步地說明。按照圖3-4所示流程對一個修改后的445節(jié)點實際系統(tǒng)算例進行了計算。445節(jié)點實際大型電力系統(tǒng)有445個節(jié)點，544條支路，含有大量的小阻抗支路。為了對各種方法進行比較，把這些小阻抗支路改為正常阻抗支路以滿足各種方法的要求。采用本發(fā)明和已有專利方法對445節(jié)點實際系統(tǒng)算例進行了計算，計算時收斂精度為0.00001。兩種潮流計算方法分別為：方法1：中國專利zl201710056328.9方法，未采用稀疏矩陣技術；方法2：本發(fā)明方法。兩種方法的潮流計算的計算時間見表1，計算時間不包括數(shù)據(jù)讀入和輸出及支路功率計算的時間。表1兩種快速分解法潮流計算計算時間比較潮流計算方法計算時間(s)方法10.1487方法20.0575從表1可見，中國專利zl201710056328.9方法沒有采用稀疏矩陣技術，計算時間較長；本發(fā)明采用稀疏矩陣技術并對某些計算過程根據(jù)稀疏矩陣的特點進行改進可以明顯提高計算速度，同時采用matlab的稀疏矩陣技術比較簡單方便。為了進一步優(yōu)化算法，本發(fā)明在計算導納矩陣和系數(shù)矩陣b'時，采用以下不同處理方法：方法1：導納矩陣和系數(shù)矩陣b'都初始化為n×n階稀疏矩陣，所有元素初值為零；方法2：導納矩陣和系數(shù)矩陣b'都初始化為n×n階稀疏矩陣，所有元素初值為零，預留n+2nbr個存儲空間，其中nbr為支路數(shù)；方法3：導納矩陣和系數(shù)矩陣b'都初始化為n×n階全矩陣，所有元素初值為零，先計算補償設備導納，再計算支路導納，形成導納矩陣和系數(shù)矩陣b'后再轉化為稀疏矩陣；方法4：導納矩陣和系數(shù)矩陣b'都初始化為n×n階全矩陣，所有元素初值為零，先計算支路導納，再計算補償設備導納，形成導納矩陣和系數(shù)矩陣b'后再轉化為稀疏矩陣。幾種不同處理導納矩陣和系數(shù)矩陣b'方法的潮流計算的計算時間見表2，計算時間不包括數(shù)據(jù)讀入和輸出及支路功率計算的時間。表2幾種快速分解法潮流計算計算時間比較潮流計算方法計算時間(s)方法10.0714方法20.0693方法30.0616方法40.0575從表2可見，直接把矩陣初始化為稀疏矩陣雖然比采用全矩陣速度快，但并不是最好的方法，先按全矩陣形成矩陣元素，再轉化為稀疏矩陣是比較好的處理方法。在計算導納矩陣時，先計算支路導納，再用無功補償設備參數(shù)修正導納矩陣對角元比較快。本發(fā)明可以在任何版本的matlab編程語言實現(xiàn)，但建議使用較新版本的matlab語言。本發(fā)明不局限于本實施例，任何在本發(fā)明披露的技術范圍內的等同構思或者改變，均列為本發(fā)明的保護范圍。當前第1頁12