1、目錄第一章課程設計概述11.1 設計目的11.2 設計要求11.3 設計題目11.4 設計內容21.5設計時間2第二章（針對該課程設計題目的）設計思路32.1 潮流計算題目32.2 對課題的分析及求解思路3第三章電力系統潮流計算概述53.1電力系統敘述53.2潮流計算簡介53.3潮流計算的意義及其發展6第四章導納矩陣的原理及計算方法84.1自導納和互導納的確定方法84.2潮流計算的基本方程124.3電力系統節點分類154.4潮流計算的約束條件16第五章牛頓拉夫遜法及程序代碼185.1牛頓拉夫遜法18牛頓拉夫遜法概要18牛頓拉夫遜法的求解過程19程序框圖如下：215.2 潮流計算程序代碼21第六

2、章運行結果分析306.1 正常運行結果分析306.1.1 MATLAB運行結果306.1.2 PowerWorld運行結果356. 1. 3 綜合分析366.2 非正常運行結果分析376.2.1 若某發電站發電量減半時對系統供電影響376.2.2 若某發電站發電量為零時對系統供電影響396.2.3 若某負荷切掉時對系統供電影響406.2.4 若系統環網解開對系統的影響42第七章總結44參考文獻45第一章 課程設計概述1.1 設計目的1. 掌握電力系統潮流計算的基本原理和電力系統運行方式的變化2. 讓學生結合課程教學內容了解實際電力系統分析中的某些最基本的計算、分析過程、常規典型計算方法和結果分

3、析等方面知識，加深電力系統分析課程的理解和掌握3. 掌握并能熟練運用一門計算機語言（MATLAB語言或C語言或C+語言）；4. 采用計算機語言對潮流計算進行計算機編程計算。5. 培養和提高學生對電力系統分析與設計能力。 1.2 設計要求1. 程序源代碼；2. 選定算例的輸入，輸出文件；3. 程序說明；4. 選定算例的程序計算過程；5. 選定算例的手算過程（至少迭代2次）（可選）。學生開展課程設計可以7人選擇同一題目，但每個人必須獨立完成所有給定任務，學生必須按計劃開展課程設計, 每一階段完成情況必須接受教師檢查。編程實現部分，學生可以根據自己掌握的編程工具情況選擇，如MATLAB、C、VC等，

4、 最后必須提供程序代碼，課程設計必須完成對系統的分析與設計任務,編寫相應的分析與設計報告。 1.3 設計題目1.高壓輸電網潮流的計算機算法程序設計（PQ分解法、牛頓-拉夫遜法）2.中壓配電網潮流的計算機算法程序設計（前推后代法、同倫延拓法等）3.電力系統短路故障的計算機算法程序設計（要求不限）本設計選擇為高壓輸電網潮流的計算機算法程序設計（PQ分解法、牛頓-拉夫遜法）1.4 設計內容1.根據電力系統網絡推導電力網絡數學模型，寫出節點導納矩陣；2.賦予各節點電壓變量（直角坐標系形式）初值后，求解不平衡量；3.形成雅可比矩陣；4.求解修正量后，重新修改初值，從2開始重新循環計算；5.求解的電壓變量

5、達到所要求的精度時，再計算各支路功率分布、功率損耗和平衡節點功率; 6.上機編程調試；連調；7.計算分析給定系統潮流分析。8.準備計算機演示答辯，書寫該課程設計說明書（必須計算機打印）。1.5設計時間2013年春季17周至18周第二章（針對該課程設計題目的）設計思路2.1 潮流計算題目G 1 2 3 -10Mvar 180MW 100MvarG 5 4 40MW 30Mvar50MW 30Mvar如圖所示，變壓器參數為：Zt1=3+j110, Zt2=0.8+j23,變比分別為231/110kV,231/121kV;220kV線路參數為：Z12=5.9+j31.5;110kV線路：Z45=65

6、+j100; 34段，Z=65+j100，所有阻抗都已經歸算至220kV側。2.2 對課題的分析及求解思路等值電路的計算電壓是衡量電力系統電能質量的標準之一。電壓過高或過低，都將對人身及其用電設備產生重大的影響。保證用戶的電壓接近額定值是電力系統調度的基本任務之一。當系統的電壓偏離允許值時，電力系統必須應用電壓調節技術調節系統電壓的大小，使其維持在允許值范圍內。本文經過手算形成了等值電路圖，并編寫好了程序得出節點電壓標幺值，使其滿足所要求的調整范圍。我們首先對給定的程序輸入部分作了簡要的分析，程序開始需要我們確定輸入節點數、支路數、平衡母線號、支路參數矩陣、節點參數矩陣。（1）為了保證整個系統

7、潮流計算的完整性，我們把凡具有母線及發電機處均選作節點，這樣，可以確定有5個節點，節點號見等值電路圖。（2）確定完節點及編號后，各條支路也相應確定了，網絡中總計有5條支路，我們對各支路參數進行了計算。根據所給實際電路圖和題中的已知條件，有以下公式計算各輸電線路的阻抗和對地支路電容的標幺值和變壓器的阻抗標幺值。選擇電壓基準值為和功率基準值所以數據的輸入該系統中，母線1為平衡節點，保持定值電壓Ù=220+j0kV,5節點為PV節點，注入有功功率為40MW，無功功率為 30Mvar，其余節點均為PQ節點，注入功率分別為2=-j10MVA，3=-180-j100MVA,4=-50-j30MV

