電力系統(tǒng)潮流計算22450

1、電力系統(tǒng)課程設(shè)計題 目: 電力系統(tǒng)潮流計算 院系名稱： 電氣工程學(xué)院 專業(yè)班級： 電氣F1206班 學(xué)生姓名： 學(xué) 號： 指導(dǎo)教師： 張孝遠(yuǎn) 成績：指導(dǎo)老師簽名： 日期： 目錄原始資料21 概述42 潮流計算節(jié)點介紹42.1 變量的分類52.2 節(jié)點的分類53 計算方法簡介63.1 牛頓拉夫遜法原理63.1.1 牛頓拉夫遜法概要63.1.2 牛頓法的框圖及求解過程83.2 MATLAB簡介94 潮流分布計算104.1 系統(tǒng)的一次接線圖104.2 參數(shù)計算104.3 豐大及枯大下地潮流分布情況144.3.1 該地區(qū)變壓器的有功潮流分布數(shù)據(jù)154.3.2 重、過載負(fù)荷元件統(tǒng)計表175 設(shè)計心得17

2、參考文獻18附錄:程序19原始資料一、系統(tǒng)接線圖見附件1。二、系統(tǒng)中包含發(fā)電廠、變電站、及其間的聯(lián)絡(luò)線路。500kV變電站以外的系統(tǒng)以一個等值發(fā)電機代替。各元件的參數(shù)見附件2。設(shè)計任務(wù)1、手動畫出該系統(tǒng)的電氣一次接線圖，建立實際網(wǎng)絡(luò)和模擬網(wǎng)絡(luò)之間的聯(lián)系。2、根據(jù)已有資料，先手算出各元件的參數(shù)，后再用Matlab表格核算出各元件的參數(shù)。3、潮流計算1）對兩種不同運行方式進行潮流計算，注意110kV電網(wǎng)開環(huán)運行。2）注意將電壓調(diào)整到合理的范圍110kV母線電壓控制在106kV117kV之間;220kV母線電壓控制在220 kV242kV之間。附件一：4附件二：1、變壓器：兩個220kV變電站均采用

3、參數(shù)一致的三繞組變壓器，具體參數(shù)如下。220kV變電站參數(shù)表高壓側(cè)繞組中壓側(cè)繞組低壓側(cè)繞組容量12012060（120）電壓220110(121)10.5110kV及以下的變電站的變壓器省略，即可將負(fù)荷直接掛在110kV母線上。而110kV升壓變只計及以下參數(shù)。110kV變電站參數(shù)表序號變電站名容量X1X01A162B283C1205D636E207G31.58F12.5+31.59水電站130MW、2*20MVA0.82250 0.52500 10水電站272MW、3*31.5MVA0.34444 0.33333 11水電站324MW、31.5MVA1.03333 0.33333 12水電站

4、54*7.5MW、2*20MVA1.03333 0.33333 13水電站418MW、30MVA1.28333 1.05000 14火電廠50MW、31.5+40MVA0.48285 0.26250 2、線路：具體參數(shù)如下。220kV線路參數(shù)表序號線路名稱導(dǎo)線牌號線路長度km1ML2x240102ML2x240103MH2x3005110kV線路參數(shù)表序號線路名稱導(dǎo)線牌號線路長度km1水電站1水電站2150302AB95803BC9514水電站4C150655BD95636CD240607水電站2C240758ED24014.79水電站3D15010.510水電站5T節(jié)點705.511T節(jié)點D

5、70512T節(jié)點G240513H F1503.914H D185715水電站2L240/213016LG2401017DF150418火電廠D24013、發(fā)電機各發(fā)電機的參數(shù)如下：XdXq裝機容量功率因數(shù)水電站10.1740.174300.85水電站20.1740.174720.85水電站30.1740.174240.85水電站50.1740.174300.85水電站40.1740.174180.85華鑫電廠0.1740.17450.85出力情況：水力發(fā)電機豐大出力70%,枯大出力20%?；鹆Πl(fā)電機豐大出力80%,枯大出力80%。4、負(fù)荷各110kV變電站豐大負(fù)荷按該站變電容量的50%估算，枯

