電力系統(tǒng)穩(wěn)態(tài)分析-牛頓拉夫遜法教材

1、0引言潮流是配電網絡分析的基礎，用于電網調度、運行分析、操作模擬和設計 規(guī)劃，同時也是電壓優(yōu)化和網絡接線變化所要參考的內容。潮流計算通過數(shù)值 仿真的方法把電力系統(tǒng)的詳細運行情況呈現(xiàn)給工作人員，從而便于研究系統(tǒng)在 給定條件下的穩(wěn)態(tài)運行特點。隨著市場經濟的發(fā)展，經濟利益是企業(yè)十分看重 的，而線損卻是現(xiàn)階段阻礙企業(yè)提高效益的一大因素。及時、準確的潮流計算 結果，可以給出配電網的潮流分布、理論線損及其在網絡中的分布，從而為配 電網的安全經濟運行提供參考。從數(shù)學的角度來看，牛頓-拉夫遜法能有效進行非線性代數(shù)方程組的計算且具有二次收斂的特點，具有收斂快、精度高的特 點，在輸電網中得到廣泛應用。隨著現(xiàn)代計算

2、機技術的發(fā)展，利用編程和相關 軟件，可以更好、更快地實現(xiàn)配電網功能，本文就是結合牛頓-拉夫遜法的基本原理，利用C+程序進行潮流計算，計算結果表明該方法具有良好的收斂性、 可靠性及正確性。1牛頓-拉夫遜法基本介紹1.1潮流方程對于N個節(jié)點的電力網絡（地作為參考節(jié)點不包括在內），如果網絡結構 和元件參數(shù)已知，則網絡方程可表示為：Y0 二卩（1-1）式中，Y為N*N階節(jié)點導納矩陣；V為N*1維節(jié)點電壓列向量；I為N*1維節(jié) 點注入電流列向量。如果不計網絡元件的非線性，也不考慮移相變壓器，則Y為對稱矩陣。電力系統(tǒng)計算中，給定的運行變量是節(jié)點注入功率，而不是節(jié)點注入電 流，這兩者之間有如下關系：詡=S?

3、（ 1-2）式中，S為節(jié)點的注入復功率，是N*1維列矢量；S?為S的共軛；E = diag V?是由節(jié)點電壓的共軛組成的N*N階對角線矩陣。由(1-1 )和(1-2),可得：S?二 EYv上式就是潮流方程的復數(shù)形式，是N維的非線性復數(shù)代數(shù)方程組。將其展開，有:P -jQiYjV j=1,2 ,.,N( 1-3)j電式中，j i表示所有和i相連的節(jié)點j ,包括j二i。將節(jié)點電壓用極坐標表示，即令 V，代入式(1-3)中則有：P - jQi = V -日i 瓦(Gj +jBj V z耳 j專= V2：V(G +jBj cosj - j sin 8j i故有：P =VW V(G cos Qij +

4、Bj sin 壞)j転i=1,2,N(1-4)Q =VV(G sin eij Bj cos q )L.j丟式(1-4)是用極坐標表示的潮流方程。而節(jié)點功率誤差：AR = pSP V 瓦 Vj (Gj cosij + Bj sin 州)(1-5)J -iQ =QiSP VVj(Gj cos Bj sin )(1-6)J -i式中：PSP，QiSP為節(jié)點i給定的有功功率及無功功率。1.2牛頓-拉夫遜法基本原理1.2.1牛拉法的一般描述牛拉法是把非線性方程式的求解過程變成反復對相應的線性方程式的求解 過程，即非線性問題通過線性化逐步近似，這就是牛拉法的核心。下面以非線 性方程式的求解過程來進行說明。

5、設電力網絡的節(jié)點功率方程一般形式如下:(1-7)SPy = y x式中，ySP為節(jié)點注入功率給定值；y為ySP對應的物理量和節(jié)點電壓之間的函 數(shù)表達式；x為節(jié)點電壓。寫成功率偏差的形式：SP(1-8)f x 二 y -y x = 0應用牛拉法求解如下。在給定的初值 x 0處將式（1-8）作一階泰勒展開:f (xC) +T .xAxx 0定義JJ為潮流方程的雅克比矩陣，J。為J在x 0處的值，則有：X 二-J 0Jf x 0用厶X修正x 0就得到x的新值。如果用k表示迭代次數(shù)，寫成一般的表達式， 有：xf） = -J （X（k） f （x（k）（1-9）X k 1 二 y-x k對于潮流收斂的情

6、況，x k 1應比x k更接近于解點。收斂條件為: maxfj （xF 兒 c e由簡單迭代法收斂性分析的結論知，越接近解點，牛頓-拉夫遜法收斂越快，它具有二階收斂速度。由圖1.1可以直觀地了解牛拉法的步驟：圖1.1牛頓-拉夫遜法的幾何解釋1.2.2極坐標的牛頓-拉夫遜法在極坐標中，f x有如下的形式：PSP P(B,V )-SP- - fgj心嚴Qv)lqs” _Qv)_”(1-10)-r共2n-r個方程，狀態(tài)變量為：xTV冋馬III M V2川仏共2n-r個待求量。r個PV節(jié)點的電壓幅值給定，不需求解。潮流雅克比矩陣的 維數(shù)是(2n-r) *(2n-r)，結構如下:；：P : P “Cf