8、A.首先計算節點數據矩陣如下，%(bus#)(volt) (ang) (p) (q) (bus type)bus= 1 1.00 0.00 0.00 0.00 3; 2 1.00 0.00 0.00 -0.10 1; 3 1.00 0.00 -1.80 -1.00 1; 4 1.00 0.00 -0.50 -0.30 1; 5 1.00 0.00 0.40 0.30 2;然后計算線路數據矩陣如下， % b#1 b#2 (R) (X) (G) (B) (K)line= 1 2 5.9 31.5 0.0 0.0 0; 1 5 3 110 0.0 0.0 1.0; 2 3 0.8 23 0.0 0.

9、0 1.0; 3 4 65 100 0.0 0.0 0; 4 5 65 100 0.0 0.0 0;這就形成了計算潮流的原始數據，把它代入程序中就可以計算潮流分布。第三章 電力系統潮流計算概述3.1電力系統敘述 電力工業發展初期，電能是直接在用戶附近的發電站（或稱發電廠）中生產的，各發電站孤立運行。隨著工農業生產和城市的發展，電能的需要量迅速增加，而熱能資源和水能資源豐富的地區又往往遠離用電比較集中的城市和工礦區，為了解決這個矛盾，就需要在動力資源豐富的地區建立大型發電站，然后將電能遠距離輸送給電力用戶。同時，為了提高供電的可靠性以及資源利用的綜合經濟性，又把許多分散的各種形式的發電站，通過送

10、電線路和變電所聯系起來。這種由發電機、升壓和降壓變電所，送電線路以及用電設備有機連接起來的整體，即稱為電力系統。現代電力系統提出了“靈活交流輸電和新型直流輸電”的概念。靈活交流輸電技術是指運用固態電子器件與現代自動控制技術對交流電網的電壓、相位角、阻抗、功率以及電路的通斷進行實時閉環控制，從而提高高壓輸電線路的訴訟能力和電力系統的穩態水平。新型直流輸電技術是指應用現電力電子技術的最新成果，改善和簡化變流站的造價等。運營方式管理中，潮流是確定電網運行方式的基本出發點：在規劃領域，需要進行潮流分析驗證規劃方案的合理性；在實時運行環境，調度員潮流提供了電網在預想操作預想下的電網的潮流分布以及校驗運行

11、的可靠性。在電力系統調度運行的多個領域都涉及到電網潮流計算。潮流是確定電力網咯運行狀態的基本因素，潮流問題是研究電力系統穩態問題的基礎和前提。3.2潮流計算簡介電力系統潮流計算是研究電力系統穩態運行情況的一種計算，它根據給定的運行條件及系統接線情況確定整個電力系統各部分的運行狀態：各母線的電壓。各元件中流過的功率，系統的功率損耗等等。在電力系統規劃的設計和現有電力系統運行方式的研究中，都需要利用潮流計算來定量的分析比較供電方案或運行方式的合理性。可靠性和經濟性。此外，電力系統的潮流計算也是計算機系統動態穩定和靜態穩定的基礎，所以潮流計算是研究電力系統的一種和重要和基礎的計算。電力系統潮流計算也

12、分為離線計算和在線計算兩種，前者主要用于系統規劃設計和安排系統的運行方式，后者則用于正在運行系統的經常監視及實時控制。利用電子數字計算機進行潮流計算從50年代中期就已經開始了。在這20年內，潮流計算曾采用了各種不同的方法，這些方法的發展主要圍繞著對潮流計算的一些基本要求進行的，對潮流計算的要求可以歸納為以下幾點：（1） 計算方法的可靠性或收斂性；（2） 對計算機內存量的要求；（3） 計算速度；（4） 計算的方便性和靈活性。3.3潮流計算的意義及其發展電力系統潮流計算是電力系統分析中的一種最基本的計算，是對復雜電力系統正常和故障條件下穩態運行狀態的計算。潮流計算的目標是求取電力系統在給定運行狀態

13、的計算，即節點電壓和功率分布，用以檢查系統各元件是否過負荷。各點電壓是否滿足要求，功率的分布和分配是否合理以及功率損耗等。對現有的電力系統的運行和擴建，對新的電力系統進行規劃設計以及對電力系統進行靜態和穩態分析都是以潮流計算為基礎。潮流計算結果可用如電力系統穩態研究，安全估計或最優潮流等對潮流計算的模型和方法有直接影響。實際電力系統的潮流技術那主要采用牛頓拉夫遜法。運行方式管理中，潮流是確定電網運行方式的基本出發點；在規劃領域，需要進行潮流分析驗證規劃方案的合理性；在實時運行環境，調度員潮流提供了多個在預想操作情況下電網的潮流分布以及校驗運行可靠性。在電力系統調度運行的多個領域問題是研究電力系

14、統穩態問題的基礎和前提。在用數字見算計算機解電力系統潮流問題的開始階段，普遍采取以節點導納矩陣為基礎的逐次代入法。這個方法的原理比較簡單，要求的數字計算機內存量比較差下，適應50年代電子計算機制造水平和當時電力系統理論水平，但它的收斂性較差，當系統規模變大時，迭代次數急劇上升，在計算中往往出現迭代不收斂的情況。這就迫使電力系統的計算人員轉向以阻抗矩陣為基礎的逐次代入法。阻抗法改善了系統潮流計算問題的收斂性，解決了導納無法求解的一些系統的潮流計算，在60年代獲得了廣泛的應用，阻抗法德主要缺點是占用計算機內存大，每次迭代的計算量大。當系統不斷擴大時，這些缺點就更加突出，為了克服這些缺點，60年代中

15、期發展了以阻抗矩陣為基礎的分塊阻抗法。這個方法把一個大系統分割為幾個小的地區系統，在計算機內只需要存儲各個地區系統的阻抗矩陣及它們之間聯絡的阻抗，這樣不僅大幅度的節省了內存容量，同時也提高了計算速度?？朔杩狗ㄈ秉c是另一個途徑是采用牛頓-拉夫遜法。這是數學中解決非線性方程式的典型方法，有較好的收斂性。在解決電力系統潮流計算問題時，是以導納矩陣為基礎的，因此，只要我們能在迭代過程中盡可能保持方程式系數矩陣的稀疏性，就可以大大提高牛頓法潮流程序的效率。自從60年代中期，牛頓法中利用了最佳順序消去法以后，牛頓法在收斂性。內存要求。速度方面都超過了阻抗法，成為了60年代末期以后廣泛采用的優秀方法。 本

16、次設計是利用牛頓拉夫遜法進行求解，用MATLAB軟件編程，可以求解系統潮流分布根據題目的不同要求對參數進行調整，通過調節變壓器變比和發電廠的電壓，求解出合理的潮流分布，最后用powerworld進行潮流分析，將兩者進行比較。第四章 導納矩陣的原理及計算方法4.1自導納和互導納的確定方法 電力網絡的節點電壓方程： 為節點注入電流列向量，注入電流有正有負，注入網絡的電流為正，流出網絡的電流為負。根據這一規定，電源節點的注入電流為正，負荷節點為負。既無電源又無負荷的聯絡節點為零，帶有地方負荷的電源節點為二者代數之和。為節點電壓列向量，由于節點電壓是對稱于參考節點而言的，因而需先選定參考節點。在電力系

17、統中一般以地為參考節點。如整個網絡無接地支路，則需要選定某一節點為參考。設網絡中節點數為（不含參考節點），則，均為n*n列向量。為n*n階節點導納矩陣。節電導納矩陣的節點電壓方程： 展開為： ：是一個n*n階節點導納矩陣，其階數就等于網絡中除參考節點外的節點數。 節點導納矩陣的對角元素 (i=1，2，n)成為自導納。自導納數值上就等于在i節點施加單位電壓，其他節點全部接地時，經節點i注入網絡的電流，因此，它可以定義為：節點i的自導納數值上就等于與節點直接連接的所有支路導納的總和。節點導納矩陣的非對角元素 (j=1，2，n;i=1，2，。，n;j=i)稱互導納，由此可得互導納數值上就等于在節點i

18、施加單位電壓，其他節點全部接地時，經節點j注入網絡的電流，因此可定義為：節點j，i之間的互導納數值上就等于連接節點j，i支路到導納的負值。顯然，恒等于。互導納的這些性質決定了節點導納矩陣是一個對稱稀疏矩陣。而且，由于每個節點所連接的支路數總有一個限度，隨著網絡中節點數的增加非零元素相對愈來愈少，節點導納矩陣的稀疏度，即零元素數與總元素的比值就愈來愈高。節點導納矩陣的性質及意義節點導納矩陣的性質：（1）為對稱矩陣，=。如網絡中含有源元件，如移相變壓器，則對稱性不再成立。（2）對無接地支路的節點，其所在行列的元素之和均為零，即 。對于有接地支路的節點，其所在行列的元素之和等于該點接地支路的導納。利

19、用這一性質，可以檢驗所形成節點導納矩陣的正確性。（3）具有強對角性：對角元素的值不小于同一行或同一列中任一元素（4）為稀疏矩陣，因節點i ，j 之間無支路直接相連時=0，這種情況在實際電力系統中非常普遍。矩陣的稀疏性用稀疏度表示，其定義為矩陣中的零元素與全部元素之比，即 ， 式中Z 為中的零元素。S 隨節點數n 的增加而增加：n=50，S可達92%；n=100，S 可達90%；n=500，S可達99%，充分利用節點導納矩陣的稀疏性可節省計算機內存，加快計算速度，這種技巧稱為稀疏技術。節點導納矩陣的意義：是n*n階方陣，其對角元素 (i=1，2，-n)稱為自導納，非對角元素(i，j=1，2，n，

20、)稱為互導納。將節點電壓方程展開為可見， 表明，自導納在數值上等于僅在節點i施加單位電壓而其余節點電壓均為零（即其余節點全部接地）時，經節點i注入網絡的電流。其顯然等于與節點i直接相連的所有支路的導納之和。同時可見。表明，互導納在數值上等于僅在節點j施加單位電壓而其余節點電壓均為零時，經節點i注入網絡的電流，其顯然等于()即=。為支路的導納，負號表示該電流流出網絡。如節點ij之間無支路直接相連，則該電流為0，從而=0。注意字母幾種不寫法的不同意義：粗體黑字表示導納矩陣，大寫字母代矩陣中的第i行第j列元素，即節點i和節點j之間的互導納。小寫字母i，j支路的導納等于支路阻抗的倒數數，。根據定義直接

21、求取節點導納矩陣時，注意以下幾點：1)。節點導納矩陣是方陣，其階數就等于網絡中除去參考節點外的節點數。參考節點一般取大地，編號為零。2)。節點導納矩陣是稀疏矩陣，其各行非零非對角元素就等于與該行相對應節點所連接的不接地支路數。3)。節點導納矩陣的對角元素就等于各該節點所連接導納的總和。因此，與沒有接地支路的節點對應的行或列中，對角元素為非對角元素之和的負值。4)。節點導納矩陣的非對角元素等于連接節點i，j支路導納的負值。因此，一般情況下，節點導納矩陣的對角元素往往大于非對角元素的負值。5)。節點導納矩陣一般是對稱矩陣，這是網絡的互易特性所決定的。從而，一般只要求求取這個矩陣的上三角或下三角部分