6、大負(fù)荷按該站變電容量的60%估算。兩個220kV變電站的低壓側(cè)上各掛10MW的負(fù)荷，中壓側(cè)各掛20MW負(fù)荷。功率因素均為0.95。5、并聯(lián)電容器兩個220kV變電站的低壓側(cè)上均裝設(shè)并聯(lián)補償。補償總量按該站變電容量的20%裝設(shè)，分組原則以每組電容器的容量不超過10MVar且經(jīng)濟性較好為準(zhǔn)。1 概述潮流計算是電力系統(tǒng)最基本最常用的計算。根據(jù)系統(tǒng)給定的運行條件，網(wǎng)絡(luò)接線及元件參數(shù)，通過潮流計算可以確定各母線的電壓，包括電壓的幅值和相角，各元件流過的功率，整個系統(tǒng)的功率損耗等一系列系統(tǒng)中的潮流數(shù)據(jù)。近幾年，對潮流算法的研究仍然是如何改善傳統(tǒng)的潮流算法，即高斯-塞德爾法、牛頓法和快速解耦法。牛頓法，由于

7、其在求解非線性潮流方程時采用的是逐次線性化的方法，為了進一步提高算法的收斂性和計算速度，人們考慮采用將泰勒級數(shù)的高階項或非線性項也考慮進來，于是產(chǎn)生了二階潮流算法。后來又提出了根據(jù)直角坐標(biāo)形式的潮流方程是一個二次代數(shù)方程的特點，提出了采用直角坐標(biāo)的保留非線性快速潮流算法。潮流計算在數(shù)學(xué)上是多元非線性方程組的求解問題，求解的方法有很多種，牛頓拉夫遜Newton-Raphson法是數(shù)學(xué)上解非線性方程組的有效方法，有較好的收斂性。將N-R法用于潮流計算是以導(dǎo)納矩陣為基礎(chǔ)的，由于利用了導(dǎo)納矩陣的對稱性，稀疏性及節(jié)點編號順序優(yōu)劃等技巧，使N-R法在收斂性，占用內(nèi)存，計算速度等方面的優(yōu)點都超過了阻抗法總結(jié)

8、為在電力系統(tǒng)運行方式和規(guī)劃方案的研究中都需要進行潮流計算以比較運行方式或規(guī)劃供電方案的可行性、可靠性和經(jīng)濟性。同時為了實時監(jiān)控電力系統(tǒng)的運行狀態(tài)也需要進行大量而快速的潮流計算。因此潮流計算是電力系統(tǒng)中應(yīng)用最廣泛、最基本和最重要的一種電氣運算。在系統(tǒng)規(guī)劃設(shè)計和安排系統(tǒng)的運行方式時采用離線潮流計算在電力系統(tǒng)運行狀態(tài)的實時監(jiān)控中則采用在線潮流計算。2 潮流計算節(jié)點介紹常規(guī)的電力系統(tǒng)潮流計算中一般具有三種類型的節(jié)點：PQ、PV及平衡節(jié)點。一個節(jié)點有四個變量，即注入有功功率、注入無功功率，電壓大小及相角。常規(guī)的潮流計算一般給定其中的二個變量：PQ節(jié)點（注入有功功率及無功功率），PV節(jié)點（注入有功功率及電

9、壓的大?。胶夤?jié)點（電壓的大小及相角）。2.1 變量的分類負(fù)荷消耗的有功、無功功率、電源發(fā)出的有功、無功功率、母線或節(jié)點的電壓大小和相位、 、在這十二個變量中，負(fù)荷消耗的有功和無功功率無法控制，因它們?nèi)Q于用戶，它們就稱為不可控變量或是擾動變量。電源發(fā)出的有功無功功率是可以控制的自變量，因此它們就稱為控制變量。母線或節(jié)點電壓的大小和相位角是受控制變量控制的因變量。其中, 、主要受、的控制, 、主要受、的控制。這四個變量就是簡單系統(tǒng)的狀態(tài)變量。為了保證系統(tǒng)的正常運行必須滿足以下的約束條件：對控制變量 對沒有電源的節(jié)點則為對狀態(tài)變量的約束條件則是 對某些狀態(tài)變量還有如下的約束條件 2.2 節(jié)點的

10、分類 第一類稱PQ節(jié)點。等值負(fù)荷功率、和等值電源功率、是給定的，從而注入功率、是給定的，待求的則是節(jié)點電壓的大小和相位角。屬于這類節(jié)點的有按給定有功、無功率發(fā)電的發(fā)電廠母線和沒有其他電源的變電所母線。 第二類稱PV節(jié)點。等值負(fù)荷和等值電源的有功功率、是給定的，從而注入有功功率是給定的。等值負(fù)荷的無功功率和節(jié)點電壓的大小 也是給定的。待求的則是等值電源的無功功率，從而注入無功功率和節(jié)點電壓的相位角。有一定無功功率儲備的發(fā)電廠和有一定無功功率電源的變電所母線都可以作為PV節(jié)點； 第三類平衡節(jié)點。潮流計算時一般只設(shè)一個平衡節(jié)點。等值負(fù)荷功率、是給定的，節(jié)點電壓的大小和相位也是給定的。擔(dān)負(fù)調(diào)整系統(tǒng)頻率

11、任務(wù)的發(fā)電廠母線往往被選作為平衡節(jié)點。3 計算方法簡介3.1 牛頓拉夫遜法原理3.1.1 牛頓拉夫遜法概要首先對一般的牛頓拉夫遜法作一簡單的說明。已知一個變量X函數(shù)為： 到此方程時，由適當(dāng)?shù)慕浦党霭l(fā)，根據(jù)：反復(fù)進行計算，當(dāng)滿足適當(dāng)?shù)氖諗織l件就是上面方程的根。這樣的方法就是所謂的牛頓拉夫遜法。這一方法還可以做下面的解釋，設(shè)第次迭代得到的解語真值之差，即的誤差為時，則：把在附近對用泰勒級數(shù)展開上式省略去以后部分的誤差可以近似由上式計算出來。比較兩式，可以看出牛頓拉夫遜法的休整量和的誤差的一次項相等。用同樣的方法考慮，給出個變量的個方程：對其近似解得修正量可以通過解下邊的方程來確定：式中等號右邊的

12、矩陣都是對于的值。這一矩陣稱為雅可比（JACOBI）矩陣。按上述得到的修正向量后，得到如下關(guān)系這比更接近真實值。這一步在收斂到希望的值以前重復(fù)進行，一般要反復(fù)計算滿足為預(yù)先規(guī)定的小正數(shù)，是第n次迭代的近似值。3.1.2 牛頓法的框圖及求解過程1、用牛頓法計算潮流時，有以下的步驟：（1） 給這各節(jié)點電壓初始值；（2） 將以上電壓初始值代入公式，求修正方程的常數(shù)項向量 ；（3） 將電壓初始值在帶入上述公式，求出修正方程中系數(shù)矩陣的各元素。（4） 解修正方程式；（5） 修正各節(jié)點電壓，；（6） 將，在帶入方程式，求出；（7） 檢驗是否收斂，即（8） 如果收斂，迭代到此結(jié)束，進一步計算各線路潮流和平衡

13、節(jié)點功率，并打印輸出結(jié)果。如果不收斂，轉(zhuǎn)回（2）進行下次迭代計算，直到收斂為止。2、程序框圖如下圖3.1.2牛頓法框圖3.2 MATLAB簡介 MATLAB是用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算的高級技術(shù)計算語言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分。是由美國mathworks公司發(fā)布的主要面對科學(xué)計算、可視化以及交互式程序設(shè)計的高科技計算環(huán)境。它將數(shù)值分析、矩陣計算、科學(xué)數(shù)據(jù)可視化以及非線性動態(tài)系統(tǒng)的建模和仿真等諸多強大功能集成在一個易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計以及必須進行有效數(shù)值計算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳

14、統(tǒng)非交互式程序設(shè)計語言（如C、Fortran）的編輯模式，代表了當(dāng)今國際科學(xué)計算軟件的先進水平。MATLAB是一種交互式、面向?qū)ο蟮某绦蛟O(shè)計語言廣泛應(yīng)用于工業(yè)界與學(xué)術(shù)界主要用于矩陣運算同時在數(shù)值分析、自動控制模擬、數(shù)字信號處理、動態(tài)分析、繪圖等方面也具有強大的功能。MATLAB程序設(shè)計語言結(jié)構(gòu)完整且具有優(yōu)良的移植性它的基本數(shù)據(jù)元素是不需要定義的數(shù)組。它可以高效率地解決工業(yè)計算問題特別是關(guān)于矩陣和矢量的計算。MATLAB與C語言和FORTRAN語言相比更容易被掌握。通過M語言可以用類似數(shù)學(xué)公式的方式來編寫算法大大降低了程序所需的難度并節(jié)省了時間，從而可把主要的精力集中在算法的構(gòu)思而不是編程上。目

15、前電子計算機已廣泛應(yīng)用于電力系統(tǒng)的分析計算潮流計算是其基本應(yīng)用軟件之一?，F(xiàn)有很多潮流計算方法。對潮流計算方法有五方面的要求（1）計算速度快（2）內(nèi)存需要少（3）計算結(jié)果有良好的可靠性和可信（4）適應(yīng)性好亦即能處理變壓器變比調(diào)整、系統(tǒng)元件的不同描述和與其它程序配合的能力強。4 潮流分布計算4.1 系統(tǒng)的一次接線圖圖4.1 系統(tǒng)的一次連接圖4.2 參數(shù)計算設(shè)定基準(zhǔn)值，Ub=Uav.n，則各參數(shù)如下。（1）發(fā)電機的次暫態(tài)電抗：X=X*Sb/Sn，ZB=UB2/SN發(fā)電機參數(shù) 單位（MW）電廠裝機容量枯水出力比例豐水出力比例豐大有功豐大無功枯大有功枯大無功短路X*''水電站1300.2

16、0.721.000 13.015 6.000 3.718 0.493水電站2720.20.750.400 31.235 14.400 8.924 0.205417水電站3240.20.716.800 10.412 4.800 2.975 0.61625水電站5300.20.721.000 13.015 6.000 3.718 0.493華鑫電廠500.80.840.000 24.790 40.000 24.790 0.2958（2）110KV升壓變壓器的參數(shù)：電阻：R=PK*UN2/ 1000SN2;電抗：X=UK(%)*UN2/ 100 SN;電導(dǎo)：G=P0/ 1000 UN2;電納：B=I

17、0(%)*SN/ 100 UN2;式中UN以KV為單位，SN以MVA為單位，P0、PK以KW為單位。110KV變壓器參數(shù) 變電站名容量X1X0SbR*X*G*B*A16B28C120D63E20G31.5F12.5+31.5水電站130MW、2*20MVA0.82250.5251000.037364 0.581212 0.000199 0.001445 水電站272MW、3*31.5MVA0.344440.333331000.021199 0.369023 0.000280 0.001992 水電站324MW、31.5MVA1.033330.333331000.021199 0.369023

18、0.000280 0.001992 水電站54*7.5MW、2*20MVA1.033330.333331000.037364 0.581212 0.000199 0.001445 火電廠50MW、31.5+40MVA0.482850.26251000.021199 0.369023 0.000280 0.001992 1000.013838 0.290606 0.000379 0.002529 （3）220KV三繞組變壓器的參數(shù)：電阻：R=PK*UN2/ 1000SN2;電抗：X=UK(%)*UN2/ 100 SN;電導(dǎo)：G=P0/ 1000 UN2;電納：B=I0(%)*SN/ 100 UN

19、2220kV變電站參數(shù)表SBR* X* G*B*高壓側(cè)繞組中壓側(cè)繞組低壓側(cè)繞組高壓側(cè)繞組（1）1000.002920 0.116650 0.001650 0.012000 容量120120120中壓側(cè)繞組（2）100-0.002138 0.071667 0.000499 0.003630 電壓22012110.5低壓側(cè)繞組（3）1000.003580 -0.004575 0.000004 0.000027 （4）110KV線路參數(shù)： r=/S； lg ； Xl*=x1*l*SB/U2av.n。110KV線路參數(shù)標(biāo)么值序號線路名稱導(dǎo)線牌號線路長度kmUB(KV)SB（MW)R*X*B*1水電站1

20、水電站2150301151000.047637 0.090737 0.011506 2AB95801151000.199622 0.252250 0.029730 3BC9511151000.002495 0.003153 0.000372 4水電站4C150651151000.103214 0.196597 0.024929 5BD95631151000.157202 0.198647 0.023412 6CD240601151000.059887 0.175123 0.023964 7水電站2C240751151000.074858 0.218904 0.029955 8ED24014.7

21、1151000.014672 0.042905 0.005871 9水電站3D15010.51151000.016673 0.031758 0.004027 10水電站5T節(jié)點705.51151000.018715 0.017758 0.001986 11T節(jié)點D7051151000.017013 0.016144 0.001805 12T節(jié)點G24051151000.004991 0.014594 0.001997 13H F1503.91151000.006193 0.011796 0.001496 14H 008998 0.021066 0.002740 15水

22、電站2L2401301151000.129754 0.379433 0.051921 16LG240101151000.009981 0.029187 0.003994 17D006352 0.012098 0.001534 18火電廠D24011151000.000998 0.002919 0.000399 110KV變電站負(fù)荷參數(shù) 變電站總?cè)萘控S變站比枯站變比豐大有功豐大無功枯大有功枯大無功A160.50.68.0000 2.6295 9.6000 3.1554 B280.50.614.0000 4.6016 16.8000 5.5219 C1200.50.660

23、.0000 19.7210 72.0000 23.6653 D630.50.631.5000 10.3535 37.8000 12.4243 E200.50.610.0000 3.2868 12.0000 3.9442 F440.50.622.0000 7.2311 26.4000 8.6773 G31.50.50.615.7500 5.1768 18.9000 6.2121 各變電站補償電容參數(shù) 站別運行方式ABCDEFGH低壓側(cè)L低壓側(cè)豐大34.609.805.901.92.0枯大34.40005.901.92.04.3 豐大及枯大下地潮流分布情況電壓是衡量電力系統(tǒng)電能質(zhì)量的標(biāo)準(zhǔn)之一。電壓

24、過高或過低，都將對人身及其用電設(shè)備產(chǎn)生重大的影響。保證用戶的電壓接近額定值是電力系統(tǒng)調(diào)度的基本任務(wù)之一。當(dāng)系統(tǒng)的電壓偏離允許值時，電力系統(tǒng)必須應(yīng)用電壓調(diào)節(jié)技術(shù)調(diào)節(jié)系統(tǒng)電壓的大小，使其維持在允許值范圍內(nèi)。本文經(jīng)過手算形成了等值電路圖，并編寫好了程序得出節(jié)點電壓標(biāo)幺值，使其滿足所要求的調(diào)整范圍。我們首先對給定的程序輸入部分作了簡要的分析，程序開始需要我們確定輸入節(jié)點數(shù)、支路數(shù)、平衡母線號、支路參數(shù)矩陣、節(jié)點參數(shù)矩陣。（1）為了保證整個系統(tǒng)潮流計算的完整性，我們把凡具有母線及發(fā)電機處均選作節(jié)點，這樣，我們確定電廠一母線上的發(fā)電機作為平衡節(jié)點，節(jié)點號為，其它機組作為PV節(jié)點，節(jié)點號為，其余節(jié)點均為PQ

25、節(jié)點，節(jié)點號見等值電路圖。（2）確定完節(jié)點及編號后，各條支路也相應(yīng)確定了，我們對各支路參數(shù)進行了計算。根據(jù)所給實際電路圖和題中的已知條件，有以下公式計算各輸電線路的阻抗和對地支路電容的標(biāo)幺值和變壓器的阻抗標(biāo)幺值。4.3.1 該地區(qū)變壓器的有功潮流分布數(shù)據(jù)(1)該圖為豐大潮流模型圖：圖3.3.1(1) 豐大潮流模型圖該運行方式的電壓合理，負(fù)荷分配也均勻，但是有些線路的負(fù)載率偏低。比如水電廠2-L站的負(fù)載率僅僅為21.1% ，M-L站的線路的負(fù)載率也只是8.8%。(2)該圖為枯大潮流模型圖：圖3.3.1(2) 枯大潮流模型圖該運行方式的電壓合理，負(fù)荷分配也均勻，有些線路的負(fù)載率偏低但是有些線路的負(fù)

26、載率則偏高。比如水電廠2-L站的負(fù)載率僅僅為9.4%，M-L站的線路的負(fù)載率也只是 7.6% D-H站的線路卻重載到84%。4.3.2 重、過載負(fù)荷元件統(tǒng)計表類型負(fù)載率實際是在功率額定容量建議線路0.83870.486.2限D(zhuǎn)站負(fù)荷到4MW以下5 設(shè)計心得通過這次課程設(shè)計,我發(fā)現(xiàn)自己有很多不足的地方,如基礎(chǔ)知識掌握不牢固,很多知識點都忘記了,計算速度慢及準(zhǔn)確性低,分析問題能力不夠全面等等。同時，在設(shè)計的過程中遇到很多問題，如怎樣使用WORD的工具,計算公式輸入，畫圖等。明白了有些東西看起來很簡單，但一旦做起來卻需要很多心思，要注意到很多細(xì)節(jié)問題。要做到能好好理解課本的內(nèi)容，一定要認(rèn)認(rèn)真真做一次