7、_ 胡T 討次T |陋Q陋Q喬一n_r上式右側的對電壓幅值的偏導數(shù)項中的電壓幅值的階數(shù)減少了1，為使雅克比矩陣的各部分子矩陣具有一致的形式，在實際計算中，常將該項乘以電壓幅 值，并選取也V / v = QVi / V1 沁/ V2川 U V,作為待求的修正 量，則雅克比矩陣可寫成：;fT.x：P気 I cAQ 二J；QVn（1-11）x的修正量x，用將式（1-10）和（1-11）代入式（1-9）的修正方程即可求得 它修正x直到maxfj （xf ）J c總為止。I將式（1-11）用下式表示：HiiP,莎 NHMH ii 7+ V2NH對角元素:非對角元素:MiHjNjMij：P，苗V MNii

8、竽V =VLiV, Li:P,HMHij爭 Vj 二VNijVjNij=BiV2V2V2cos 弓j - G_Gj cos 可（1-12）sin片B sin Qi.ijij（1-13）= VMiV,Mj =-N其中每個字塊的計算公式如下:31Q，Lj寸 V =VLijV，Lij =Hij2牛頓法潮流計算步驟2.1、2.2 所示。2.2表示出來。2.1程序流程圖在了解了牛拉法的原理之后，明確程序編寫思路，如圖其中圖2.1中的“計算電壓幅值和角度”步驟較多，單獨用圖圖2.1牛頓法計算潮流的程序框圖圖22電壓幅值和角度求解步驟框圖當不符合收斂的條件“ amontk1”時，即認為計算不收斂。具體程序見

9、附錄。22計算步驟下面討論的是極坐標形式的牛頓法求解過程，大致分為以下幾個步驟： 形成節(jié)點導納矩陣； 給各節(jié)點電壓設初值（V（0）用）； 根據(jù)式（1-12）、（ 1-13）生成雅克比矩陣（H、N、M、L）; 將節(jié)點電壓初值代入式（1-5）、式（1-6），求出修正方程式的常數(shù)項向量 汩，Qi ； 求解修正方程，得到電壓幅值和角度； 判斷是否收斂，若收斂，計算平衡節(jié)點和線路功率； 輸出結果，并結束。3算例3.1系統(tǒng)模型本文以圖3.1所示電力網絡為例，調用基于牛頓-拉夫遜法的C+程序其中節(jié)點4設為平衡節(jié)點，電壓標幺值為1.05,計算誤差為0.0000013.2輸入與輸出將圖3.1所示模型的相關數(shù)據(jù)放

10、在data.dat文件中S 和 ibQsoft拉夫遜法da帖.dat(B文件(E) 輛 萱看世)插入工程 組建(B)工MCD窗口應幫助凹茴)t噩鳥口國宙|鈿JJzll dl4 1.05 O.O0QOO1e1 1 2 0.1。上 0-C1&282 14 6.12 05 0,019203 2 4 C.08 0丄 0.91413 01 1 3 0 0.3 0.9090n1 1 Q 0 0.3 0.182 2 Q 0 氛劉 0.133 3 Q.5 0 Q G61 3 1 ,1 0 00圖3.2輸入節(jié)點和支路數(shù)據(jù)對各個數(shù)字含義的解釋如下:網絡模型有四個節(jié)點，四條支路，編號見圖3.1。第一個零下面三行數(shù)為

11、支路參數(shù)，分別表示三條支路的起始和終止節(jié)點編號，后面的為電阻、電抗和 電納，電導均為0,例如：1 2 0.1 0.4 0.01528第二個零下面的為變壓器支 路，各數(shù)字意義同支路參數(shù)。接下去三行均為節(jié)點參數(shù)，分別表示注入有功功 率和無功功率。調用text.cpp文件，得到運行結果，見圖3.3和圖3.4。r - *C；Use rsAd minist ratorDe s Ictopl里片牛碩?S夫遜法ID亡 bugtext.exe G1 0 0 壓e- :4電:1 :4占心點度 數(shù)節(jié)節(jié)精點 節(jié)平曲導納矩陣為；1.04209-8.243&Bj-0.588235+2.35294j0+3.&6703J-

12、0.453858+1.89107j一陽.588235+2 35294J1.06?-4,72738j 0*0j一048076?*2 40385J0i+3_&6703j0+j0-3.33333 j 關 +陰-0.453658+1.89107j-0_4807G?+2 _!0385 j 0+0j0.934&27-4-26157J第丄迭代的節(jié)點電壓幅值和角度成弧度0993590.?763111.11.05-0.00881675-0-1078130-115138 B第2迭代的節(jié)點電壓幅值和角度弧度3.9848516.9&S644丄.11.65-0-0070932-01.1104820-1190830第M迭

13、代的節(jié)點電壓幅值和角度弧度9.984G230.9&47721.11.05-越.00867344-0.1125270.1175520第4迭代的節(jié)點電壓幅值和角度（弧度G.984G190.9&47671.11.05-0.00871534-0.11257?0.117510圖3.3運行結果1* C:U se rsAd mini strato rDe s kto p童片l牛額螢天遜J 云 Debug t?xt. exeB.9846190-9647671_11-0S-fl.0R871534-0.11257?0.117510平衝節(jié)點功率為s0.367886+Jf-2648?5線路上的功率為匕S=0+J0S=

