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

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

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

3、S&（ 1-2）式中，S&為節(jié)點(diǎn)的注入復(fù)功率，是N*1維列矢量；S&為S&的共軛；E diag V&是由節(jié)點(diǎn)電壓的共軛組成的N*N階對(duì)角線矩陣。由(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é)點(diǎn)j ,包括j i o將節(jié)點(diǎn)電壓用極坐標(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é)點(diǎn)功率誤差：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é)點(diǎn)i給定的有功功率及無(wú)功功率。1.2牛頓-拉夫遜法基本原理牛拉法的一般描述牛拉法是把非線性方程式的求解過(guò)程變成反復(fù)對(duì)相應(yīng)的線性方程式的求解 過(guò)程，即非線性問(wèn)題通過(guò)線性化逐步近似，這就是牛拉法的核心。下面以非線 性方程式的求解過(guò)程來(lái)進(jìn)行說(shuō)明。設(shè)電力網(wǎng)絡(luò)的節(jié)點(diǎn)功率方程一般形式如下:SP/、y y x （

5、 1-7）式中，ySP為節(jié)點(diǎn)注入功率給定值；y為ySP對(duì)應(yīng)的物理量和節(jié)點(diǎn)電壓之間的函數(shù)表達(dá)式；x為節(jié)點(diǎn)電壓。寫成功率偏差的形式：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）對(duì)于潮流收斂的情況,應(yīng)比xk更接近于解點(diǎn)。收斂條件為:maxf i x k由簡(jiǎn)單迭代法收斂性分析的結(jié)論知，越接近解點(diǎn)，牛頓-拉夫遜法收斂越快，它具有二階收斂速度。由圖1.1可以直觀地了解

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

7、10）和（1-11）代入式它修正x直到max.將式（1-11）用下式表示:其中每個(gè)字塊的計(jì)算公式如下:Hh對(duì)角元素:非對(duì)角元素:NiiMiHijNjMjLj二 VQAVT(1-11)(1-9)為止。的修正方程即可求得x的修正量x，用Bii衛(wèi)VVQjQ-V.V j2牛頓法潮流計(jì)算步驟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中的“計(jì)算電壓幅值和

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

9、模型本文以圖3.1所示電力網(wǎng)絡(luò)為例，調(diào)用基于牛頓-拉夫遜法的C+程序其中節(jié)點(diǎn)4設(shè)為平衡節(jié)點(diǎn)，電壓標(biāo)幺值為1.05,計(jì)算誤差為0.0000013.2輸入與輸出將圖3.1所示模型的相關(guān)數(shù)據(jù)放在data.dat文件中$ Efibcsgft /isual C+ + - d牛就拉夫遜搖日上曰文他B 輛吒)童看GO捶入工程0 會(huì) IMffi 京口血 m(H匕 Qy |叮|虛錯(cuò)|恤|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é)點(diǎn)和支路數(shù)據(jù)對(duì)各個(gè)數(shù)字含義的解釋如下:網(wǎng)絡(luò)模型有四個(gè)節(jié)點(diǎn)，四條支路，編號(hào)見圖3.1。第一個(gè)零下面三行數(shù)為支路參數(shù)，分別表示三條支路的起始和終止節(jié)點(diǎn)編號(hào)，后面的為電阻、電抗和 電納，電導(dǎo)均為0,例如：1 2 0.1 0.4 0.01528第二個(gè)零下面的為變壓器支 路，各數(shù)字意義同支路參數(shù)。接下去三行均為節(jié)點(diǎn)參數(shù)，分別表示注入有功功 率和無(wú)功功率。調(diào)用text.cpp文件，得到運(yùn)行結(jié)果，見圖3.3和圖3.4。pSirk * C/- U5er&Ad minjt rdtorDesktopg片、

11、牛和拉夫遜法De b jgtext.exe計(jì)昇精 = 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é)點(diǎn)電壓幅值和角度呦度）0.?9345?0.?763111.11.05-0.00881675-0-107