22、。非標準變比變壓器等值電路變壓器型等值電路更便于計算機反復計算,更適宜于復雜網絡的潮流計算.雙繞組變壓器可用阻抗與一個理想變壓器串聯的電路表示.理想變壓器只是一個參數,那就是變比?，F在變壓器阻抗按實際變比歸算到低壓側為例,推導出變壓器型等值電路.a 雙繞組變壓器原理圖b 變壓器阻抗歸算到低壓側等值模型流入和流出理想變壓器的功率相等式中,是理想變壓器的變比,和 分別為變壓器高,低繞組的實際電壓.從圖b直接可得：從而可得: 式中,又因節點電流方程應具有如下形式:將式（1-8）與（1-9）比較，得：因此可得各支路導納為: 由此可得用導納表示的變壓器型等值電路:圖 c4.2潮流計算的基本方程在潮流問題

23、中，任何復雜的電力系統都可以歸納為以下元件（參數）組成。（1）發電機（注入電流或功率）（2）負荷（注入負的電流或功率）（3）輸電線支路（電阻，電抗）（4）變壓器支路（電阻，電抗，變比）（5）母線上的對地支路（阻抗和導納）（6）線路上的對地支路（一般為線路充電點容導納）集中了以上各類型的元件的簡單網絡如圖 (a) 潮流計算用的電網結構圖(b) 潮流計算等值網絡采用導納矩陣時，節點注入電流和節點電壓構成以下線性方程組其中 可展開如下形式 由于實際電網中測量的節點注入量一般不是電流而是功率，因此必須將式中的注入電流用節點注入功率來表示。節點功率與節點電流之間的關系為 式中，因此用導納矩陣時，PQ節點

