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

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

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

3、S&（ 1-2）式中，S&為節(jié)點的注入復(fù)功率，是N*1維列矢量；S&為S&的共軛；E diag V&是由節(jié)點電壓的共軛組成的N*N階對角線矩陣。由(1-1 )和(1-2),可得：s&EY&上式就是潮流方程的復(fù)數(shù)形式，是N維的非線性復(fù)數(shù)代數(shù)方程組。將其展開，有:pjQiV& YjV&j=1,2 ,.,N(1-3)j i式中，j i表示所有和i相連的節(jié)點j ,包括j i o將節(jié)點電壓用極坐標(biāo)表示，即令V& V i ,代入式(1-3)中則有：P jQiVi Gj ijBij V jV V Gj ijBijcos ijj sin ij故有：P V j iV G cos ijBjsin ji=1,2,N

4、(1-4)Q V V G sin ijBjcos ijj i式(1-4)是用極坐標(biāo)表示的潮流方程。而節(jié)點功率誤差：P PSP Vi Vj(Gj cos j Bj sin J (1-5)j iQi QiSP VVj(Gij cos ij Bq sin ij)(1-6)j i式中：PiSP，QiSP為節(jié)點i給定的有功功率及無功功率。1.2牛頓-拉夫遜法基本原理牛拉法的一般描述牛拉法是把非線性方程式的求解過程變成反復(fù)對相應(yīng)的線性方程式的求解 過程，即非線性問題通過線性化逐步近似，這就是牛拉法的核心。下面以非線 性方程式的求解過程來進(jìn)行說明。設(shè)電力網(wǎng)絡(luò)的節(jié)點功率方程一般形式如下:SP/、y y x （

5、 1-7）式中，ySP為節(jié)點注入功率給定值；y為ySP對應(yīng)的物理量和節(jié)點電壓之間的函數(shù)表達(dá)式；x為節(jié)點電壓。寫成功率偏差的形式：SPf x y y x 0（ 1-8）應(yīng)用牛拉法求解如下。在給定的初值x 0處將式（1-8）作一階泰勒展開:定義J耳為潮流方程的雅克比矩陣,xJ。為J在x 0處的值，則有:xJ01f x用x修正x 0就得到x的新值。如果用表示迭代次數(shù)，寫成一般的表達(dá)式,有:k 1 f x k（1-9）對于潮流收斂的情況,應(yīng)比xk更接近于解點。收斂條件為:maxf i x k由簡單迭代法收斂性分析的結(jié)論知，越接近解點，牛頓-拉夫遜法收斂越快，它具有二階收斂速度。由圖1.1可以直觀地了解

6、牛拉法的步驟:圖1.1牛頓-拉夫遜法的幾何解釋極坐標(biāo)的牛頓-拉夫遜法在極坐標(biāo)中，f x有如下的形式:,V,VpSPQSP,V,V(1-10)共2n-r個方程，狀態(tài)變量為:T vtn V1 V2 L Vn r共2n-r個待求量。r個PV節(jié)點的電壓幅值給定,不需求解。潮流雅克比矩陣的維數(shù)是(2n-r) *(2n-r)，結(jié)構(gòu)如下:npV7_QV7 n r1，為使雅克比上式右側(cè)的對電壓幅值的偏導(dǎo)數(shù)項中的電壓幅值的階數(shù)減少了 矩陣的各部分子矩陣具有一致的形式，在實際計算中，常將該項乘以電壓幅T值，并選取 V / VV,/ V,V2/ V2LVnr/ Vnr作為待求的修正量，則雅克比矩陣可寫成:將式（1-

7、10）和（1-11）代入式它修正x直到max.將式（1-11）用下式表示:其中每個字塊的計算公式如下:Hh對角元素:非對角元素:NiiMiHijNjMjLj二 VQAVT(1-11)(1-9)為止。的修正方程即可求得x的修正量x，用Bii衛(wèi)VVQjQ-V.V j2牛頓法潮流計算步驟2.1程序流程圖vnM nVMM MiiVLMV Hj Vj, HjVNjVj,NjVMjV，MjVLjV IBiBijQ疋P才P疋Qcos ijGj cosNjHjij(1-12)sinjBj sinij(1-13)在了解了牛拉法的原理之后，明確程序編寫思路，如圖2.1、2.2所示。其中圖2.1中的“計算電壓幅值和

8、角度”步驟較多，單獨用圖2.2表示出來。圖2.1牛頓法計算潮流的程序框圖圖22電壓幅值和角度求解步驟框圖當(dāng)不符合收斂的條件“ amontk1 ”時，即認(rèn)為計算不收斂。具體程序見附錄。2.2計算步驟下面討論的是極坐標(biāo)形式的牛頓法求解過程，大致分為以下幾個步驟： 形成節(jié)點導(dǎo)納矩陣； 給各節(jié)點電壓設(shè)初值（V（0）, i（0）； 根據(jù)式（1-12）、（ 1-13）生成雅克比矩陣（H、N、M、L）; 將節(jié)點電壓初值代入式（1-5）、式（1-6），求出修正方程式的常數(shù)項向量P, Qi ； 求解修正方程，得到電壓幅值和角度； 判斷是否收斂，若收斂，計算平衡節(jié)點和線路功率； 輸出結(jié)果，并結(jié)束。3算例3.1系統(tǒng)

9、模型本文以圖3.1所示電力網(wǎng)絡(luò)為例，調(diào)用基于牛頓-拉夫遜法的C+程序其中節(jié)點4設(shè)為平衡節(jié)點，電壓標(biāo)幺值為1.05,計算誤差為0.0000013.2輸入與輸出將圖3.1所示模型的相關(guān)數(shù)據(jù)放在data.dat文件中$ Efibcsgft /isual C+ + - d牛就拉夫遜搖日上曰文他B 輛吒)童看GO捶入工程0 會 IMffi 京口血 m(H匕 Qy |叮|虛錯|恤|JIdl丄1料1|1.05 &. 90999101 1 2 9.1 0.4 氛 91528? 1 4 9.12 0.5 0.017203 2 4 0.08 0.4 0.0141301130 0.3 fl.909a01 1 0 0

10、 0.3 0.1B2 2 0 0.55 U.133 3 0,5 0 0 Q013 1-19 00圖3.2輸入節(jié)點和支路數(shù)據(jù)對各個數(shù)字含義的解釋如下:網(wǎng)絡(luò)模型有四個節(jié)點，四條支路，編號見圖3.1。第一個零下面三行數(shù)為支路參數(shù)，分別表示三條支路的起始和終止節(jié)點編號，后面的為電阻、電抗和 電納，電導(dǎo)均為0,例如：1 2 0.1 0.4 0.01528第二個零下面的為變壓器支 路，各數(shù)字意義同支路參數(shù)。接下去三行均為節(jié)點參數(shù)，分別表示注入有功功 率和無功功率。調(diào)用text.cpp文件，得到運行結(jié)果，見圖3.3和圖3.4。pSirk * C/- U5er&Ad minjt rdtorDesktopg片、

11、牛和拉夫遜法De b jgtext.exe計昇精 = le-B06導(dǎo)納矩陣為.1-04209-8.24368j -0.S8823S*2-35294jB+3.6&703j-G.4S38Se*1.89107J-0.588235+2.35294J1.069-4.72738j 0*0J 0.48079+2.40385jB+3-66703J00j0-3.33333j Q+0J-0.453858+1.89107J-0-480769 *2-4B385j 0+0j0.934627-4.26159j第1迭代的節(jié)點電壓幅值和角度呦度）0.?9345?0.?763111.11.05-0.00881675-0-107

12、S130.1151360第2迭代的節(jié)點電壓幅值和角度戀度0,9848510.96S0441.1105-0.0070932-0.1104820.1190830第迭代的節(jié)點電壓幅值和角度t弧度S.9846230.9647721.1丄.05-0.0867344-0-1125270.117552 Q第4迭代的節(jié)點電壓幅值和角度那度）fl.9846190.96476?1.11.Q5-0.00871534-0-1125?0.11751圖3.3運行結(jié)果1t C ：MJse rs Ad nn i n i stratoDe sttop返片 牛頓垃天辺法 ibLigltwxt 申 f3,98461?0.?&476

13、71,11.05-0.00871534-0Al25790-117510平衡節(jié)點功率為：9.3788G*jfl.264895線路上的功率為：S=9-2459S1+J-0.S=-B.5 +J-0-029001S=-0434&S63+j-0_136187S2i)=-0.2431&+j0-0110535& =B-J0S=0*J0S=-0.312?49+j-014036&-&.50.097016soa-e*j0S=0*ja531?=B*jaI& =0.0482113 +J0.10-164S(42=0-319671+j0_160255S43=0+j0ls=0+jaPress on# key to cont

14、inue圖3.4運行結(jié)果23.3結(jié)果分析將上述仿真結(jié)果整理為表格 3.1、3.2，其中“ + ”表示節(jié)點i輸出功率給節(jié) 點j，“-”表示節(jié)點j輸出功率給i(縱向為i，橫向為j)。表3.1節(jié)點有功功率輸入與輸出節(jié)點號一 123410+0.245981-0.5-0.0465632-0.2431600-0.31294930.500040.04821430.31967100表3.2節(jié)點無功功率輸入與輸出節(jié)點號123410-0.014708-0.029001-0.13618720.011050500-0.1403630.09701600040.104640.16025500根據(jù)表格計算:節(jié)點 1 有功功

15、率：0+0.245981-0.5-0.046563=-0.300582無功功率：0-0.014708-0.029001-0.136187=-0.179896節(jié)點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é)點的注入有功、無功分別

16、為:P仁0.3，Q1= 0.18 ；P2=0.55，Q2=0.130.3005820.30.3100%0.194%潮流計算誤差：QQ0.180.179896100%0.0578%0.18巳 巳0.5561090.55100% 1.11%0.55QQ20.130.1293095100%0.531%0.13P可見，誤差均在允許范圍內(nèi)線路損耗:S120.245981j 0.014708( 0.24316j 0.0110505)0.002821j 0.0036575S130.5 j0.029001(0.5 j 0.097016)j 0.068015S140.046563j 0.136187(0.048

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

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

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

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

21、ance ;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

22、meG;double *lineConsumeB;ifstream instream;public:/下標(biāo)的轉(zhuǎn)換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);

23、 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 為理想變壓器側(cè)， para2 為歸算阻抗側(cè) void countTranformer(double para,double*G ,dou

24、ble*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=-

25、GIJ/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; /導(dǎo)納矩陣的顯示 void display()coutendl 導(dǎo)

26、納矩陣為： endl;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 ;elsecoutconductancenij+susceptance nijj ;coutendl;coutendl;/顯示任意秩矩陣void displayMatrix(double *matr,int ranki,int rankj)int nij;for(int i=1;i=ranki

27、;i+)for(int j=1;j=rankj;j+)nij=(i-1)*rankj+j;coutmatrnij ; coutendl;/coutendl;/構(gòu)造函數(shù) , 打開文件 powerFlowCalculation(string input) instream.open(input(); assert(instream.is_open();/析構(gòu)函數(shù)，關(guān)閉文件powerFlowCalculation()instream.close();/初始化，根據(jù)節(jié)點數(shù)開辟數(shù)組void initi() nodeNumber=readParameter0; / 讀取節(jié)點數(shù)PVNodeNumber=0;

28、conductance=new doublenodeNumber*nodeNumber; susceptance =new doublenodeNumber*nodeNumber; admittanceAmplitude=new doublenodeNumber*nodeNumber; admittanceAngle=new doublenodeNumber*nodeNumber;PQData=new double(nodeNumber-1)*6; / 開辟 PQData PVData=new double(nodeNumber-1)*5; / 開辟 PVData voltageAmplitu

29、de=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;a

30、=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+) / 導(dǎo)納初始化conductancej=0.0;susceptance j=0.0;admittanceAmplitudej=0.0;admittanceAnglej=0.0;lineConsumeGj=0.0;lineConsumeBj=0.0;for(in

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

32、ty=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)lineD

33、ataSerial+;/ 不同元件的計算/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 計算精度 :readParameter3endl

34、; 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() double vector

35、Angle(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; /將直角坐標(biāo)的導(dǎo)納矩陣轉(zhuǎn)換為極坐標(biāo)形式 void converAdmittanceMatrix()int n;for(int i=1;i=nodeNumber;i+)for(int j=1;

36、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(int

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

38、(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+=admittanceAmplitudeconverIndex(

39、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)*vol

40、tageAmplitudej*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);/

41、end for loop k sum*=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);return

42、 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=-s

43、um;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; b

44、eta=voltageAnglei-voltageAnglek- admittanceAngleconverIndex(i,k);sum+=admittanceAmplitudeconverIndex(i,k)*voltageAmplitudek*sin(beta);sum*=voltageAmplitudei;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=-1*voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplit

45、udej*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+nodeNumbe

46、r*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;n=nodeNumber

47、*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 的下標(biāo)轉(zhuǎn)換int converSerial(int i,int j,int rank=6)int serial;serial=(i-1)*rank+j;return serial;/生成常數(shù)向量void creatConstantVector()int i; /i 為 P、 Q 的下標(biāo)double sum=0.0;double beta=0.0;for(int n=1;n=(nodeNumber-1)*2;n+) /n 為向量的下標(biāo)i=(n+1)/2; /i 為節(jié)點號sum=0.0