12、S130.1151360第2迭代的節(jié)點(diǎn)電壓幅值和角度戀度0,9848510.96S0441.1105-0.0070932-0.1104820.1190830第迭代的節(jié)點(diǎn)電壓幅值和角度t弧度S.9846230.9647721.1丄.05-0.0867344-0-1125270.117552 Q第4迭代的節(jié)點(diǎn)電壓幅值和角度那度）fl.9846190.96476?1.11.Q5-0.00871534-0-1125?0.11751圖3.3運(yùn)行結(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é)點(diǎn)功率為：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運(yùn)行結(jié)果23.3結(jié)果分析將上述仿真結(jié)果整理為表格 3.1、3.2，其中“ + ”表示節(jié)點(diǎn)i輸出功率給節(jié) 點(diǎn)j，“-”表示節(jié)點(diǎn)j輸出功率給i(縱向?yàn)閕，橫向?yàn)閖)。表3.1節(jié)點(diǎn)有功功率輸入與輸出節(jié)點(diǎn)號(hào)一 123410+0.245981-0.5-0.0465632-0.2431600-0.31294930.500040.04821430.31967100表3.2節(jié)點(diǎn)無(wú)功功率輸入與輸出節(jié)點(diǎn)號(hào)123410-0.014708-0.029001-0.13618720.011050500-0.1403630.09701600040.104640.16025500根據(jù)表格計(jì)算:節(jié)點(diǎn) 1 有功功

15、率：0+0.245981-0.5-0.046563=-0.300582無(wú)功功率：0-0.014708-0.029001-0.136187=-0.179896節(jié)點(diǎn)2有功功率：-0.24316+0+0-0.312949=-0.556109無(wú)功功率：0.0110505+0+0-0.14036=-0.1293095節(jié)點(diǎn)3有功功率：0.5+0+0+0=0.5無(wú)功功率：0.097016+0+0+0=0.097016節(jié)點(diǎn)4有功功率：0.0482143+0.319671+0+0=0.3678853無(wú)功功率:0.10464+0.160255+0+0=0.264895根據(jù)已知條件，兩個(gè)PQ節(jié)點(diǎn)的注入有功、無(wú)功分別

16、為:P仁0.3，Q1= 0.18 ；P2=0.55，Q2=0.130.3005820.30.3100%0.194%潮流計(jì)算誤差：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é)論通過(guò)上面的分析與計(jì)算，驗(yàn)證了程序的正確性。由于編寫過(guò)程的不足，線路損耗沒能直接計(jì)算出來(lái)，而是需要手算，比較遺憾。程序在運(yùn)行過(guò)程中，需 要區(qū)分三種不同的節(jié)點(diǎn)，這由子程序保證實(shí)現(xiàn)。相比于快速分解法，牛拉法程 序較為復(fù)雜，但更精確一點(diǎn)，潮流誤差較小。4總結(jié)本文基于牛頓-拉夫遜潮流算法的基本原理，利用 C+編程計(jì)算了一個(gè)4節(jié) 點(diǎn)簡(jiǎn)單電力網(wǎng)絡(luò)的潮流，并驗(yàn)證了運(yùn)行成果，誤差在允許范圍之內(nèi)。因?yàn)榕＠?法計(jì)算過(guò)程中要不斷

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

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

21、ance ;double *admittanceAmplitude;double *admittanceAngle;double *lineData;int PVNodeNumber; /PV 節(jié)點(diǎn)數(shù)double *jacobiMatrix; / 雅克比矩陣double *PQData; / 矩陣，存放 PQ 節(jié)點(diǎn)的數(shù)據(jù)double *PVData; / 矩陣，存放 PV 節(jié)點(diǎn)的數(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;/線路的計(jì)算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);/變壓器的計(jì)算，規(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); /對(duì)地支路的計(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é)點(diǎn)數(shù)開辟數(shù)組void initi() nodeNumber=readParameter0; / 讀取節(jié)點(diǎn)數(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;/計(jì)算導(dǎo)納矩陣主程序void countAdmittance()int lineDataSerial=1;for(int facility=0;facilityreadParameter0;if (readParameter00.5) break; / 序號(hào)為 0 時(shí)退出循環(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+;/ 不同元件的計(jì)算/end for if(facility=2) /end for loop iswitch(facility) case 0:initi();balanceNodeindex=(int)readParameter1;balanceNodeVoltage=readParameter2; calculationAccuracy=readParameter3; cout 節(jié)點(diǎn)數(shù) :nodeNumberendl 平衡節(jié)點(diǎn) :readParameter1endl 平衡節(jié)點(diǎn)電壓 :readParameter2endl 計(jì)算精度 :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é)點(diǎn)和PV節(jié)點(diǎn)的數(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é)點(diǎn)的值已經(jīng)在初始化中賦值了，這里給平衡節(jié)點(diǎn)和 PV 節(jié)點(diǎn)賦值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é)點(diǎn)號(hào)sum=0.0