(完滿word版)電力系統潮流計算課程設計(完滿word版)電力系統潮流計算課程設計(完滿word版)電力系統潮流計算課程設計一、問題重述 課程設計要求1、在讀懂程序的基礎上畫出潮流計算根本流程圖2、經過輸入數據，進行潮流計算輸出結果3、對不同樣的負荷變化，解析潮流分布，寫出解析說明。4、對不同樣的負荷變化，進行潮流的調治控制，并說明調治控制的方法，并列表表示調治控制的參數變化。5、打印利用DDRTS進行潮流解析繪制的系統圖，以及潮流分布圖。 課程設計題目1、系統圖：兩個發電廠分別經過變壓器和輸電線路與四個變電所相連。變電所1變電所2變電所4變電所310kV母線35kV母線10kV母線35kV母線一次側電壓220kV一次側電壓220kV線路長為60km線路長為80km線路長為100km線路長為80km線路長為100km線路長為80km母線1母線2。。。。。。。。。。。。。。。。。。。。。。。。。。母線32*QFQ-50-22*QFS-50-22*TQN-100-22*TQN-100-2電廠一電廠二2、發電廠資料：母線1和2為發電廠高壓母線，發電廠一總裝機容量為〔400MW〕，母線3為機壓母線，機壓母線上裝機容量為〔100MW〕，最大負荷和最小負荷分別為50MW和30MW；發電廠二總裝機容量為〔200MW〕。3、變電所資料：〔1〕變電所1、2、3、4低壓母線的電壓等級分別為：10KV35KV10KV35KV〔2〕變電所的負荷分別為：50MW 40MW 50MW 60MW1〔3〕每個變電所的功率因數均為 cosφ=0.85；〔4〕變電所2和變電所4分別配有兩臺容量為75MVA的變壓器，短路耗費414KW，短路電壓〔%〕=16.7；變電所1和變電所3分別配有兩臺容量為63MVA的變壓器，短路耗費為245KW，短路電壓〔%〕=10.5；4、輸電線路資料：發電廠和變電所之間的輸電線路的電壓等級及長度標于圖中， 單位長度的電阻為 ，單位長度的電抗為 ，單位長度的電納為 2.78*10-6S。 課程設計根本內容對給定的網絡查找潮流計算所需的各元件等值參數，畫出等值電路圖。輸入各支路數據，各節點數據利用給定的程序進行在變電所在某一負荷情況下的潮流計算，并對計算結果進行解析。隨從變電所負荷按必然比率發生變化，進行潮流計算解析。4個變電所的負荷同時以2%的比率增大；4個變電所的負荷同時以2%的比率下降1和4號變電所的負荷同時以2%的比率下降，而2和3號變電所的負荷同時以2%的比率上升；在不同樣的負荷情況下，解析潮流計算的結果，若是各母線電壓不滿足要求，進行電壓的調整。〔變電所低壓母線電壓 10KV要求調整范圍在 之間；電壓35KV要求調整范圍在35-36之間〕輪流斷開環網一回線，解析潮流的分布。利用DDRTS軟件，進行繪制系統圖進行上述各種情況潮流的解析，并進行結果的比較。最后形成課程設計成品說明書。二、問題解析 節點設置及分類依照系統圖可知此系統為兩端供電網路，將母線1，2設為節點1，10，將變電所1、2、3、4的上下壓側分別設為節點2、3、4、5、6、7、8、9。并且，將節點1設為平衡節點，將節點10設為PV節點，其他節點設為PQ節點。2.2參數求取設定基準值SB100MVA，UBSB484依照題目220KV，因此ZB2UB原始資料，計算發電廠、變壓器及線路的參數。〔1〕運用以下公式計算變壓器參數：22PkUNRTUK%UN2XT1000S100SN2ZTRTXTZTZT〔2〕計算線路參數ZBZRjX(rjx)L〔3〕變電所負荷分別為：變電所1SL變電所2SL變電所3SL變電所4SL將參數整理，見下表：首端號尾端號阻抗有名值阻抗標幺值電納有名值電納標幺值12141623004500466700688900810〔4〕計算變壓器分接頭變比變壓器有5個抽頭，電壓調治范圍為UN2*2.5%，UN對應的分接頭開始時設變壓器高壓側接主接頭 ,降壓變壓器5個分接頭時的非標準變比 k*以備調壓時采用k*1220(15%)1011220k*2220(12.5%)1011220k*32201011220k*4220(12.5%)1011220k*5220(15%)1011220對變電所低壓母線為 35KV時，非標準變比與 10KV時同樣。 計算方法利用牛頓拉夫遜法進行求解，用MATLAB軟件編程，能夠求解系統潮流分布依照題目的不同樣要求對參數進行調整，經過調治變壓器變比和發電廠的電壓，3求解出合理的潮流分布，最后用 PSAT進行潮流解析，將兩者進行比較。 牛頓—拉夫遜法1、牛頓—拉夫遜法大綱第一對一般的牛頓—拉夫遜法作一簡單的說明。一個變量 X函數為：f(X) 0(0)到此方程時，由合適的近似值 X 出發，依照：(n1)f(XXX(n)f(X

(n)(n)

)(n 1,2,......))(n)屢次進行計算，當X滿足合適的收斂條件就是上面方程的根。這樣的方法就是所謂的牛頓—拉夫遜法。(n)這一方法還可以夠做下邊的講解，設第n次迭代獲取的解語真值之差，即X的誤差為時，那么：f(X(n))0把f(X(n))在X(n)周邊對用泰勒級數張開2f(X(n))f(X(n))f(X(n))f(X(n))......02!上式省略去 2今后局部f(X(n)) f(X(n)) 0(n)的誤差能夠近似由上式計算出來。(Xf(X

(n)(n)

))(n)比較兩式，能夠看出牛頓—拉夫遜法的休整量和 X 的誤差的一次項相等。用同樣的方法考慮，給出 n個變量的n個方程：f1(X1,X2, ,Xn) 0f2(X1,X2, ,Xn) 0fn(X1,X2, ,Xn) 0對其近似解X1得修正量 X1能夠經過解下邊的方程來確定：4f1(X1,X2,,Xn)f1f1f1x1x2xnX1f2(X1,X2,,Xn)f2f2f2X2x1x2xnfn(X1,X2,,Xn)fnfnfnXnx1x2xn式中等號右邊的矩陣fn都是對于X1,X2,,Xn的值。這一矩陣稱為雅可比xn〔JACOBI〕矩陣。按上述獲取的修正向量X1,X2,,Xn后，獲取以下關系XnXnXn這比X1,X2, ,Xn更湊近真實值。這一步在收斂到希望的值以前重復進行，一般要屢次計算滿足maxX1n1X1n1,X2n1X2n1,,Xnn1Xnn1為起初規定的小正數，Xnn1是第n次迭代Xn的近似值。2、牛頓法的框圖及求解過程1〕用牛頓法計算潮流時，有以下的步驟：①給這各節點電壓初始值e(0),f(0)；②將以上電壓初始值代入公式，求修正方程的常數項向量P(0), Q(0),(V2)(0)；③將電壓初始值在帶入上述公式，求出修正方程中系數矩陣的各元素。④解修正方程式e(0),f(0)；⑤修正各節點電壓e(1)e(0)e(0)，f(1)f(0)f(0)；⑥將e(1)，f(1)在帶入方程式，求出P(1),Q(1),(V2)(1)；(k) (k)⑦檢驗可否收斂，即 max Pi , Qi若是收斂，迭代到此結束，進一步計算各線路潮流和平衡節點功率，并打印輸出結果。若是不收斂，轉回②進行下次迭代計算，直到收斂為止。〔2〕程序框圖以下：5啟動輸入原始數據形成節點導納矩陣分解各節點初始電壓的實部和虛部迭代次數 K=0(k)(k)(k)(k)求PQ節點的Pi,Qi，求PV節點的Pi，Ui置節點號 i=0是雅克比矩陣可否形成，i>n否求得雅克比矩陣各元素(k)(k)Jij(k)(k)(k)Sij(k)HijNijLijRij增大節點號i=i+1把雅克比矩陣單位化求解修正方程，得 ei(k)，fi(k)回帶各電壓新值， K=K+1求解最大修正量| ei(k)|max，| fi(k)|max否可否收斂是計算輸出電壓大小及相角，節電功率及支路耗費停止三、問題求解3.1等值電路的計算電壓是衡量電力系統電能質量的標準之一。電壓過高或過低，都將對人身及其用電設備產生重要的影響。保證用戶的電壓湊近額定值是電力系統調換的根本6任務之一。當系統的電壓偏離贊同值時，電力系統必定應用電壓調治技術調治系統電壓的大小，使其保持在贊同值范圍內。本文經過手算形成了等值電路圖，并編寫好了程序得出節點電壓標幺值，使其滿足所要求的調整范圍。我們第一對給定的程序輸入局部作了簡要的解析，程序開始需要我們確定輸入節點數、支路數、平衡母線號、支路參數矩陣、節點參數矩陣。〔1〕為了保證整個系統潮流計算的完滿性，我們把凡擁有母線及發電機處均選作節點，這樣，能夠確定有10個節點，節點號見等值電路圖。〔2〕確定完節點及編號后，各條支路也相應確定了，網絡中總計有13條支路，我們對各支路參數進行了計算。依照所給實質電路圖和題中的條件，有以下公式計算各輸電線路的阻抗和對地支路電容的標幺值和變壓器的阻抗標幺值。選擇電壓基準值為UB220KV和功率基準值SB100MVA因此ZBUB484SB3.2依照題目要求及原始資料畫出系統等值電路圖 ：35791:k1:k1:k1:k12342468jjj110U=231V發電廠一δ=0P=200MW發電廠二U=231V圖1系統等值電路利用psat軟件可得計算所需系統圖如圖 2所示。7圖23.3潮流計算3.3.1B1、B2矩陣的形成：依照所求參數，以及B1矩陣的含義，列寫B1矩陣以下:B1=[1210;1410;1610;23011;45011;4610;67011;6810;89011;81010];由各個變電所負荷功率能夠計算出總功率為200MW，而發電廠一、二的總裝機容量分別為400MW和200MW。令發電廠二的功率為200MW，為了減小線路上的耗費，令發電機的電壓為額定電壓的1.05倍。依照前面表達的節點分類，形成B2矩陣以下：B2=[0023123101;00220002;0220002;800220002;0220002;00220002;0220002;00220002;0220002;200023123103]; 編寫程序并運行從結果中能夠觀察到各個節點電壓標幺值分別為 :節點12345678910電壓折算到有名值今后，能夠觀察到節點3、5、7、9的電壓都與題目給定的范圍對照偏高，因此調治變壓器分接頭和發電廠電壓，經過對系統進行屢次調整，最后獲取合理結果，調整結果以下:調治方法：調治方法電廠1電壓電廠2電壓分接頭1分接頭2分接頭3分接頭4未調整1111調整后調治結果：節點號12345678910未調整調整后由上表觀察到，進行最后一次調治節氣點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費由此能夠看出，隨著變壓器分接頭的調高有功耗費逐漸增大，吻合實質情況，而最后調治使得電壓在規定范圍內時，有功耗費增加了，相對與實質耗費小很多，能夠認為是合理的。詳盡的潮流分布見下表9各條支路的首端功率 Si 各條支路的尾端功率 Sj 各節點的功率S為(節點號從小到大排列)： 獲取的各節點的電壓標幺值， 電壓角度，節點注入無功，支路首端無功如圖所示：10值度5幺1角標壓壓0電電510-551000節點號節點號1功功無無入端注0首點路節支34567891002345678910121圖33.3.3 4個變電所負荷同時以 2%的比率增大由于變電所負荷以 2%的比率增大，因此各個變電所負荷變為：10變電所1SL變電所2SL變電所3SL變電所4SL相應地，B2矩陣變為：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];將程序中的B2矩陣代替后，運行并進行調治，此系統經過調治后到達要求，對應的調治方法和結果見下表：調治方法：調治方法電廠1電壓電廠2電壓分接頭1分接頭2分接頭3分接頭4調整前1111調整后調治結果：節點電壓12345678910調整前調整后由上表觀察到，進行最后一次調治節氣點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費11由電壓和功率耗費能夠觀察到，最后一次的電壓和有功耗費均吻合題目要求，能夠為是合理的。詳盡的潮流分布以下：各條支路的首端功率 Si 各條支路的尾端功率 Sj 各節點的功率S為(節點號從小到大排列)：獲取的圖像以下：10值度5幺1角標壓壓0電電510-551000節點號節點號1功功無無入端注0首點路節支34567891002345678910121圖43.3.4 4個變電所負荷同時以 2%的比率下降12由于變電所負荷以2%的比率增大，因此各個變電所負荷變為：變電所1SL變電所2SL變電所3SL變電所4SL相應地，B2矩陣變為：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];將程序中的B2矩陣代替后，運行并進行調治，對應的調治方法和結果見下表：調治方法：調治方法電廠1電壓電廠2電壓分接頭1分接頭2分接頭3分接頭4調整前1111調整后調治結果見下表：節點號12345678910調整前調整后由上表觀察到，進行最后一次調治節氣點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路13(8,10)支路總耗費詳盡的潮流分布以下表各條支路的首端功率Si各條支路的尾端功率Sj 各節點的功率S為(節點號從小到大排列)： 獲取的各節點的電壓標幺值，電壓角度，節點注入無功，支路首端無功情況見圖7。10值度5幺1角標壓壓0電電510-551000節點號節點號1功功無無入端注0首點路節支34567891002345678910121圖5141、4號變電所負荷同時以2%的比率下降，2、3號變電所負荷同時以2%的比率增加當1、4號變電所負荷同時以2%的比率下降，2、3號變電所負荷同時以2%的比率增加，各個變電所負荷變為：變電所1SL變電所2SL變電所3SL變電所4SL相應地，B2矩陣變為：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];將程序中的B2矩陣代替后，運行并進行調治，對應的調治方法和結果見下表：調治方法：調治方法電廠1電壓電廠2電壓分接頭1分接頭2分接頭3分接頭4調整前1111調整后調治結果見下表：節點號12345678910調整前調整后由上表觀察到，進行最后一次調治后節點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路15(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費由電壓和功率耗費能夠觀察到，最后一次的電壓和有功耗費均吻合題目要求，能夠為是合理的。詳盡的潮流分布以下表各條支路的首端功率 Si 各條支路的尾端功率 Sj 各節點的功率S為(節點號從小到大排列)： 獲取的各節點的電壓標幺值，電壓角度，節點注入無功，支路首端無功情況見圖6。1610值度5幺1角標壓壓0電電510-551000節點號節點號1功功無無入端注0首點路節支567891001234圖6 斷開1、4支路時的潮流解析當斷開1、4支路時，以致支路數變為 9和B1矩陣改變，，其他參數不變。用MATLAB軟件進行潮流計算進行調整。B1矩陣變為B1=[1210;1610;2301;4501;4610;6701;6810;8901;81010];節點號12345678910調整前調整后由上表觀察到，進行最后一次調治后節點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費 未調整 調整后(1,2)支路 (1,6)支路 (2,3)支路 (4,5)支路 17(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費由電壓和功率耗費能夠觀察到，最后一次的電壓和有功耗費均吻合題目要求，能夠為是合理的。詳盡的潮流分布以下表各條支路的首端功率 Si 各條支路的尾端功率 Sj各節點的功率S為(節點號從小到大排列)：獲取的調治后各節點的電壓標幺值，電壓角度，節點注入無功，支路首端無功情況見圖7。10值5度幺1角標壓壓0電電510-551000節點號節點號11功功無無入端注0首0點路節支34567891023456789121圖718 斷開1、6支路時的潮流解析當斷開1、6支路時，以致支路數變為 9和B1矩陣改變，，其他參數不變。用MATLAB軟件進行潮流計算進行調整。B1矩陣變為B1=[1210;1410;2301;4501;4610;6701;6810;8901;81010];節點號12345678910調整前調整后由上表觀察到，進行最后一次調治后節點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費由電壓和功率耗費能夠觀察到，最后一次的電壓和有功耗費均吻合題目要求，能夠為是合理的。詳盡的潮流分布以下表各條支路的首端功率 Si 各條支路的尾端功率 Sj 19 各節點的功率S為(節點號從小到大排列)：獲取的調治后各節點的電壓標幺值，電壓角度，節點注入有功，支路首端無功情況見圖8。15值10度幺角5標壓壓1電電0510-551000節點號節點號1功功無無入端注0首點路0節支34567891023456789121圖83.3.8斷開4、6支路時的潮流解析當斷開4、6支路時，以致支路數變為9和B1矩陣改變，，其他參數不變。用MATLAB軟件進行潮流計算進行調整。B1矩陣變為B1=[1210;1410;1610;2301;4501;6701;6810;208901;81010];表17節點號12345678910調整前調整后由上表觀察到，進行最后一次調治后節點3，5，7，9的電壓均在題目贊同的范圍內，對線路耗費進行解析，統計每次調整后各個支路的有功耗費，記錄于下表：功率耗費未調整調整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路總耗費由電壓和功率耗費能夠觀察到，最后一次的電壓和有功耗費均吻合題目要求，能夠為是合理的。詳盡的潮流分布以下表各條支路的首端功率Si各條支路的尾端功率Sj各節點的功率S為(節點號從小到大排列)：獲取的調治后各節點的電壓標幺值，電壓角度，節點注入無功，支路首端無功情況見圖9。2115值10度幺角5標壓壓1電電0510-551000節點號節點號1功功無無入端注首點0路節支345678910023456789121圖9四、誤差解析以系統在正常運行的情況下為例。利用PSAT軟件，比對MATLAB中的B1、B2矩陣，獲取各個節點的電壓情況。將PSAT軟件運行的結果對MATLAB運行的結果進行校驗。同時利用該軟件對圖形顏色的變化的描繪和數據對結果進行呈現。PSAT軟件運行的結果如圖 10所示。獲取各節點的電壓標幺值，并與 MATLAB程序運行獲取的結果進行比較。見下表節點12345678910PSATMatlab誤差0022VoltageMagnitudeProfile1.u.p[V012345678910Bus#圖10經過上述的比較可知:應用DDRTS仿真和MATLAB計算的到的結果誤差幾乎為0-0.02之間，相差很小，大體同樣，能夠認為此潮流計算是合理、正確的。23五、心得領悟及總結此次課程設計使我在潮流計算、MATLAB的使用方面均有所提高，但也裸露出了一些問題；理論知識的積蓄還是缺乏，對MATLAB的性能和特點還不能夠有一個全面的掌握，對MATLAB中PSAT工具箱使用不夠熟悉，相信經過今后的學習能填充這些缺乏，從而到達一個新的層次。潮流計算是電力系統的最根本、最常用的解析計算。用以研究系統規劃和運行中提出的各種問題。對規劃中的電力系統，經過潮流計算能夠檢驗所提出的電力系統規劃方案可否滿足各種運行方式的要求；對運行中的電力系統，經過潮流計算能夠預知各種負荷變化和網絡結構的改變會不會危及系統的安全，系統中所有母線的電壓可否在贊同的范圍以內，系統中各種元件(線路、變壓器等 )可否會出現過負荷，以及可能出現過負荷時應起初采用哪些預防措施等。實質的情況遠比我們計算的情況復雜，這讓我深刻認識了潮流計算的重要性。精準的潮流計算不但能夠使電網處于牢固且平衡的狀態，這對電網的安全運行起到要點性的作用，還可以夠實現運行的經濟性。因此我認為學好電力系統的要點在于學好潮流計算，電力系統的安全運行不是兒戲，一個很小的失誤即可能造成整個電網的崩潰，這給國家和社會帶來的危害和損失之大是可想而知的。學好電力系統這門課、學好潮流計算是我們電專業學生的當務之急。沒有牢固的基礎就不能能建筑出高樓大廈，因此，要想在今后的工作崗位上有所成就，就必定打好基礎，努力提高自己的專業涵養。參照文件1、?電力系統解析基礎?，韋剛，中國電力初版社 ;2、?電力系統穩態解析?，陳珩，中國電力初版社 ;3、?MATLAB合用指南?，蘇金明，王永利，電子工業初版社。24附錄一、正常情況下的參數二、MATLAB程序%本程序的功能是用牛頓——拉夫遜法進行潮流計算B1矩陣：1、支路首端號；2、尾端號；3、支路阻抗；4、支路對地電納5、支路的變比；6、支路首端處于K側為1，1側為0B2矩陣：1、該節點發電機功率；2、該節點負荷功率；3、節點電壓初始值4、PV節點電壓V的給定值；5、節點所接的無功補償設備的容量6、節點分類標號clear;n=10;%input(' 請輸入節點數：n=');nl=10;%input(' 請輸入支路數：nl=');isb=1;%input(' 請輸入平衡母線節點號： isb=');pr=0.00001;%input(' 請輸入誤差精度：pr=');B1=[128.5+20.1i0.000556i10;1413.6+32.16i0.0002224i10;1613.6+32.16i0.0002224i10;231.495+40.335i01.051;451.78+53.885i01.0251;4610.2+24.12i0.0001668i10;25671.495+40.335i01.0251;686.8+16.08i0.0004448i10;891.78+53.885i01.0251;8108.5+20.1i0.000556i10];%input('請輸入由支路參數形成的矩陣： B1=');B2=[0022922901;00220002;050+30.987i220002;00220002;040+27.79i220002;00220002;050+30.987i220002;00220002;060+37.18i220002;200022922903];%input('請輸入各節點參數形成的矩陣：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);%－－－－－－－更正局部－－－－－－－－－－－－ym=1;SB=100;UB=220;%ym=input('您輸入的參數是標么值？〔假設不是那么輸入一個不為零的數值〕');ifym~=0%SB=input(' 請輸入功率基準值:SB=');%UB=input(' 請輸入電壓基準值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;fori=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1 矩陣B1=');disp(B1)fori=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2 矩陣B2=');disp(B2)end%%%---------------------------------------------------26fori=1:nl %支路數ifB1(i,6)==0 %左節點處于低壓側p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非對角元Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %對角元K側Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %對角元1側end%求導納矩陣disp(' 導納矩陣Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出導納陣的實部和虛部fori=1:n %給定各節點初始電壓的實部和虛部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV節點電壓給定模值endfori=1:n %給定各節點注入功率S(i)=B2(i,1)-B2(i,2); %i節點注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i節點無功補償量end%===================================================================P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;whileIT2~=0IT2=0;a=a+1;fori=1:nifi~=isb %非平衡節點C(i)=0;D(i)=0;forj1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);% Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);% Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);% 節 點 功 率 P 計 算 ei Σ(Gij*ej-Bij*fj)+fi Σ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);% 節 點 功 率 Q 計 算 fi Σ(Gij*ej-Bij*fj)-ei Σ(Gij*fj+Bij*ej)%求P',Q'27V2=e(i)^2+f(i)^2; %電壓模平方%=========以下針對非PV節點來求取功率差及Jacobi矩陣元素=========ifB2(i,6)~=3 %非PV節點DP=P(i)-P1; %節點有功功率差DQ=Q(i)-Q1; % 節點無功功率差%===============以上為除平衡節點外其他節點的功率計算=================%=================求取Jacobi矩陣===================forj1=1:nifj1~=isb&j1~=i %非平衡節點&非對角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); %dP/df=dQ/deX3=X2; %X2=dp/dfX3=dQ/deX4=-X1; %X1=dP/deX4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseifj1==i&j1~=isb %非平衡節點&對角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;% 擴展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;% 擴展列△PJ(m,q)=X2;endendelse%===============下面是針對 PV節點來求取 Jacobi 矩陣的元素===========DP=P(i)-P1; %PV節點有功誤差DV=V(i)^2-V2; %PV節點電壓誤差forj1=1:nifj1~=isb&j1~=i %非平衡節點&非對角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); %dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;28J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseifj1==i&j1~=isb%非平衡節點&對角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendendendend%=========以上為求雅可比矩陣的各個元素 =====================fork=3:N0 %N0=2*n 〔從第三行開始，第一、二行是平衡節點〕k1=k+1;N1=N; % N=N0+1即N=2*n+1擴展列△P、△Qfork2=k1:N1 %擴展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); %非對角元規格化endJ(k,k)=1; % 對角元規格化ifk~=3 % 不是第三行%============================================================k4=k-1;fork3=3:k4 %用k3行從第三行開始到當前行前的 k4行消去fork2=k1:N1 %k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去運算endJ(k3,k)=0;endifk==N0break;end%==========================================fork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去運算endJ(k3,k)=0;29endelsefork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去運算endJ(k3,k)=0;endendend%==== 上面是用線性變換方式將 Jacobi矩陣化成單位矩陣=====fork=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); % 更正節點電壓實部k1=k+1;f(L)=f(L)-J(k1,N); % 更正節點電壓虛部end%------ 更正節點電壓-----------fork=3:N0DET=abs(J(k,N));ifDET>=pr % 電壓誤差量可否滿足要求IT2=IT2+1;% 不滿足要求的節點數加 1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp(' 迭代次數：');disp(ICT1);disp(' 沒有到達精度要求的個數： ');disp(ICT2);fork=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%===============計算各輸出量===========================disp(' 各節點的實質電壓標幺值 E為(節點號從小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp(' 各節點的電壓大小 V為(節點號從小到大排列)：');30disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp(' 各節點的電壓相角 sida為(節點號從小到大排列)：');disp(sida);forp=1:nC(p)=0;forq=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp(