1、Companynumber:0089WT-8898YT-W8CCB-BUUT-202108電力系統課程設計題目：電力系統潮流計算院系名稱：電氣工程學院專業班級：電氣F1206班學生姓名：學號：指導教師：張孝遠成績：目錄原始資料指導老師簽名-21概述日期:42潮流計算節點介紹4變量的分類5節點的分類53計算方法簡介6牛頓拉夫遜法原理6牛頓拉夫遜法概要6牛頓法的框圖及求解過程8MATLAB簡介94潮流分布計算10系統的一次接線圖10參數計算10豐大及枯大下地潮流分布情況14該地區變壓器的有功潮流分布數據15重、過載負荷元件統計表175設計心得17參考文獻18附錄:程序19原始資料一、系統接線圖見附

2、件1。二、系統中包含發電廠、變電站、及其間的聯絡線路。500kV變電站以外的系統以一個等值發電機代替。各元件的參數見附件2。設計任務1、手動畫出該系統的電氣一次接線圖，建立實際網絡和模擬網絡之間的聯系。2、根據已有資料，先手算出各元件的參數，后再用Matlab表格核算出各元件的參數。3、潮流計算1）對兩種不同運行方式進行潮流計算，注意110kV電網開環運行。2）注意將電壓調整到合理的范圍110kV母線電壓控制在106kV117kV之間；220kV母線電壓控制在220kV242kV之間。附件一：A"一一一:二GEFL：x150H12.5+31.53x40D.2x31.5課程設計地理接線

3、示意圖1x31.590+120火電廠110kV線路B20+8水電站1302x10O-2x80水電站272水電站3口24E3水電站火電廠鬥水電站54臼牽引站曰I-220kV線路llOkV變電站_220kV變電站附件二：1、變壓器：兩個220kV變電站均采用參數一致的三繞組變壓器，具體參數如下。220kV變電站參數表咼壓側繞組中壓側繞組低壓側繞組容量12012060(120)電壓220110(121)110kV及以下的變電站的變壓器省略，即可將負荷直接掛在110kV母線上。而110kV升壓變只計及以下參數。110kV變電站參數表序號變電站名容量X1X01A162B283C1205D636E207G

4、8F+9水電站130MW、2*20MVA10水電站272MW、3*11水電站324MW、12水電站54*、2*20MVA13水電站418MW、30MVA14火電50MW廠、+40MVA2、線路：具體參數如下。220kV線路參數表序號線路名稱導線牌號線路長度km1ML2x240102ML2x240103MH2x3005110kV線路參數表序號線路名稱導線牌號線路長度km1水電站1水電站2150302AB95803BC9514水電站4C150655BD95636CD240607水電站2C240758ED2409水電站3D15010水電站5T節點7011T節點D70512T節點G240513HF15

5、014HD185715水電站2L240/213016LG240103、發電機各發電機的參數如下：XdttXqtt裝機容量功率因數水電站130水電站272水電站324水電站530水電站418華鑫電廠5出力情況：水力發電機豐大出力70%,枯大出力20%火力發電機豐大出力80%,枯大出力80%。4、負荷各110kV變電站豐大負荷按該站變電容量的50%估算，枯大負荷按該站變電容量的60%估算。兩個220kV變電站的低壓側上各掛10MW的負荷，中壓側各掛20MW負荷。功率因素均為。5、并聯電容器兩個220kV變電站的低壓側上均裝設并聯補償。補償總量按該站變電容量的20%裝設，分組原則以每組電容器的容量不

6、超過10MVar且經濟性較好為準。1 概述潮流計算是電力系統最基本最常用的計算。根據系統給定的運行條件，網絡接線及元件參數，通過潮流計算可以確定各母線的電壓，包括電壓的幅值和相角，各元件流過的功率，整個系統的功率損耗等一系列系統中的潮流數據。近幾年，對潮流算法的研究仍然是如何改善傳統的潮流算法，即高斯-塞德爾法、牛頓法和快速解耦法。牛頓法，由于其在求解非線性潮流方程時采用的是逐次線性化的方法，為了進一步提高算法的收斂性和計算速度，人們考慮采用將泰勒級數的高階項或非線性項也考慮進來，于是產生了二階潮流算法。后來又提出了根據直角坐標形式的潮流方程是一個二次代數方程的特點，提出了采用直角坐標的保留非

7、線性快速潮流算法。潮流計算在數學上是多元非線性方程組的求解問題，求解的方法有很多種，牛頓拉夫遜Newton-Raphson法是數學上解非線性方程組的有效方法，有較好的收斂性。將N-R法用于潮流計算是以導納矩陣為基礎的，由于利用了導納矩陣的對稱性，稀疏性及節點編號順序優劃等技巧，使N-R法在收斂性，占用內存，計算速度等方面的優點都超過了阻抗法總結為在電力系統運行方式和規劃方案的研究中都需要進行潮流計算以比較運行方式或規劃供電方案的可行性、可靠性和經濟性。同時為了實時監控電力系統的運行狀態也需要進行大量而快速的潮流計算。因此潮流計算是電力系統中應用最廣泛、最基本和最重要的一種電氣運算。在系統規劃設

8、計和安排系統的運行方式時采用離線潮流計算在電力系統運行狀態的實時監控中則采用在線潮流計算。2 潮流計算節點介紹常規的電力系統潮流計算中一般具有三種類型的節點：PQ、PV及平衡節點。一個節點有四個變量，即注入有功功率、注入無功功率，電壓大小及相角。常規的潮流計算一般給定其中的二個變量:PQ節點（注入有功功率及無功功率），PV節點（注入有功功率及電壓的大小），平衡節點（電壓的大小及相角）。變量的分類負荷消耗的有功、無功功率一一p、Q、p、QL1L1L2L2電源發出的有功、無功功率一一P、Q、p、QG1G1G2G2母線或節點的電壓大小和相位U、U、5、81212在這十二個變量中，負荷消耗的有功和無功

9、功率無法控制，因它們取決于用戶，它們就稱為不可控變量或是擾動變量。電源發出的有功無功功率是可以控制的自變量，因此它們就稱為控制變量。母線或節點電壓的大小和相位角一是受控制變量控制的因變量。其中，U、U主要受Q、Q的控制，5、12G1G215主要受p、p的控制。這四個變量就是簡單系統的狀態變量。2G1G2為了保證系統的正常運行必須滿足以下的約束條件：對控制變量對沒有電源的節點則為對狀態變量U的約束條件則是i對某些狀態變量5還有如下的約束條件i節點的分類第一類稱PQ節點。等值負荷功率P、Q和等值電源功率P、Q是LiLiGiGi給定的，從而注入功率P、Q是給定的，待求的則是節點電壓的大小U和相iii

10、位角5。屬于這類節點的有按給定有功、無功率發電的發電廠母線和沒有其他i電源的變電所母線。(第二類稱PV節點。等值負荷和等值電源的有功功率P、P是給定LiGi的，從而注入有功功率P是給定的。等值負荷的無功功率Q和節點電壓的大小iLiU也是給定的。待求的則是等值電源的無功功率Q，從而注入無功功率QiGii和節點電壓的相位角5。有一定無功功率儲備的發電廠和有一定無功功率電源i的變電所母線都可以作為PV節點；第三類平衡節點。潮流計算時一般只設一個平衡節點。等值負荷功率P、Q是給定的，節點電壓的大小和相位也是給定的。擔負調整系統頻率任LsLs務的發電廠母線往往被選作為平衡節點。3 計算方法簡介牛頓拉夫遜

11、法原理牛頓拉夫遜法概要首先對一般的牛頓一拉夫遜法作一簡單的說明。已知一個變量X函數為：到此方程時，由適當的近似值X出發，根據：反復進行計算，當X(n)滿足適當的收斂條件就是上面方程的根。這樣的方法就是所謂的牛頓拉夫遜法這一方法還可以做下面的解釋，設第n次迭代得到的解語真值之差，即X(n)的誤差為£時，貝,把f(X(n)+£)在X(n)附近對£用泰勒級數展開上式省略去£2以后部分X(n)的誤差可以近似由上式計算出來。比較兩式，可以看出牛頓一拉夫遜法的休整量和X(n)的誤差的一次項相等。用同樣的方法考慮，給出n個變量的n個方程：對其近似解X'得修正量

12、AX'可以通過解下邊的方程來確定：11式中等號右邊的矩陣%都是對于X:X'，,X'的值。這一矩陣稱為雅可dx12nn比(JACOBI)矩陣。按上述得到的修正向量AX:AX'，,AX'后，得到如下關12n系這比X:X'，,X'更接近真實值。這一步在收斂到希望的值以前重復進12n行，一般要反復計算滿足£為預先規定的小正數，Xn+1是第n次迭代X的近似值。nn牛頓法的框圖及求解過程1、用牛頓法計算潮流時，有以下的步驟：(1) 給這各節點電壓初始值e(0),f；(2) 將以上電壓初始值代入公式，求修正方程的常數項向量AP(0),AQ(0

13、),(AV2)(0)；（3）將電壓初始值在帶入上述公式，求出修正方程中系數矩陣的各元素。（4）解修正方程式2（o）,曲（o）；（5）修正各節點電壓e=e（o）+Ae（o），/=/+Af；（6）將e，/在帶入方程式，求出AP（i）,AQ（i）,（AV2）（1）；（7）檢驗是否收斂，即maxIp（k）,AQ（K）LsII（8）如果收斂，迭代到此結束，進一步計算各線路潮流和平衡節點功率，并打印輸出結果。如果不收斂，轉回（2）進行下次迭代計算，直到收斂為止。2、程序框圖如下MATLAB簡介MATLAB是用于算法開發、數據可視化、數據分析以及數值計算的高級技術計算語言和交互式環境，主要包括MATLAB和