24、可以表示為把這個關系代入式中 ，得就是電力系統潮流計算的數學模型-潮流方程。它具有如下特點：（1）它是一組代數方程，因而表征的是電力系統的穩定運行特性。（2）它是一組非線性方程，因而只能用迭代方法求其數值解。（3）由于方程中的電壓和導納既可以表為直角坐標，又可表為極坐標，因而潮流方程有多種表達形式-極坐標形式，直角坐標形式和混合坐標形式。a。取 ，得到潮流方程的極坐標形式：b。 取 ， ，得到潮流方程的直角坐標形式：c。取，得到潮流方程的極坐標形式：不同坐標形式的潮流方程適用于不同的迭代解法。例如：利用牛頓-拉夫遜迭代法求解，以直角坐標和極坐標形式的潮流方程為方便；而P-Q解耦法是在極坐標形式

25、的基礎上發展而成，本次設計采用牛頓-拉夫遜迭代法求解且采用極坐標形式。（4）它是一組n個復數方程，因而實數方程數為2n個但方程中共含4n個變量：P，Q，U和，i=1，2，n，故必須先指定2n個變量才能求解。4.3電力系統節點分類用一般的電路理論求解網絡方程，目的是給出電壓源(或電流源)研究網絡內的電流(或電壓)分布，作為基礎的方程式，一般用線性代數方程式表示。然而在電力系統中，給出發電機或負荷連接母線上電壓或電流(都是向量)的情況是很少的，一般是給出發電機母線上發電機的有功功率(P)和母線電壓的幅值(U)，給出負荷母線上負荷消耗的有功功率(P)和無功功率(Q)。主要目的是由這些已知量去求電力系

