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

0 引言潮流是配電網(wǎng)絡(luò)分析的基礎(chǔ)，用于電網(wǎng)調(diào)度、運(yùn)行分析、操作模擬和設(shè)計規(guī)劃，同時也是電壓優(yōu)化和網(wǎng)絡(luò)接線變化所要參考的內(nèi)容。潮流計算通過數(shù)值仿真的方法把電力系統(tǒng)的詳細(xì)運(yùn)行情況呈現(xiàn)給工作人員，從而便于研究系統(tǒng)在給定條件下的穩(wěn)態(tài)運(yùn)行特點(diǎn)。隨著市場經(jīng)濟(jì)的發(fā)展，經(jīng)濟(jì)利益是企業(yè)十分看重的，而線損卻是現(xiàn)階段阻礙企業(yè)提高效益的一大因素。及時、準(zhǔn)確的潮流計算結(jié)果，可以給出配電網(wǎng)的潮流分布、理論線損及其在網(wǎng)絡(luò)中的分布，從而為配電網(wǎng)的安全經(jīng)濟(jì)運(yùn)行提供參考。從數(shù)學(xué)的角度來看，牛頓-拉夫遜法能有效進(jìn)行非線性代數(shù)方程組的計算且具有二次收斂的特點(diǎn)，具有收斂快、精度高的特點(diǎn)，在輸電網(wǎng)中得到廣泛應(yīng)用。隨著現(xiàn)代計算機(jī)技術(shù)的發(fā)展，利用編程和相關(guān)軟件，可以更好、更快地實現(xiàn)配電網(wǎng)功能，本文就是結(jié)合牛頓-拉夫遜法的基本原理，利用C+程序進(jìn)行潮流計算，計算結(jié)果表明該方法具有良好的收斂性、可靠性及正確性。1 牛頓-拉夫遜法基本介紹1.1 潮流方程對于N個節(jié)點(diǎn)的電力網(wǎng)絡(luò)（地作為參考節(jié)點(diǎn)不包括在內(nèi)），如果網(wǎng)絡(luò)結(jié)構(gòu)和元件參數(shù)已知，則網(wǎng)絡(luò)方程可表示為： （1-1）式中，為N*N階節(jié)點(diǎn)導(dǎo)納矩陣；為N*1維節(jié)點(diǎn)電壓列向量；為N*1維節(jié)點(diǎn)注入電流列向量。如果不計網(wǎng)絡(luò)元件的非線性，也不考慮移相變壓器，則為對稱矩陣。電力系統(tǒng)計算中，給定的運(yùn)行變量是節(jié)點(diǎn)注入功率，而不是節(jié)點(diǎn)注入電流，這兩者之間有如下關(guān)系： （1-2）式中，為節(jié)點(diǎn)的注入復(fù)功率，是N*1維列矢量；為的共軛；是由節(jié)點(diǎn)電壓的共軛組成的N*N階對角線矩陣。由（1-1）和（1-2），可得：上式就是潮流方程的復(fù)數(shù)形式，是N維的非線性復(fù)數(shù)代數(shù)方程組。將其展開，有： j=1,2，.，N （1-3）式中， 表示所有和相連的節(jié)點(diǎn)，包括。將節(jié)點(diǎn)電壓用極坐標(biāo)表示，即令，代入式（1-3）中則有：故有： i=1,2,N （1-4）式（1-4）是用極坐標(biāo)表示的潮流方程。而節(jié)點(diǎn)功率誤差： (1-5) (1-6)式中：，為節(jié)點(diǎn)給定的有功功率及無功功率。 1.2 牛頓-拉夫遜法基本原理1.2.1 牛拉法的一般描述牛拉法是把非線性方程式的求解過程變成反復(fù)對相應(yīng)的線性方程式的求解過程，即非線性問題通過線性化逐步近似，這就是牛拉法的核心。下面以非線性方程式的求解過程來進(jìn)行說明。設(shè)電力網(wǎng)絡(luò)的節(jié)點(diǎn)功率方程一般形式如下： （1-7）式中，為節(jié)點(diǎn)注入功率給定值；為對應(yīng)的物理量和節(jié)點(diǎn)電壓之間的函數(shù)表達(dá)式；為節(jié)點(diǎn)電壓。寫成功率偏差的形式： （1-8）應(yīng)用牛拉法求解如下。在給定的初值處將式（1-8）作一階泰勒展開：定義為潮流方程的雅克比矩陣，為在處的值，則有：用修正就得到的新值。如果用k表示迭代次數(shù)，寫成一般的表達(dá)式，有： （1-9）對于潮流收斂的情況，應(yīng)比更接近于解點(diǎn)。收斂條件為： 由簡單迭代法收斂性分析的結(jié)論知，越接近解點(diǎn)，牛頓-拉夫遜法收斂越快，它具有二階收斂速度。由圖1.1可以直觀地了解牛拉法的步驟：圖1.1 牛頓-拉夫遜法的幾何解釋1.2.2 極坐標(biāo)的牛頓-拉夫遜法在極坐標(biāo)中，有如下的形式： （1-10）共2n-r個方程，狀態(tài)變量為：共2n-r個待求量。r 個PV節(jié)點(diǎn)的電壓幅值給定，不需求解。潮流雅克比矩陣的維數(shù)是（2n-r）*(2n-r)，結(jié)構(gòu)如下：上式右側(cè)的對電壓幅值的偏導(dǎo)數(shù)項中的電壓幅值的階數(shù)減少了1，為使雅克比矩陣的各部分子矩陣具有一致的形式，在實際計算中，常將該項乘以電壓幅值，并選取作為待求的修正量，則雅克比矩陣可寫成： （1-11）將式（1-10）和（1-11）代入式（1-9）的修正方程即可求得x的修正量，用它修正x直到為止。將式（1-11）用下式表示：其中每個字塊的計算公式如下：對角元素： （1-12）非對角元素： （1-13）2 牛頓法潮流計算步驟2.1 程序流程圖在了解了牛拉法的原理之后，明確程序編寫思路，如圖2.1、2.2所示。 其中圖2.1中的“計算電壓幅值和角度”步驟較多，單獨(dú)用圖2.2表示出來。圖2.1 牛頓法計算潮流的程序框圖 圖2.2 電壓幅值和角度求解步驟框圖當(dāng)不符合收斂的條件“amontk1”時,即認(rèn)為計算不收斂。具體程序見附錄。2.2 計算步驟下面討論的是極坐標(biāo)形式的牛頓法求解過程，大致分為以下幾個步驟： 形成節(jié)點(diǎn)導(dǎo)納矩陣； 給各節(jié)點(diǎn)電壓設(shè)初值（）；根據(jù)式（1-12）、（1-13）生成雅克比矩陣（H、N、M、L）； 將節(jié)點(diǎn)電壓初值代入式(1-5)、式(1-6)，求出修正方程式的常數(shù)項向量； 求解修正方程，得到電壓幅值和角度；判斷是否收斂，若收斂，計算平衡節(jié)點(diǎn)和線路功率； 輸出結(jié)果，并結(jié)束。3 算例3.1 系統(tǒng)模型本文以圖3.1所示電力網(wǎng)絡(luò)為例，調(diào)用基于牛頓-拉夫遜法的C+程序。圖3.1 系統(tǒng)模型其中節(jié)點(diǎn)4設(shè)為平衡節(jié)點(diǎn)，電壓標(biāo)幺值為1.05，計算誤差為0.。3.2 輸入與輸出將圖3.1所示模型的相關(guān)數(shù)據(jù)放在data.dat文件中圖3.2 輸入節(jié)點(diǎn)和支路數(shù)據(jù)對各個數(shù)字含義的解釋如下：網(wǎng)絡(luò)模型有四個節(jié)點(diǎn)，四條支路，編號見圖3.1。第一個零下面三行數(shù)為支路參數(shù)，分別表示三條支路的起始和終止節(jié)點(diǎn)編號，后面的為電阻、電抗和電納，電導(dǎo)均為0，例如：1 2 0.1 0.4 0.01528。第二個零下面的為變壓器支路，各數(shù)字意義同支路參數(shù)。接下去三行均為節(jié)點(diǎn)參數(shù)，分別表示注入有功功率和無功功率。調(diào)用text.cpp文件，得到運(yùn)行結(jié)果，見圖3.3和圖3.4。圖3.3 運(yùn)行結(jié)果1圖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(縱向為i，橫向為j)。表3.1節(jié)點(diǎn)有功功率輸入與輸出節(jié)點(diǎn)號123410+0.-0.5-0.2-0.2431600-0.30.500040.0.00表3.2節(jié)點(diǎn)無功功率輸入與輸出節(jié)點(diǎn)號123410-0.-0.-0.20.00-0.1403630.00040.104640.00根據(jù)表格計算：節(jié)點(diǎn)1有功功率：0+0.-0.5-0.=-0. 無功功率：0-0.-0.-0.=-0.節(jié)點(diǎn)2有功功率：-0.24316+0+0-0.=-0.無功功率：0.+0+0-0.14036=-0.節(jié)點(diǎn)3有功功率：0.5+0+0+0=0.5無功功率：0.+0+0+0=0.節(jié)點(diǎn)4有功功率：0.+0.+0+0=0.無功功率：0.10464+0.+0+0=0.根據(jù)已知條件，兩個PQ節(jié)點(diǎn)的注入有功、無功分別為：P1=0.3，Q1=0.18；P2=0.55，Q2=0.13潮流計算誤差：可見，誤差均在允許范圍內(nèi)。線路損耗： 3.4 結(jié)論通過上面的分析與計算，驗證了程序的正確性。由于編寫過程的不足，線路損耗沒能直接計算出來，而是需要手算，比較遺憾。程序在運(yùn)行過程中，需要區(qū)分三種不同的節(jié)點(diǎn)，這由子程序保證實現(xiàn)。相比于快速分解法，牛拉法程序較為復(fù)雜，但更精確一點(diǎn)，潮流誤差較小。4 總結(jié)本文基于牛頓-拉夫遜潮流算法的基本原理，利用C+編程計算了一個4節(jié)點(diǎn)簡單電力網(wǎng)絡(luò)的潮流，并驗證了運(yùn)行成果，誤差在允許范圍之內(nèi)。因為牛拉法計算過程中要不斷生成新的雅各布矩陣，所以相對來說占用內(nèi)存較多，但收斂速度快，這在程序運(yùn)行過程中可以體現(xiàn)出來。本文程序并不是特別實用，因為真正的電力網(wǎng)絡(luò)不可能只有幾個節(jié)點(diǎn)，而且各種電力設(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-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.31附錄/*/#includeiostream.h#includemath.h#include#include#include#define nodeNumber readParameter6using namespace std;/class powerFlowCalculation /導(dǎo)納矩陣的計算類private: int readDataAmount; /存放元件參數(shù)的讀取個數(shù)int balanceNodeindex; /平衡結(jié)點(diǎn)號double balanceNodeVoltage; /平衡節(jié)點(diǎn)電壓幅值double balanceNodePowerP; /平衡節(jié)點(diǎn)功率double balanceNodePowerQ;double calculationAccuracy; /計算精度double readParameter7; /宏定義double *conductance;double *susceptance ; 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 *lineConsumeG;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); 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,double*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=-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)納矩陣為：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 ; else coutconductancenij+susceptance nijj ; coutendl;coutendl;/顯示任意秩矩陣void displayMatrix(double *matr,int ranki,int rankj) int nij; for(int i=1;i=ranki;i+) for(int j=1;j=rankj;j+) nij=(i-1)*rankj+j; coutmatrnij ; coutendl;/coutendl;/構(gòu)造函數(shù),打開文件powerFlowCalculation(string inputfilename) instream.open(inputfilename.data();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; 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; /開辟PVDatavoltageAmplitude=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=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(int 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(facility=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) lineDataSerial+; /end for if(facility=2) /end for loop i switch(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計算精度:readParameter3endl; 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 vectorAngle(double a,double b)if(fabs(a)0. & fabs(b)0.) return 0.0; double alpha; alpha=atan(fabs(b)/fabs(a);if(a=0)alpha=3.-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;j=nodeNumber;j+) n=converIndex(i,j);admittanceAmplituden=sqrt(conductancen*conductancen+susceptancen*susceptancen);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 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=balanceNodeVoltage;voltageAnglebalanceNodeindex=0.0;for(int i=1;i=PVNodeNumber;i+)n=PVDataconverIndex(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(i,k)*voltageAmplitudek*sin(beta);/end for loop ksum*=voltageAmplitudei;sum=voltageAmplitudei*voltageAmplitudei*admittanceAmplitudeconverIndex(i,i)*sin(alpha);sum=-sum;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*sin(beta); return sum;/生成N矩陣double creatN(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(i,k)*voltageAmplitudek*cos(beta);/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 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=-sum;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*cos(beta);/end for ifreturn sum;/生成H矩陣double creatH(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*sin(beta);sum*=voltageAmplitudei;elsebeta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum=-1*voltageAmplitudei*admittanceAmplitudeconverIndex(i,j)*voltageAmplitudej*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+nodeNumber*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*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 nserial=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)號 sum=0.0;if(n+1)%2=1) /求Q點(diǎn)for(int j=1;j=nodeNumber;j+)beta=voltageAnglei-voltageAnglej-admittanceAngleconverIndex(i,j);sum+=voltageAmplitudei*voltageAmplitudej*admittanceAmplitudeconverIndex(i,j)*sin(beta);/end for loop jif(n=(nodeNumber-PVNodeNumber-1)*2)constantVectorn=-(PQDataconverSerial(i,4