14、Simulink兩大部分。是由美國mathworks公司發布的主要面對科學計算、可視化以及交互式程序設計的高科技計算環境。它將數值分析、矩陣計算、科學數據可視化以及非線性動態系統的建模和仿真等諸多強大功能集成在一個易于使用的視窗環境中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一種全面的解決方案，并在很大程度上擺脫了傳統非交互式程序設計語言（如C、Fortran）的編輯模式，代表了當今國際科學計算軟件的先進水平。MATLAB是一種交互式、面向對象的程序設計語言廣泛應用于工業界與學術界主要用于矩陣運算同時在數值分析、自動控制模擬、數字信號處理、動態分析、繪圖等方面也具有強大

15、的功能。MATLAB程序設計語言結構完整且具有優良的移植性它的基本數據元素是不需要定義的數組。它可以高效率地解決工業計算問題特別是關于矩陣和矢量的計算。MATLAB與C語言和FORTRAN語言相比更容易被掌握。通過M語言可以用類似數學公式的方式來編寫算法大大降低了程序所需的難度并節省了時間，從而可把主要的精力集中在算法的構思而不是編程上。目前電子計算機已廣泛應用于電力系統的分析計算潮流計算是其基本應用軟件之一。現有很多潮流計算方法。對潮流計算方法有五方面的要求（1）計算速度快（2）內存需要少（3）計算結果有良好的可靠性和可信（4）適應性好亦即能處理變壓器變比調整、系統元件的不同描述和與其它程序

16、配合的能力強。潮:系統的一次接線圖圖系統的一次連接圖參數計算設定基準值5100MVA，Ub=，則各參數如下。B（1）發電機的次暫態電抗：X=X*Sb/Sn,Zb=Ub2/Sn發電機參數單位（MW）電廠裝機容量枯水出力比例豐水出力比例豐大有功豐大無功枯大有功枯大無功短路X*''水電站130水電72電阻：R=Pk*Un2/1000Sn2;電抗：X=UK(%)*UN2/100SN;電導：G=P0/1000UN2;電納：B=I0(%)*SN/100UN2;式中Un以KV為單位，Sn以MVA為單位，p0、Pk以KW為單位110KV變壓器參數變電站名容量XIX0SbR*X*G*B*A16B

17、28C120D63E20GF+水電站130MW、2*20MVA100水電站272MW、3*100水電站324MW、100水電站54*、2*20MVA100火電廠50MW、+40MVA100100220KV三繞組變壓器的參數：電阻：R=Pk*Un2/1000Sn2;電抗：X=UK(%)*UN2/100SN;電導：G=P0/1000UN2;電納：B=I0(%)*SN/100UN2220kV變電站參數表SBR*X*G*B*高壓側繞組中壓側繞組低壓側繞組高壓側繞組(1)100容量120120120中壓側繞組(2)100電壓220121低壓側繞組(3)1004)110KV線路參數：r二p/$二=5；X1

18、*=x1*1*SB/。110KV線路參數標么值序號線路名稱導線牌號線路長度kmUB(KV)SB(MW)R*X*B*1水電站1水電站2150301151002AB95801151003BC9511151004水電站4C150651151005BD95631151006CD240601151007水電站2C240751151008ED2401151009水電站3電站5T節點7011510011T節點D70511510012T節點G240511510013HD185711510015水電站2L24013011510016LG240101151001

19、7DF150411510018火電D2401115100110KV變電站負荷參數變電站總容曰豐變各站卜卜比變枯電站泳-A瑙償電容參豐大無T功枯大有T功枯大無功量站匸比變比功功A16B28C120D63E20F44G站別運行方式ABCDEFGH低壓側L低壓側豐大3000枯大30000豐大及枯大下地潮流分布情況電壓是衡量電力系統電能質量的標準之一。電壓過高或過低，都將對人身及其用電設備產生重大的影響。保證用戶的電壓接近額定值是電力系統調度的基本任務之一。當系統的電壓偏離允許值時，電力系統必須應用電壓調節技術調節系統電壓的大小，使其維持在允許值范圍內。本文經過手算形成了等值電路圖，并編寫好了程序得出

20、節點電壓標幺值，使其滿足所要求的調整范圍。我們首先對給定的程序輸入部分作了簡要的分析，程序開始需要我們確定輸入節點數、支路數、平衡母線號、支路參數矩陣、節點參數矩陣。（1）為了保證整個系統潮流計算的完整性，我們把凡具有母線及發電機處均選作節點，這樣，我們確定電廠一母線上的發電機作為平衡節點，節點號為，其它機組作為PV節點，節點號為，其余節點均為PQ節點，節點號見等值電路圖。（2）確定完節點及編號后，各條支路也相應確定了，我們對各支路參數進行了計算。根據所給實際電路圖和題中的已知條件，有以下公式計算各輸電線路的阻抗和對地支路電容的標幺值和變壓器的阻抗標幺值。該地區變壓器的有功潮流分布數據（1）該

21、圖為豐大潮流模型圖：豐大潮流模型圖該運行方式的電壓合理，負荷分配也均勻，但是有些線路的負載率偏低。比如水電廠2L站的負載率僅僅為%，ML站的線路的負載率也只是%。(2)該圖為枯大潮流模型圖：)枯大潮流模型圖該運行方式的電壓合理，負荷分配也均勻，有些線路的負載率偏低但是有些線路的負載率則偏高。比如水電廠2L站的負載率僅僅為,ML站的線路的負載率也只是%D-H站的線路卻重載到84%。重、過載負荷元件統計表類型負載率實際是在功率額定容量建議線路限D站負荷到4MW以下5設計心得通過這次課程設計,我發現自己有很多不足的地方,如基礎知識掌握不牢固,很多知識點都忘記了,計算速度慢及準確性低,分析問題能力不夠

22、全面等等。同時,在設計的過程中遇到很多問題,如怎樣使用WORD的工具,計算公式輸入,畫圖等。明白了有些東西看起來很簡單,但一旦做起來卻需要很多心思,要注意到很多細節問題。要做到能好好理解課本的內容,一定要認認真真做一次計算。因此,完成課程設計使我對課本的內容加深了理解。總體來說,這次的課程設計不單在專業基礎方面反映了我的學習還要加倍努力，還在對一些軟件的應用需要加強。計算在各種情況下的潮流分布，對于豐大和枯大情況下的潮流分布有了明確的認識。在本次課程設計過程當中，鍛煉了自己實際操作分析能力，理論聯系實際，對運行中的電力系統，通過潮流計算可以預知各種負荷變化和網絡結構的改變會不會危及系統的安全，

23、系統中所有母線的電壓是否在允許的范圍以內，系統中各種元件（線路、變壓器等）是否會出現過負荷，以及可能出現過負荷時應事先采取哪些預防措施等。同時采用MATLAB進行計算機的計算，在計算時采用特殊算法使得潮流計算的過程更快，效率更高。而采用計算機的運算應該是未來的一種趨勢，所以我會學習一定的編輯語言如C,C+等，以提高運算準確性和快速性。總體而言，這次的課程設計對我們運用所學知識，發現、提出、分析和解決實際問題、鍛煉實踐能力的考察，使我們更清楚地知道不足之出，從而提高我們。參考文獻1 于永源主編電力系統分析湖南師范大學出版社M.1992年7月2 陳珩編電力系統穩態分析水利電力出版社M.1995年1

24、月第二版3 邱曉燕劉天琪電力系統分析的計算機算法北京中國電力出版2004 李光琦電力系統暫態分析M.北京：水利電力出版社，5 陸敏政主編電力系統習題集水利電力出版社M.199O年附錄:程序%潮流計算MATLAB粗略程序？%creatanew_datat=0;s=0;r=0;w=0;number=input('Howmanynodearethere=');%ConvertPqtoanewarrayforii=1:numberifdata(ii,4)=1t=t+1;forjj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;s=s+1;%

25、recordthenumberofthePQnodeend;end;%Convertpvtoanewarrayforii=1:numberifdata(ii,4)=2t=t+1;forjj=1:14new_data1(t,jj)=data(ii,jj);a(1,t)=ii;r=r+1;%recordthenumberofthePVnodeend;end;%Convertset_vtoanewarrayforii=1:numberifdata(ii,4)=3t=t+1;forjj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;w=w+1;end;en