14、0.245981+j-0.0147088=-0-5+j-a-029001S=-.046563+j-0.136187S(21=-0.24316+j0_0110505S=0+j0S=0+j0S(24=-0-312949+j-0.1403&S=0.5+JB.09?B16S=0+j0S(33=0+j0SO4=0*j0S=0+J0i044S=0.319671+j0_160255S=0*j08=0+j0Press an ke to cont inue圖3.4運行結果23.3結果分析將上述仿真結果整理為表格3.1、3.2，其中“ + ”表示節(jié)點i輸出功率給節(jié)點j，“-”表示節(jié)點j輸出功率給i（縱向為i，橫向

15、為j）。表3.1節(jié)點有功功率輸入與輸出節(jié)點號一 123410+0.245981-0.5-0.046563i2-0.2431600-0.31294930.5000140.04821430.31967100表3.2節(jié)點無功功率輸入與輸出節(jié)點號123410-0.014708-0.029001-0.136187i20.011050500-0.1403630.097016000i40.104640.16025500根據(jù)表格計算:節(jié)點 1 有功功率：0+0.245981-0.5-0.046563=-0.300582無功功率：0-0.014708-0.029001-0.136187=-0.179896節(jié)點

16、2 有功功率：-0.24316+0+0-0.312949=-0.556109無功功率：0.0110505+0+0-0.14036=-0.1293095節(jié)點3有功功率：0.5+0+0+0=0.5無功功率：0.097016+0+0+0=0.097016節(jié)點 4 有功功率：0.0482143+0.319671+0+0=0.3678853無功功率：0.10464+0.160255+0+0=0.264895根據(jù)已知條件，兩個PQ節(jié)點的注入有功、無功分別為：P1=0.3，Q1= 0.18 ；P2=0.55，Q2=0.13潮流計算誤差：碼 0.300582 - 0.3P 一 0.3100%= 0.194%0

17、.18 - 0.179896078100%= 0.0578%P2 _ 0.556109 - 0.55可 _055100% = 1.11%oT0.13 - 0.1293095073100%= 0.531%可見，誤差均在允許范圍內線路損耗：S12= 0.245981 -j 0.014708( -0.24316j 0.0110505) = 0.002821 -j 0.0036575S13= -0.5 - j 0.029001(0.5 j 0.097016) - j 0.068015S14- -0.046563 -j 0.136187(0.0482143 j 0.10464)0.0016513 - j

18、 0.031547S24 - -0.312949 -j 0.14036(0.319671 j 0.160255) = 0.006722 j 0.0198953.4結論通過上面的分析與計算，驗證了程序的正確性。由于編寫過程的不足，線 路損耗沒能直接計算出來，而是需要手算，比較遺憾。程序在運行過程中，需 要區(qū)分三種不同的節(jié)點，這由子程序保證實現(xiàn)。相比于快速分解法，牛拉法程 序較為復雜，但更精確一點，潮流誤差較小。4總結本文基于牛頓-拉夫遜潮流算法的基本原理，利用 C+編程計算了一個4節(jié) 點簡單電力網絡的潮流，并驗證了運行成果，誤差在允許范圍之內。因為牛拉 法計算過程中要不斷生成新的雅各布矩陣，所以

19、相對來說占用內存較多，但收 斂速度快，這在程序運行過程中可以體現(xiàn)出來。本文程序并不是特別實用，因 為真正的電力網絡不可能只有幾個節(jié)點，而且各種電力設備的情況也會復雜很 多，因此程序會變得非常大，占用極大內存。但是，我還是通過這次練習，進 一步鞏固了書本上的理論知識，了解了實際操作的過程步驟。最后，感謝楊偉 老師的悉心指導！參考文獻1 朱紅，趙琦，王慶寶 C+程序設計教程M.北京：清華大學出版社，2010.2 張伯明，陳壽孫，嚴正.高等電力網絡分析M.北京：清華大學出版社， 2007.3 吳明波.牛頓-拉夫遜法在潮流計算中的應用J.內蒙古科技與經濟，2011,21： 111-112，115.4

20、明日科技.VisualC+從入門到精通M.北京：清華大學出版社，2012.5 顧潔，陳章潮，徐禧.一種新的配電網潮流算法-改進牛頓拉夫遜法J.華東電 力，2000，5:10-12.6 朱文強.牛頓-拉夫遜法在配電網中的應用J.水利科技，2004，3:55-56,58.7 劉明波，謝敏，趙維新.大電網最優(yōu)潮流計算M.北京：科學出版社，2010.附錄/*/#includeiostream.h #includemath.h #include #include #include#define nodeNumber readParameter6 using namespace std;/class po

21、werFlowCalculation / 導納矩陣的計算類private:int readDataAmount; / 存放元件參數(shù)的讀取個數(shù)int balanceNodeindex;/ 平衡結點號double balanceNodeVoltage; / 平衡節(jié)點電壓幅值 double balanceNodePowerP; / 平衡節(jié)點功率 double balanceNodePowerQ;double calculationAccuracy; / 計算精度 double readParameter7; / 宏定義 double *conductance;double *susceptance

22、;double *admittanceAmplitude;double *admittanceAngle;double *lineData;int PVNodeNumber; /PV 節(jié)點數(shù)double *jacobiMatrix; / 雅克比矩陣double *PQData; / 矩陣，存放 PQ 節(jié)點的數(shù)據(jù) double *PVData; / 矩陣，存放 PV 節(jié)點的數(shù)據(jù) double *voltageAmplitude; / 電壓幅值 double *voltageAngle; / 電壓角度 double *constantVector; / 常數(shù)向量 double *lineConsu

23、meG;double *lineConsumeB;ifstream instream;public:/下標的轉換int converIndex(int i,int j)int serial;serial=(i-1)*nodeNumber+j;return serial;/線路的計算void countLineBranch(double para,double*G,double*B)double GIJ,BIJ;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4);

24、nij=converIndex(para1,para2);nji=converIndex(para2,para1);nii=converIndex(para1,para1);njj=converIndex(para2,para2);Gnji=-GIJ;Gnij=-GIJ;Bnji=-BIJ;Bnij=-BIJ;Gnii+=GIJ;Gnjj+=GIJ;Bnii+= (BIJ+para5);Bnjj+= (BIJ+para5);/變壓器的計算，規(guī)定 para1 為理想變壓器側， para2 為歸算阻抗側 void countTranformer(double para,double*G ,doub

25、le*B)double GIJ,BIJ,b_i,b_j;int nii,njj,nij,nji;GIJ=para3/(para3*para3+para4*para4);BIJ=-para4/(para3*para3+para4*para4); b_i=BIJ*(1-para5)/para5/para5; b_j=BIJ*(para5-1)/para5;nij=converIndex(para1,para2);nji=converIndex(para2,para1);nii=converIndex(para1,para1);njj=converIndex(para2,para2);Gnij=-G

26、IJ/para5;Bnij=-BIJ/para5;Gnji=-GIJ/para5;Bnji=-BIJ/para5;Gnii+= GIJ/para5;Bnii+= (BIJ/para5+b_i);Gnjj+= GIJ/para5;Bnjj+= (BIJ/para5+b_j);/對地支路的計算void countGroundBranch(double para,double*G ,double*B)int nii;nii=converIndex(para1,para1);Gnii+=para2;Bnii+=para3;/導納矩陣的顯示void display()coutendl 導納矩陣為： en

27、dl;int nij;for(int i=1;i=nodeNumber;i+)for(int j=1;j=nodeNumber;j+)nij=converIndex(i,j);if(susceptance nij0) coutconductancenijsusceptance nijj ;else coutconductancenij+susceptance nijj coutendl;coutendl;/顯示任意秩矩陣void displayMatrix(double *matr,int ranki,int rankj)int nij;for(int i=1;i=ranki;i+)for(i

28、nt j=1;j=rankj;j+)nij=(i-1)*rankj+j;coutmatrnij ;coutendl;/coutendl;/構造函數(shù) , 打開文件 powerFlowCalculation(string inputfilename)instream.open(inputfilename.data(); assert(instream.is_open();/析構函數(shù)，關閉文件powerFlowCalculation()instream.close();/初始化，根據(jù)節(jié)點數(shù)開辟數(shù)組void initi()nodeNumber=readParameter0; / 讀取節(jié)點數(shù)PVNode

29、Number=0;conductance=new doublenodeNumber*nodeNumber; susceptance =new doublenodeNumber*nodeNumber; admittanceAmplitude=new doublenodeNumber*nodeNumber; admittanceAngle=new doublenodeNumber*nodeNumber;PQData=new double(nodeNumber-1)*6; / 開辟 PQDataPVData=new double(nodeNumber-1)*5; / 開辟 PVDatavoltage

30、Amplitude=new doublenodeNumber;voltageAngle=new doublenodeNumber;jacobiMatrix=new doublenodeNumber*nodeNumber*4;constantVector=new doublenodeNumber*2;lineData=new doublenodeNumber*(nodeNumber-1)/2*5;lineConsumeG=new doublenodeNumber*nodeNumber;lineConsumeB=new doublenodeNumber*nodeNumber;for(int a=1

31、;a=nodeNumber*(nodeNumber-1)/2*5;a+)lineDataa=0.0;for(int i=1;i=nodeNumber;i+) / 電壓向量的初始化voltageAmplitudei=1.0;voltageAnglei=0.0;for(int j=1;j=(nodeNumber*nodeNumber);j+) / 導納初始化conductancej=0.0;susceptance j=0.0;admittanceAmplitudej=0.0;admittanceAnglej=0.0;lineConsumeGj=0.0;lineConsumeBj=0.0;for(i

32、nt k=0;k=nodeNumber*nodeNumber*4;k+) / 雅克比矩陣初始化 jacobiMatrixk=0.0;for(int n=1;n=nodeNumber*2;n+) / 常數(shù)向量初始化constantVectorn=0.0;/計算導納矩陣主程序void countAdmittance()int lineDataSerial=1;for(int facility=0;facilityreadParameter0;if (readParameter00.5) break; / 序號為 0 時退出循環(huán)for(int i=1;ireadParameteri;if(facil

33、ity=1)lineDatalineDataSerial=readParameteri;lineDataSerial+;if(facility=2)lineDatalineDataSerial=readParameteri;if(i=5)lineDatalineDataSerial=0.0; lineDataSerial-; lineDatalineDataSerial*=readParameteri; lineDataSerial-;lineDatalineDataSerial*=readParameteri; lineDataSerial+=2;/end for if(i=5)lineDa

34、taSerial+;/end for if(facility=2) /end for loop iswitch(facility) / 不同元件的計算 case 0:initi();balanceNodeindex=(int)readParameter1; balanceNodeVoltage=readParameter2; calculationAccuracy=readParameter3; cout 節(jié)點數(shù) :nodeNumberendl 平衡節(jié)點 :readParameter1endl 平衡節(jié)點電壓 :readParameter2endl 計算精度 :readParameter3end

35、l;break;case 1: countLineBranch(readParameter,conductance,susceptance ); break;case 2:countTranformer(readParameter,conductance,susceptance ); break;case 3:countGroundBranch(readParameter,conductance,susceptance ); break;/end switch/end for loop /end for (facility) /end countAdmittance()/*雅克比矩陣*/dou

36、ble vectorAngle(double a,double b)if(fabs(a)0.0000001 & fabs(b)0.0000001) return 0.0; double alpha;alpha=atan(fabs(b)/fabs(a);if(a=0)alpha=3.1415926-alpha;else if(a=0 & b=0 & b=0) alpha=-alpha;return alpha; /將直角坐標的導納矩陣轉換為極坐標形式 void converAdmittanceMatrix()int n;for(int i=1;i=nodeNumber;i+)for(int j=

37、1;j=nodeNumber;j+) n=converIndex(i,j);admittanceAmplituden=sqrt(conductancen*conductancen+susceptancen*susceptan cen);admittanceAnglen=vectorAngle(conductancen,susceptancen);/讀取PQ節(jié)點和PV節(jié)點的數(shù)據(jù)void readNodePowerData()for(int j=1;jPQDataj;double dustbin; / 將 PQData 與 PVData 之間的 0 去掉instreamdustbin;for(in

38、t i=1;iPVDatai; if(i-1)%5=0) if(PVDatai=0) break;PVNodeNumber+=1;/end for loop read PVData/電壓幅值和角度賦初值void vectorAssignment() /PQ 節(jié)點的值已經在初始化中賦值了，這里給平衡節(jié)點和 PV 節(jié)點賦值int n;voltageAmplitudebalanceNodeindex=balanceNodeV oltage; voltageAnglebalanceNodeindex=0.0;for(int i=1;i=PVNodeNumber;i+)n=PVDataconverInd

39、ex(i,2); voltageAmplituden=PVDataconverIndex(i,3); voltageAnglen=0.0;/生成 L 矩陣double creatL(int i,int j)double sum=0.0;double alpha=0.0;double beta=0.0;if(i=j) alpha=admittanceAngleconverIndex(i,i); for(int k=1;k=nodeNumber;k+) beta=-admittanceAngleconverIndex(i,k); sum+=admittanceAmplitudeconverInde

40、x(i,k)*voltageAmplitudek*sin(beta);/end for loop k sum*=voltageAmplitudei;sum=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex(i,i)*sin(alpha );sum=-sum;elsebeta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*v

41、oltageAmplitudej*sin(beta); return sum;/生成N矩陣double creatN(i nt i,i nt j)double sum=0.0;double alpha=0.0;double beta=0.0;if(i=j)alpha=admittanceAngleconverIndex(i,i);for(int k=1;k=nodeNumber;k+) beta=-admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta

42、); /end for loop ksum*=voltageAmplitudei;sum+=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex(i,i) *cos(alpha);sum=-sum;elsebeta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=- voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*cos(beta);ret

43、urn sum;/生成 J 矩陣double creatJ(int i,int j)double sum=0.0;double beta;if(i=j)for(int k=1;k=nodeNumber;k+)if(k=i) continue; beta=voltageAnglei-voltageAnglek- admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*cos(beta); /end for loop ksum*=voltageAmplitudei; sum

44、=-sum;else beta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej* cos(beta);/end for ifreturn sum;/生成H矩陣double creatH(i nt i,i nt j)double sum=0.0; double beta;if(i=j) for(int k=1;k=nodeNumber;k+)if(k=i) continue

45、; beta=voltageAnglei-voltageAnglek- admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*sin(beta);sum*=voltageAmplitudei;else beta=voltageAnglei-voltageAnglej- admittanceAngleconverIndex(i,j);sum=- 1*voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltage

46、Amplitudej*sin(beta);/end for ifreturn sum;/生成雅各比矩陣void creatJacobiMatrix()int serial;int serialH;int serialN;int serialJ;int serialL;for(int i=1;i=2*(nodeNumber-1);i=i+2)for(int j=1;j=2*(nodeNumber-1);j=j+2)serialH=(i-1)*nodeNumber*2+j; serialN=(i-1)*nodeNumber*2+j+1; serialJ=(i-1)*nodeNumber*2+j+n

47、odeNumber*2; serialL=(i-1)*nodeNumber*2+j+nodeNumber*2+1; jacobiMatrixserialH=creatH(i+1)/2,(j+1)/2); jacobiMatrixserialN=creatN(i+1)/2,(j+1)/2); jacobiMatrixserialJ=creatJ(i+1)/2,(j+1)/2); jacobiMatrixserialL=creatL(i+1)/2,(j+1)/2);/end for loop j/end for loop ifor(int n=(nodeNumber-PVNodeNumber)*2

48、;n=nodeNumber*2;n=n+2)for(int m=1;m=nodeNumber*2;m+)if(n=m)serial=(m-1)*nodeNumber*2+n; jacobiMatrixserial=1.0;elseserial=(m-1)*nodeNumber*2+n; jacobiMatrixserial=0.0; serial=(n-1)*nodeNumber*2+m; jacobiMatrixserial=0.0;/end for if/end for loop m/end for loop n serial=nodeNumber*nodeNumber*4-nodeNumber*2-1; jacobiMatrixserial=1.0;/ rank 為 6 的下標轉換int converSerial(int i,int j,int rank=6)int serial;serial=(i-1)*rank+j;return serial;/生成常數(shù)向量void creatConstantVector(