26、統內的各種電氣量。所以，根據電力系統中各節點性質的不同，很自然地把節點分成三類： PQ節點對這一類點，事先給定的是節點功率(P，Q)，待求的未知量是節點電壓向量(U，)，所以叫PQ節點。通常變電所母線都是PQ節點，當某些發電機的輸出功率P。Q給定時，也作為PQ節點。PQ節點上的發電機稱之為PQ機(或PQ給定型發電機)。在潮流計算中，系統大部分節點屬于PQ節點。 PU節點這類節點給出的參數是該節點的有功功率P及電壓幅值U，待求量為該節點的無功功率Q及電壓向量的相角。這類節點在運行中往往要有一定可調節的無功電源。用以維持給定的電壓值。通常選擇有一定無功功率儲備的發電機母線或者變電所有無功補償設備的

27、母線做PU節點處理。PU節點上的發電機稱為PU機(或PU給定型發電機) 平衡節點在潮流計算中，這類節點一般只設一個。對該節點，給定其電壓值，并在計算中取該節點電壓向量的方向作為參考軸，相當于給定該點電壓向量的角度為零。也就是說，對平衡節點給定的運行參數是U和，因此有城為U節點，而待求量是該節點的P。Q，整個系統的功率平衡由這一節點承擔。關于平衡節點的選擇，一般選擇系統中擔任調頻調壓的某一發電廠(或發電機)，有時也可能按其他原則選擇，例如，為提高計算的收斂性?？梢赃x擇出線數多或者靠近電網中心的發電廠母線作平衡節點。以上三類節點4個運行參數P。Q。U。中，已知量都是兩個，待求量也是兩個，只是類型不

28、同而已。4.4潮流計算的約束條件電力系統運行必須滿足一定技術和經濟上的要求。這些要求夠成了潮流問題中某些變量的約束條件，常用的約束條件如下：1. 節點電壓應滿足2.從保證電能質量和供電安全的要求來看，電力系統的所有電氣設備都必須運行在額定電壓附近。PU節點電壓幅值必須按上述條件給定。因此，這一約束條件對PQ節點而言。3. 節點的有功功率和無功功率應滿足4.PQ節點的有功功率和無功功率，以及PU節點的有功功率，在給定是就必須滿足上述條件，因此，對平衡節點的P和Q以及PU節點的Q應按上述條件進行檢驗。5. 節點之間電壓的相位差應滿足 為了保證系統運行的穩定性，要求某些輸電線路兩端的電壓相位不超過一

29、定的數值。這一約束的主要意義就在于此。 因此，潮流計算可以歸結為求解一組非線性方程組，并使其解答滿足一定的約束條件。常用的方法是迭代法和牛頓法，在計算過程中，或得出結果之后用約束條件進行檢驗。如果不能滿足要求，則應修改某些變量的給定值，甚至修改系統的運行方式，重新進行計算。第五章牛頓拉夫遜法及程序代碼5.1牛頓拉夫遜法牛頓拉夫遜法概要首先對一般的牛頓拉夫遜法作一簡單的說明。已知一個變量X函數為：到此方程時，由適當的近似值出發，根據：反復進行計算，當滿足適當的收斂條件就是上面方程的根。這樣的方法就是所謂的牛頓拉夫遜法。這一方法還可以做下面的解釋，設第次迭代得到的解語真值之差，即的誤差為時，則：把

30、在附近對用泰勒級數展開上式省略去以后部分的誤差可以近似由上式計算出來。比較兩式，可以看出牛頓拉夫遜法的休整量和的誤差的一次項相等。用同樣的方法考慮，給出個變量的個方程：對其近似解得修正量可以通過解下邊的方程來確定：式中等號右邊的矩陣都是對于的值。這一矩陣稱為雅可比（JACOBI）矩陣。按上述得到的修正向量后，得到如下關系這比更接近真實值。這一步在收斂到希望的值以前重復進行，一般要反復計算滿足為預先規定的小正數，是第n次迭代的近似值。牛頓拉夫遜法的求解過程（1）用牛頓法計算潮流時，有以下的步驟：給這各節點電壓初始值；將以上電壓初始值代入公式，求修正方程的常數項向量 ；將電壓初始值在帶入上述公式，

