牛拉法潮流計(jì)算程序(附3機(jī)9節(jié)點(diǎn)結(jié)果對(duì)比)(共11頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上 摘要電力系統(tǒng)潮流計(jì)算是研究電力系統(tǒng)穩(wěn)態(tài)運(yùn)行的一種重要方法,它根據(jù)給定的運(yùn)行條件及系統(tǒng)接線情況確定整個(gè)電力系統(tǒng)各部分的運(yùn)行狀態(tài),包括各母線的電壓、線路的功率分布以及功率損耗等等。潮流計(jì)算主要用于電網(wǎng)規(guī)劃和靜態(tài)安全分析，它可為擴(kuò)建電力網(wǎng)絡(luò)，以達(dá)到規(guī)劃周期內(nèi)所需要的輸電能力提供依據(jù)；也可以對(duì)預(yù)想事故進(jìn)行模擬和分析，校核預(yù)想事故下的電力系統(tǒng)安全性。本文簡(jiǎn)單介紹了牛頓-拉夫遜潮流計(jì)算的原理、模型與算法，然后用具體的實(shí)例，利用MATLAB對(duì)牛頓-拉夫遜法的算法進(jìn)行了驗(yàn)證。關(guān)鍵詞：電力系統(tǒng)潮流計(jì)算 牛頓-拉夫遜法 MATLAB一、牛拉法的數(shù)學(xué)模型對(duì)一個(gè)N節(jié)點(diǎn)的電力網(wǎng)路，列寫節(jié)點(diǎn)電

2、壓方程，即I=Y V(1.1)式中，I為節(jié)點(diǎn)注入電流列相量，Y為節(jié)點(diǎn)導(dǎo)納矩陣，V為節(jié)點(diǎn)電壓列相量。由于異地測(cè)量的兩個(gè)電流缺少時(shí)間同步信息，以注入功率替換注入電流作為已知量。即(1.2)其中,Yij=Gij+jBij，帶入上式，得到有功功率和無(wú)功功率方程Pi=Vij=1nVjGijcosij+Bijsinij(1.3)Qi=Vij=1nVjGijsinij-Bijcosij(1.4)大部分情況下，已知PQ，求解V。考慮到電網(wǎng)的功率平衡，至少選擇一臺(tái)發(fā)電機(jī)來(lái)平衡全網(wǎng)有功功率，即至少有一個(gè)平衡節(jié)點(diǎn)，常選擇調(diào)頻或出線較多的發(fā)電機(jī)作為平衡節(jié)點(diǎn)。具有無(wú)功補(bǔ)償?shù)哪妇€能保持電壓幅值恒定，這類節(jié)點(diǎn)可作為PV節(jié)點(diǎn)

3、。 潮流計(jì)算中節(jié)點(diǎn)分類總結(jié)如下：表 1.1 節(jié)點(diǎn)分類PQ節(jié)點(diǎn)負(fù)荷母線、變電站母線絕大部分是PQ節(jié)點(diǎn)已知PQ，求VPV節(jié)點(diǎn)有無(wú)功儲(chǔ)備的發(fā)電機(jī)母線很小部分是PV節(jié)點(diǎn)已知PV（V=1），求Q平衡節(jié)點(diǎn)選擇一個(gè)發(fā)電機(jī)母線已知V（V=1,=0）,求PQ電壓相量作為參考軸=0已知電力系統(tǒng)有m個(gè)PQ節(jié)點(diǎn)，r個(gè)PV節(jié)點(diǎn)和1個(gè)平衡節(jié)點(diǎn)，則可以提取m+r個(gè)有功功率方程和m個(gè)無(wú)功功率方程，從而求解出m+r個(gè)和m個(gè)V，其余節(jié)點(diǎn)的有功和無(wú)功可通過(guò)式（1.3）、（1.4）求得，這樣就完成了潮流計(jì)算。二、潮流計(jì)算流程用于潮流計(jì)算的有功和無(wú)功功率方程的階數(shù)一般很高，常選擇牛頓-拉夫遜迭代法進(jìn)行求解高維的非線性方程組。結(jié)束max

4、F(Xk)<？XkXk+1=Xk+XkmaxXk<？Xk=-J-1XkFXk牛頓-拉夫遜法的迭代過(guò)程如下：圖2.1 牛拉法迭代流程圖將功率方程改成Fx=0的形式Pi=Pis-Vij=1nVjGijcosij+Bijsinij=0(2.1)Qi=Qis-Vij=1nVjGijsinij-Bijcosij=0(2.2)可以得到j(luò)Vj=-PijPiVjQijQiVj-1PiQi(2.3)為了便于計(jì)算，將上式改為jVjVj=-PijPiVjVjQijQiVjVj-1PiQi(2.4)求取功率方程的雅克比矩陣ij時(shí)，Pij=-ViVjGijsinij-Bijcosij(2.5)PiVjVj=