26、d;%creatanew_data2x,y=size(data2)forii=1:xforjj=1:2formm=1:numberifdata2(ii,jj)=a(1,mm)new_data2(ii,jj)=mm;end;end;end;end;forii=1:xforjj=3:14new_data2(ii,jj)=data2(ii,jj);end;end;%creataYY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,number);forii=1:x%forjj=1:14iii=new_data2(ii,1);j

27、jj=new_data2(ii,2);ifnew_data2(ii,5)=2sub=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)*i;Y(iii,jjj)=-sub./new_data2(ii,14);YY(iii,jjj)=sub./new_data2(ii,14);Y(jjj,iii)=-sub/new_

28、data2(ii,14);YY(jjj,iii)=sub./new_data2(ii,14);yy(iii,jjj)=(ii,14)./(new_data2(ii,14).*new_data2(ii,14).*sub;yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14).*sub;elseY(iii,jjj)=-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).

29、*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).*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

30、_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(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)*i;yy

31、(iii,jjj)=new_data2(ii,8)./2.*i;yy(jjj,iii)=new_data2(ii,8)./2.*i;end;%end;end;foriii=1:numberY(iii,iii)=0;end;%forii=1:x%forjj=1:14foriii=1:numberforjj=1:number%ifiii=jjY(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);%end;end;end;%creatB,Gforii=1:numberforjj=1:numberG(ii,jj)=real(Y(ii,jj);B(ii,jj)=imag

32、(Y(ii,jj);end;end;%creatInitial_PInitial_QInitial_Vforii=1:(s+r)set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7)./100;end;forii=1:s;set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8)./100;end;forii=1:rset_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%trytomodifyforsikeofcorrectingend;Initial_p_q_v=set_P;se

33、t_Q;set_V;disp(Initial_p_q_v);%creatInitial_e,Initial_fforii=1:number-1e(ii,1)=1;f(ii,1)=;%changeftotestusedtobeend;e(number,1)=new_data1(number,12);f(number,1)=0;%e(64,1)=;%test118ieee%f(14,1)=0;%e(10,1)=;%e(11,1)=;%e(12,1)=;%e(13,1)=;%StartNEWTOWNCALULATIONfortry_time=1:25%CreateverynodeconsumePQa

34、ndUn=s;m=r;forii=1:(n+m)sum1=0;forjj=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;forii=1:nsum2=0;forjj=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

35、);end;q(ii,1)=sum2;end;disp('q=');disp(q);u=zeros(n+m),1);forii=(n+1):(n+m)u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;forii=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);%creatJacobiandisp(n);disp(m);forii=1:(n+m)forj

36、j=1:(n+m)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);elsess=0;qq=0;fornum=1:(n+m+1)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;PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1PE(ii,jj)=-qq-G

37、(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1end;end;end;copy=;disp('=copy=')forii=1:nforjj=1:m+nif(ii=jj)QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1elsess=0;qq=0;fornum=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G

38、(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1end;end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros(n+m),(n+m);UF=zeros(n+m),(n+m);forii=n+1:n+mforjj=1:(n+m)UE(ii,

39、jj)=0;UF(ii,jj)=0;elsess=0;qq=0;fornum=1:(n+m+1)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;UF(ii,jj)=-2.*f(ii,1);UE(ii,jj)=-2.*e(ii,1);end;end;end;forii=(n+1):(n+m)forjj=1:(n+m)extra_UE(ii-n),jj)=UE(ii,jj);extra_UF(ii-n),jj)=UF(ii,jj);end;end;%dis

40、p('extra_UE');%disp(extra_UE);%disp('extra_Uf');%disp(extra_UF);Jacobian=PF,PE;QF,QE;extra_UF,extra_UE;%disp('Jacobian=');%disp(Jacobian);%creatsubstractresultsubstract_result=Initial_p_q_v-sum;%disp('substract_result');%disp(substract_result);%calculatedelta_f_edelta

41、_f_e=-inv(Jacobian)*substract_result;%disp(delta_f_e);forii=1:number-1;f(ii,1)=f(ii,1)+delta_f_e(ii,1);e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1);end;ifmax(substract_result)<1e-4break;end;end;%disp('substract_result');%disp(substract_result);%disp('e=');%disp(e);%disp('f=');

42、%disp(f);forii=1:numberuuu(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1);end;forii=1:numberforjj=1:numberifii=a(1,jj)Old_Uresult(ii,1)=U_RESULT(jj,1)end;end;end;forii=1:numberOld_Uresult(ii,2)=ii;end;%disp('U_result');%disp(U_RESULT);disp('=');disp('Thelastresultis:'