27、計算。因此，完成課程設(shè)計使我對課本的內(nèi)容加深了理解?？傮w來說，這次的課程設(shè)計不單在專業(yè)基礎(chǔ)方面反映了我的學(xué)習(xí)還要加倍努力，還在對一些軟件的應(yīng)用需要加強。計算在各種情況下的潮流分布，對于豐大和枯大情況下的潮流分布有了明確的認(rèn)識。在本次課程設(shè)計過程當(dāng)中，鍛煉了自己實際操作分析能力，理論聯(lián)系實際，對運行中的電力系統(tǒng)，通過潮流計算可以預(yù)知各種負(fù)荷變化和網(wǎng)絡(luò)結(jié)構(gòu)的改變會不會危及系統(tǒng)的安全，系統(tǒng)中所有母線的電壓是否在允許的范圍以內(nèi)，系統(tǒng)中各種元件(線路、變壓器等)是否會出現(xiàn)過負(fù)荷，以及可能出現(xiàn)過負(fù)荷時應(yīng)事先采取哪些預(yù)防措施等。同時采用MATLAB進行計算機的計算，在計算時采用特殊算法使得潮流計算的過程更快

28、，效率更高。而采用計算機的運算應(yīng)該是未來的一種趨勢，所以我會學(xué)習(xí)一定的編輯語言如C,C+等，以提高運算準(zhǔn)確性和快速性?？傮w而言，這次的課程設(shè)計對我們運用所學(xué)知識，發(fā)現(xiàn)、提出、分析和解決實際問題、鍛煉實踐能力的考察，使我們更清楚地知道不足之出，從而提高我們。參考文獻1于永源 主編.電力系統(tǒng)分析 湖南師范大學(xué)出版社M.1992年7月2陳珩編.電力系統(tǒng)穩(wěn)態(tài)分析 水利電力出版社M.1995年1月第二版3邱曉燕.劉天琪電力系統(tǒng)分析的計算機算法北京中國電力出版 2004李光琦.電力系統(tǒng)暫態(tài)分析M.北京：水利電力出版社，2002.5  5陸敏政 主編.電力系統(tǒng)習(xí)題集 水利電力出版社M.1

29、990年附錄:程序%潮流計算MATLAB 粗略程序              %creat a new_datat=0;s=0;r=0;w=0;number=input('How many node are there=');% Convert Pq to a new arrayfor ii=1:number    if data(ii,4)=1     

30、60;  t=t+1;            for jj=1:14                new_data1(t,jj)=data(ii,jj);            end;  

31、0;         a(1,t)=ii;            s=s+1;                           

32、0;            %record the number of the PQ node    end;end;%Convert pv to a new arrayfor ii=1:number    if data(ii,4)=2        t=t+1;        

33、;    for jj=1:14                new_data1(t,jj)=data(ii,jj);            end;            a(1,t)=ii;&#

34、160;           r=r+1;                                      &

35、#160; %record the number of the PV node    end;end;%Convert set_v to a new arrayfor ii=1:number    if data(ii,4)=3        t=t+1;            for jj=1:14     

36、;           new_data1(t,jj)=data(ii,jj);            end;            a(1,t)=ii;          &#

37、160; w=w+1;    end;end;%creat a new_data2x,y=size(data2)             for ii=1:x    for jj=1:2        for mm=1:number         

38、   if data2(ii,jj)=a(1,mm)                new_data2(ii,jj)=mm;            end;        end;    end;end;for ii=1

39、:x    for jj=3:14        new_data2(ii,jj)=data2(ii,jj);    end;end;%creat a YY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,number);for ii=1:x   % for jj=1:14        iii=new_

40、data2(ii,1);        jjj=new_data2(ii,2);        if new_data2(ii,5)=2            sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_dat

