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

0 引言潮流是配電網絡分析的基礎，用于電網調度、運行分析、操作模擬和設計規(guī)劃，同時也是電壓優(yōu)化和網絡接線變化所要參考的內容。潮流計算通過數值仿真的方法把電力系統的詳細運行情況呈現給工作人員，從而便于研究系統在給定條件下的穩(wěn)態(tài)運行特點。隨著市場經濟的發(fā)展，經濟利益是企業(yè)十分看重的，而線損卻是現階段阻礙企業(yè)提高效益的一大因素。及時、準確的潮流計算結果，可以給出配電網的潮流分布、理論線損及其在網絡中的分布，從而為配電網的安全經濟運行提供參考。從數學的角度來看，牛頓-拉夫遜法能有效進行非線性代數方程組的計算且具有二次收斂的特點，具有收斂快、精度高的特點，在輸電網中得到廣泛應用。隨著現代計算機技術的發(fā)展，利用編程和相關軟件，可以更好、更快地實現配電網功能，本文就是結合牛頓-拉夫遜法的基本原理，利用C+程序進行潮流計算，計算結果表明該方法具有良好的收斂性、可靠性及正確性。1 牛頓-拉夫遜法基本介紹1.1 潮流方程對于N個節(jié)點的電力網絡（地作為參考節(jié)點不包括在內），如果網絡結構和元件參數已知，則網絡方程可表示為： （1-1）式中，為N*N階節(jié)點導納矩陣；為N*1維節(jié)點電壓列向量；為N*1維節(jié)點注入電流列向量。如果不計網絡元件的非線性，也不考慮移相變壓器，則為對稱矩陣。電力系統計算中，給定的運行變量是節(jié)點注入功率，而不是節(jié)點注入電流，這兩者之間有如下關系： （1-2）式中，為節(jié)點的注入復功率，是N*1維列矢量；為的共軛；是由節(jié)點電壓的共軛組成的N*N階對角線矩陣。由（1-1）和（1-2），可得：上式就是潮流方程的復數形式，是N維的非線性復數代數方程組。將其展開，有： j=1,2，.，N （1-3）式中， 表示所有和相連的節(jié)點，包括。將節(jié)點電壓用極坐標表示，即令，代入式（1-3）中則有：故有： i=1,2,N （1-4）式（1-4）是用極坐標表示的潮流方程。而節(jié)點功率誤差： (1-5) (1-6)式中：，為節(jié)點給定的有功功率及無功功率。 1.2 牛頓-拉夫遜法基本原理1.2.1 牛拉法的一般描述牛拉法是把非線性方程式的求解過程變成反復對相應的線性方程式的求解過程，即非線性問題通過線性化逐步近似，這就是牛拉法的核心。下面以非線性方程式的求解過程來進行說明。設電力網絡的節(jié)點功率方程一般形式如下： （1-7）式中，為節(jié)點注入功率給定值；為對應的物理量和節(jié)點電壓之間的函數表達式；為節(jié)點電壓。寫成功率偏差的形式： （1-8）應用牛拉法求解如下。在給定的初值處將式（1-8）作一階泰勒展開：定義為潮流方程的雅克比矩陣，為在處的值，則有：用修正就得到的新值。如果用k表示迭代次數，寫成一般的表達式，有： （1-9）對于潮流收斂的情況，應比更接近于解點。收斂條件為： 由簡單迭代法收斂性分析的結論知，越接近解點，牛頓-拉夫遜法收斂越快，它具有二階收斂速度。由圖1.1可以直觀地了解牛拉法的步驟：圖1.1 牛頓-拉夫遜法的幾何解釋1.2.2 極坐標的牛頓-拉夫遜法在極坐標中，有如下的形式： （1-10）共2n-r個方程，狀態(tài)變量為：共2n-r個待求量。r 個PV節(jié)點的電壓幅值給定，不需求解。潮流雅克比矩陣的維數是（2n-r）*(2n-r)，結構如下：上式右側的對電壓幅值的偏導數項中的電壓幅值的階數減少了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中的“計算電壓幅值和角度”步驟較多，單獨用圖2.2表示出來。圖2.1 牛頓法計算潮流的程序框圖 圖2.2 電壓幅值和角度求解步驟框圖當不符合收斂的條件“amontk1”時,即認為計算不收斂。具體程序見附錄。2.2 計算步驟下面討論的是極坐標形式的牛頓法求解過程，大致分為以下幾個步驟： 形成節(jié)點導納矩陣； 給各節(jié)點電壓設初值（）；根據式（1-12）、（1-13）生成雅克比矩陣（H、N、M、L）； 將節(jié)點電壓初值代入式(1-5)、式(1-6)，求出修正方程式的常數項向量； 求解修正方程，得到電壓幅值和角度；判斷是否收斂，若收斂，計算平衡節(jié)點和線路功率； 輸出結果，并結束。3 算例3.1 系統模型本文以圖3.1所示電力網絡為例，調用基于牛頓-拉夫遜法的C+程序。圖3.1 系統模型其中節(jié)點4設為平衡節(jié)點，電壓標幺值為1.05，計算誤差為0.。3.2 輸入與輸出將圖3.1所示模型的相關數據放在data.dat文件中圖3.2 輸入節(jié)點和支路數據對各個數字含義的解釋如下：網絡模型有四個節(jié)點，四條支路，編號見圖3.1。第一個零下面三行數為支路參數，分別表示三條支路的起始和終止節(jié)點編號，后面的為電阻、電抗和電納，電導均為0，例如：1 2 0.1 0.4 0.01528。第二個零下面的為變壓器支路，各數字意義同支路參數。接下去三行均為節(jié)點參數，分別表示注入有功功率和無功功率。調用text.cpp文件，得到運行結果，見圖3.3和圖3.4。圖3.3 運行結果1圖3.4 運行結果23.3 結果分析將上述仿真結果整理為表格3.1、3.2，其中“+”表示節(jié)點i輸出功率給節(jié)點j，“-”表示節(jié)點j輸出功率給i(縱向為i，橫向為j)。表3.1節(jié)點有功功率輸入與輸出節(jié)點號123410+0.-0.5-0.2-0.2431600-0.30.500040.0.00表3.2節(jié)點無功功率輸入與輸出節(jié)點號123410-0.-0.-0.20.00-0.1403630.00040.104640.00根據表格計算：節(jié)點1有功功率：0+0.-0.5-0.=-0. 無功功率：0-0.-0.-0.=-0.節(jié)點2有功功率：-0.24316+0+0-0.=-0.無功功率：0.+0+0-0.14036=-0.節(jié)點3有功功率：0.5+0+0+0=0.5無功功率：0.+0+0+0=0.節(jié)點4有功功率：0.+0.+0+0=0.無功功率：0.10464+0.+0+0=0.根據已知條件，兩個PQ節(jié)點的注入有功、無功分別為：P1=0.3，Q1=0.18；P2=0.55，Q2=0.13潮流計算誤差：可見，誤差均在允許范圍內。線路損耗： 3.4 結論通過上面的分析與計算，驗證了程序的正確性。由于編寫過程的不足，線路損耗沒能直接計算出來，而是需要手算，比較遺憾。程序在運行過程中，需要區(qū)分三種不同的節(jié)點，這由子程序保證實現。相比于快速分解法，牛拉法程序較為復雜，但更精確一點，潮流誤差較小。4 總結本文基于牛頓-拉夫遜潮流算法的基本原理，利用C+編程計算了一個4節(jié)點簡單電力網絡的潮流，并驗證了運行成果，誤差在允許范圍之內。因為牛拉法計算過程中要不斷生成新的雅各布矩陣，所以相對來說占用內存較多，但收斂速度快，這在程序運行過程中可以體現出來。本文程序并不是特別實用，因為真正的電力網絡不可能只有幾個節(jié)點，而且各種電力設備的情況也會復雜很多，因此程序會變得非常大，占用極大內存。但是，我還是通過這次練習，進一步鞏固了書本上的理論知識，了解了實際操作的過程步驟。最后，感謝楊偉老師的悉心指導！參考文獻1 朱紅，趙琦，王慶寶. C+程序設計教程M. 北京：清華大學出版社，2010.2張伯明，陳壽孫，嚴正.高等電力網絡分析M.北京：清華大學出版社，2007.3吳明波. 牛頓-拉夫遜法在潮流計算中的應用J. 內蒙古科技與經濟，2011, 21：111-112，115. 4明日科技. VisualC+從入門到精通M. 北京：清華大學出版社，2012.5 顧潔，陳章潮，徐蓓. 一種新的配電網潮流算法-改進牛頓拉夫遜法J.華東電力，2000，5:10-12.6 朱文強. 牛頓-拉夫遜法在配電網中的應用J. 水利科技，2004，3:55-56,58.7 劉明波，謝敏，趙維新. 大電網最優(yōu)潮流計算M. 北京：科學出版社，2010.31附錄/*/#includeiostream.h#includemath.h#include#include#include#define nodeNumber readParameter6using namespace std;/class powerFlowCalculation /導納矩陣的計算類private: int readDataAmount; /存放元件參數的讀取個數int balanceNodeindex; /平衡結點號double balanceNodeVoltage; /平衡節(jié)點電壓幅值double balanceNodePowerP; /平衡節(jié)點功率double balanceNodePowerQ;double calculationAccuracy; /計算精度double readParameter7; /宏定義double *conductance;double *susceptance ; double *admittanceAmplitude;double *admittanceAngle;double *lineData;int PVNodeNumber; /PV節(jié)點數double *jacobiMatrix; /雅克比矩陣double *PQData; /矩陣，存放PQ節(jié)點的數據 double *PVData; /矩陣，存放PV節(jié)點的數據double *voltageAmplitude; /電壓幅值double *voltageAngle; /電壓角度double *constantVector; /常數向量double *lineConsumeG;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); 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,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;/導納矩陣的顯示void display() coutendl導納矩陣為：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;/構造函數,打開文件powerFlowCalculation(string inputfilename) instream.open(inputfilename.data();assert(instream.is_open();/析構函數，關閉文件 powerFlowCalculation() instream.close();/初始化，根據節(jié)點數開辟數組void initi() nodeNumber=readParameter0; /讀取節(jié)點數 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+) /導納初始化 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+) /常數向量初始化 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(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é)點數:nodeNumberendl平衡節(jié)點:readParameter1endl平衡節(jié)點電壓: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;/將直角坐標的導納矩陣轉換為極坐標形式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é)點和PV節(jié)點的數據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é)點的值已經在初始化中賦值了，這里給平衡節(jié)點和PV節(jié)點賦值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的下標轉換int converSerial(int i,int j,int rank=6) int serial;serial=(i-1)*rank+j;return serial;/生成常數向量void creatConstantVector() int i; /i為P、Q的下標double sum=0.0;double beta=0.0;for(int n=1;n=(nodeNumber-1)*2;n+) /n為向量的下標i=(n+1)/2; /i為節(jié)點號 sum=0.0;if(n+1)%2=1) /求Q點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