5、-ViVjGijcosij+Bijsinij(2.6)Qij=ViVjGijcosij+Bijsinij(2.7)QiVjVj=-ViVjGijsinij-Bijcosij(2.8)i=j時(shí)，Pii=Vi2Bii+Qi(2.9)PiViVi=-Vi2Gii-Pi(2.10)Qii=Vi2Gii-Pi(2.11)QiViVi=Vi2Bii-Qi(2.12)牛拉法求解潮流的迭代過(guò)程如下圖2.2 牛拉法潮流迭代流程圖三、MATLAB編程實(shí)例實(shí)例采用王錫凡主編的現(xiàn)代電力系統(tǒng)分析p326-327的例題。圖3.1 某小型電力網(wǎng)絡(luò)圖例題中，1號(hào)發(fā)電機(jī)為平衡節(jié)點(diǎn)，2號(hào)和3號(hào)發(fā)電機(jī)為PV節(jié)點(diǎn)，其余為PQ節(jié)點(diǎn)。為

6、了便于編程，可將PV節(jié)點(diǎn)和平衡節(jié)點(diǎn)放在節(jié)點(diǎn)描述圖的最后，故要重新對(duì)節(jié)點(diǎn)進(jìn)行編號(hào)。圖3.2 修改節(jié)點(diǎn)編號(hào)后的電力網(wǎng)絡(luò)圖支路表如下：%首節(jié)點(diǎn)末節(jié)點(diǎn)電阻電抗 容納之半變比120.0100.0850.0881130.0170.0920.0791240.0320.1610.1531360.0390.1700.1791450.00850.0720.07451560.01190.10080.104519100.0576017400.0625018600.058601節(jié)點(diǎn)表如下：%節(jié)點(diǎn)電壓幅值 電壓相角 節(jié)點(diǎn)有功節(jié)點(diǎn)無(wú)功節(jié)點(diǎn)類型 /“1”為PQ節(jié)點(diǎn)，“2”為PV節(jié)點(diǎn)，“3”為平衡節(jié)點(diǎn)11.0 000121.0

7、 0-1.2500-0.5131.0 0-0.9-0.3141.0 000151.0 0-1-0.35161.0 000171.025 01.63000281.025 00.85000291.040 0003程序清單：取=10-5a=load('zhilu.txt');b=load('jiedian.txt');N=size(a,1);%支路數(shù)Nbus=size(b,1);%節(jié)點(diǎn)數(shù)pq=0;for k=1:Nbus %PQ節(jié)點(diǎn)個(gè)數(shù) if b(k,6)=1 pq=pq+1; endendY=zeros(Nbus);%節(jié)點(diǎn)導(dǎo)納矩陣for k=1:N t1=a(k,1

8、);t2=a(k,2);r=a(k,3);x=a(k,4);ban=a(k,5);K=a(k,6); Y(t1,t1)=Y(t1,t1)+1/(r+j*x)+j*ban; Y(t1,t2)=Y(t1,t2)-1/(K*(r+j*x); Y(t2,t1)=Y(t2,t1)-1/(K*(r+j*x); Y(t2,t2)=Y(t2,t2)+1/(K*K*(r+j*x)+j*ban;endG=real(Y);B=imag(Y); precision=1;t=0；%存儲(chǔ)迭代次數(shù)while precision>0.00001 P=zeros(Nbus,1);%存儲(chǔ)所有節(jié)點(diǎn)的有功 Q=zeros(Nb

9、us,1);%存儲(chǔ)所有節(jié)點(diǎn)的無(wú)功 for m=1:Nbus %求pv、pq和平衡節(jié)點(diǎn)的有功 for n=1:Nbus P(m,1)=P(m,1)+b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); end end for m=1:Nbus %求pq、pv和平衡節(jié)點(diǎn)的無(wú)功 for n=1:Nbus Q(m,1)=Q(m,1)+b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); end end deltp=b(1:Nbus-1,4)-P(1

10、:Nbus-1,1);%pq和pv節(jié)點(diǎn)有功差 deltq=b(1:pq,5)-Q(1:pq,1);%pq節(jié)點(diǎn)無(wú)功差 deltPQ=deltp;deltq; H=zeros(Nbus-1); for m=1:Nbus-1 %求H矩陣 for n=1:Nbus-1 if m=n H(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); else H(m,m)=b(m,2)*b(m,2)*B(m,m)+Q(m,1); end end end N=zeros(Nbus-1,pq); for m=1:Nbus-1

11、%求N矩陣 for n=1:pq if m=n N(m,n)=-b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); else N(m,m)=-b(m,2)*b(m,2)*G(m,m)-P(m,1); end end end J=zeros(pq,Nbus-1); for m=1:pq %求J矩陣 for n=1:Nbus-1 if m=n J(m,n)=b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3)+B(m,n)*sin(b(m,3)-b(n,3); else J(m,m)=b(m,

12、2)*b(m,2)*G(m,m)-P(m,1); end end end L=zeros(pq,pq); for m=1:pq %求L矩陣 for n=1:pq if m=n L(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3)-B(m,n)*cos(b(m,3)-b(n,3); else L(m,m)=b(m,2)*b(m,2)*B(m,m)-Q(m,1); end end end Jacobi=H N;J L;%雅克比矩陣 Correction=-JacobideltPQ;%計(jì)算電壓相角和幅值的修正量 for m=1:Nbus-1 b(m,3)=b(m,3)+Correction(m); end for m=1:pq b(m,2)=b(m,2)+Correction(Nbus-1+m)*b(m,2); endprecision=max(abs(deltPQ); t=t+1 bendb(Nbus,4)=P(Nbus,1);b(Nbus,