41、a2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            Y(iii,jjj)=-sub./new_data2(ii,14);            YY(iii,jjj)=sub./new_data2(ii,14);   &

42、#160;        Y(jjj,iii)=-sub/new_data2(ii,14);             YY(jjj,iii)=sub./new_data2(ii,14);            yy(iii,jjj)=(1.-new_data2(ii,14)./(new_data2(

43、ii,14).*new_data2(ii,14).*sub;            yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14).*sub;        else            Y(iii,jjj)=-new_data2(ii,6)./(ne

44、w_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;             YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6

45、).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(ii,7

46、).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)

47、*i;            yy(iii,jjj)=new_data2(ii,8)./2.*i;            yy(jjj,iii)=new_data2(ii,8)./2.*i;        end;    %end;end;for iii=1:number

48、0;   Y(iii,iii)=0;end;%for ii=1:x  %  for jj=1:14    for     iii=1:number        for jj=1:number   % if iii=jj    Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);    

49、60; % end;   end;end;%creat B, Gfor ii=1:number        for jj=1:number            G(ii,jj)= real(Y(ii,jj);             B(ii,jj)= imag(Y(ii,j

50、j);    end;end;%creat Initial_P Initial_Q Initial_Vfor ii=1:(s+r)    set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7)./100;end;for ii=1:s;    set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8)./100;end;for ii=1:r    set_V(ii,1)=new_data1(ii+s,12).*n

51、ew_data1(ii+s,12);%try to modify for sike of correctingend;Initial_p_q_v=set_P;set_Q;set_V;disp(Initial_p_q_v);%creat Initial_e,Initial_ffor ii=1:number-1    e(ii,1)=1;    f(ii,1)=0.0;%change f to test used to be 1.0end;  e(number,1)=new_data1(number,12);  f(n

52、umber,1)=0; % e(64,1)=0.88;%test 118ieee % f(64,1)=0.39395826829394;   % f(14,1)=0;   % e(10,1)=1.045;   %e(11,1)=1.01;    %e(12,1)=1.07;    %e(13,1)=1.09;% Start NEWTOWN CALULATIONfor try_time=1:25%Creat every node consume P Q an

53、d Un=s;m=r;for ii=1:(n+m)    sum1=0;    for jj=1:(n+m+1)        sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1);    end;    p(ii,1)=sum1;end;for ii=1:n

54、    sum2=0;    for jj=1:(n+m+1)        sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1);    end;    q(ii,1)=sum2;end;disp('q=');disp(q);u=zeros(

55、n+m),1);for ii=(n+1):(n+m)    u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;for ii=n+1:(n+m)     extra_u(ii-n),1)=u(ii,1);end;disp('extra_u=');disp(extra_u);sum=p;q;extra_u;disp(sum)disp(s);disp(p);%creat Jacobiandisp(n);disp(m);for ii=1:(n+m)  &#

56、160; for jj=1:(n+m)        if (ii=jj)            PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);             PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);

57、         else            ss=0;            qq=0;             for num=1:(n+m+1)    

58、;            ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);        &

59、#160;   end;                       PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1            PE(ii,jj)=-qq-G(ii,jj

60、).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1        end;    end;end;來源：( - Matlab 潮流計算程序（改進）_小周_新浪博客 copy=3.14159;disp('=copy=')for ii=1:n    for jj=1:m+n        if (ii=jj)    

61、;                 QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1                       QF(ii,jj)=G(ii,j

62、j).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1                   else            ss=0;            qq=0; 

63、                     for num=1:(n+m+1)                ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);    &

64、#160;           qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);            end;                   

65、60;              QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1            QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1      &#

66、160;           end;    end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros(n+m),(n+m);UF=zeros(n+m),(n+m);for ii=n+1:n+m    for jj=1:(n+m)        if (ii=jj)

67、0;                     UE(ii,jj)=0;            UF(ii,jj)=0;        else      

68、60;     ss=0;            qq=0;                       for num=1:(n+m+1)       

69、60;        ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);            end;

70、             UF(ii,jj)=-2.*f(ii,1);            UE(ii,jj)=-2.*e(ii,1);        end;    end;end;for ii=(n+1):(n+m)    for jj=1:(n+m)    extra_UE(ii-n),jj)=UE(ii,jj);    extra_UF(ii-n),jj)=UF(ii,jj);  &#