31、求出修正方程中系數矩陣的各元素。解修正方程式；修正各節點電壓，；將，在帶入方程式，求出；檢驗是否收斂，即如果收斂，迭代到此結束，進一步計算各線路潮流和平衡節點功率，并打印輸出結果。如果不收斂，轉回進行下次迭代計算，直到收斂為止。程序框圖如下：5.2 潮流計算程序代碼主程序：clc;clear allformat long;global nSW nPQ nPV; global nb; global nl; global bus; global line; global Y; global nodenum; global lPQ; global myf; openfile; change; ybu

32、s; NR; PQ; flow; ret; output; 打開bus line的數據文件的子程序：function openfileglobal nb nl;global line;global bus;global myf;dfile,pathname=uigetfile('*.m','Select Data File'); if pathname = 0 error(' you must select a valid data file')else lfile =length(dfile);% strip off .m eval(dfile

33、(1:lfile-2); endnl,ml=size(line); nb,mb=size(bus); outfile=''for I=1:length(dfile)if dfile(I)='.'break;end outfile=strcat(outfile,dfile(I); endoutfile=strcat(outfile,'.txt'); myf=fopen(outfile,'w'); 節點重新編號：function changeglobal nb;global nl;global nPQ;global bus;globa

34、l line;global nodenum;nSW = 0; nPV = 0; nPQ = 0; for I = 1:nb, type= bus(I,6);if type = 3, nSW = nSW + 1; SW(nSW,:)=bus(I,:);elseif type = 2, nPV = nPV +1; PV(nPV,:)=bus(I,:);else nPQ = nPQ + 1; PQ(nPQ,:)=bus(I,:);endendbus=PQ;PV;SW; newbus=1:nb'nodenum=newbus bus(:,1); bus(:,1)=newbus; for I=1:

35、nl for J=1:2for k=1:nbif line(I,J)=nodenum(k,2) line(I,J)=nodenum(k,1);breakendendendend建立節點導納矩陣：functionY = ybus(bus,line)global nl;global nb;global bus;global line;global Y;global myf;Y=zeros(nb,nb); for k=1:nlI=line(k,1); J=line(k,2); Zt=line(k,3)+j*line(k,4);if J=0 Yt=1/Zt; end%Ym=line(k,5)+j*li

36、ne(k,6); K=line(k,7); if (K=0)&(J=0) Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J);endif (K=0)&(J=0) Y(I,I)=Y(I,I)+Ym;endif K>0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/K/K; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J);endif K<0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+K*K*Yt;

37、Y(I,J)=Y(I,J)+K*Yt; Y(J,I)=Y(I,J);endend牛拉法解方程組：function NRglobal nb;global nPQ;global bus;global line;global Y;global myf;max1=100; eps1=1.0e-10; eps2=1.0e-10;fprintf(myf, '-節點導納矩陣Y-n'); for I=1:nbfor J=1:nb fprintf(myf, '%10f+j*(%10f) ', real(Y(I,J),imag(Y(I,J);end fprintf(myf, 

38、9;n');endfor i=1:max1 angl(:,1)=bus(1:nb-1,3); u(:,1)=bus(1:nPQ,2); x=angl;u; Jac=form_jac(bus,Y); del=dPQ(Y,bus); dx=Jacdel; fprintf(myf,'-第%d次迭代結果-n',i); fprintf(myf,'-第%d次迭代的雅比矩陣J-n',i);for I=1:nb+nPQ-1for J=1:nb+nPQ-1 fprintf(myf, '%10f ', Jac(I,J); end fprintf(myf, &

39、#39;n');end fprintf(myf,'-第%d次迭代的功率偏差dP和dQ-n',i);for I=1:nb+nPQ-1 fprintf(myf, '%10e ', del(I,1); fprintf(myf, 'n');end fprintf(myf,'-第%d次迭代的節點相角和電壓的偏差dx-n',i);for I=1:nb+nPQ-1 fprintf(myf, '%10e ', dx(I,1); fprintf(myf, 'n');endfor I=nb:nb+nPQ-1 d

40、x(I,1)=dx(I,1)*x(I,1); end x=x-dx; fprintf(myf,'-第%d次迭代的節點相角delta和電壓U-n',i); angl=x(1:nb-1,1); u=x(nb:nb+nPQ-1,1);for I=1:nb-1 bus(I,3)=angl(I,1); fprintf(myf, 'ang%d %10f n', I,angl(I,1);endfor I=1:nPQ % bus(I,2)=u(I,1); fprintf(myf, 'U%d %10f n', I,u(I,1);endif(max(abs(dx)&

41、lt;eps1)&(max(abs(dPQ(Y,bus)<eps2) breakendendif i=max1 error('超過最大迭代次數，不收斂停機！');end計算節點注入功率：global nb;global nPQ;global bus;global line;global Y;for I=nPQ+1:nb if bus(I,6)=3 sum=0; for J=1:nb ang=bus(I,3)-bus(J,3); A=real(Y(I,J)*cos(ang)+imag(Y(I,J)*sin(ang); sum=sum+bus(I,2)*bus(J,2

42、)*(A);end bus(I,4)=sum; end sum=0;for J=1:nb ang=bus(I,3)-bus(J,3); B=real(Y(I,J)*sin(ang)-imag(Y(I,J)*cos(ang); sum=sum+bus(I,2)*bus(J,2)*B;endbus(I,5)=sum; end計算線路功率及損耗：global nl;global lPQ;global bus;global line;for k=1:nlI=line(k,1); J=line(k,2); lPQ(k,1)=I; lPQ(k,2)=J;if J=0 Zt=line(k,3)+j*line

43、(k,4); Yt=1/Zt;end Ym=line(k,5)+j*line(k,6); K=line(k,7); Ui=bus(I,2)*(cos(bus(I,3)+j*sin(bus(I,3); if J=0 Uj=bus(J,2)*(cos(bus(J,3)+j*sin(bus(J,3); endif (K=0)&(J=0) Iij=Ui*(Yt+Ym)-Uj*Yt; Iji=Uj*(Yt+Ym)-Ui*Yt;endif (K=0)&(J=0) Iij=Ui*Ym; Iji=0;endif K>0 Iij=(Ui-Uj)*Yt/K+Ui*(Ym+Yt*(K-1)/K

44、); Iji=(Uj-Ui)*Yt/K+Uj*Yt*(1-K)/K2;endif K<0 K=-K; Iij=(Ui-Uj)*Yt*K+Ui*(Ym+Yt*(1-K); Iji=(Uj-Ui)*Yt*K+Uj*Yt*K*(K-1);end Sij=Ui*conj(Iij); Sji=Uj*conj(Iji);delS=Sij+Sji; lPQ(k,3:5)=Sij Sji delS; end恢復原節點編號：global nl;global nb;global bus;global line;global Y;global nodenum;global lPQ;for I=1:nb for

45、 J=1:nb if nodenum(J,2)=Ibreak;endendtem(I,:)=bus(nodenum(J,1),:); tem(I,1)=I;endbus=tem;for I=1:nl for J=1:2for k=1:nbif lPQ(I,J)=nodenum(k,1) lPQ(I,J)=nodenum(k,2);breakendendendend輸出：function output global nb;global nl;global bus;global lPQ;global myf;fprintf(myf, 'n-牛頓拉夫遜法潮流計算結果-n');fprintf(myf, ' 節點計算結果：n'); fprintf(myf, '節點 節點電壓 節點相角（角度） 節點注入功率n');for I=1:nb, fprintf(myf, ' %2d %10f %10f %10f+j%10f n', bus(I,1),bus(