三機九節(jié)點電力系統(tǒng)仿真matlab

1、電力系統(tǒng)仿真作業(yè)三機九節(jié)點電力系統(tǒng)暫態(tài)仿真學(xué)院：能源與動力工程學(xué)院專業(yè)：電力系統(tǒng)及其自動化學(xué)號：姓名：于永生導(dǎo)師：授課教師：目錄一、 概述.1二、 課程主要任務(wù).11.系統(tǒng)數(shù)據(jù).12.潮流計算.23.負荷等效和支路簡化.44.求解電磁功率.55.求解運動方程.56.程序活單.7(1) .主程序：.7(2) .極坐標轉(zhuǎn)換成直角坐標函數(shù)pol2rect(V ,del).16(3) .直角坐標轉(zhuǎn)換成極坐標函數(shù)rect2pol(Z) . 16(4) .求解微分方程所用的得到微分量的函數(shù)Gen_fw(t,X,Y_Gen,E,Pm0,Tj) . 16三、 課程總結(jié)及心得體會 .16四、 參考文獻.17于永

2、生電力系統(tǒng)仿真作業(yè)1一、概述在動態(tài)穩(wěn)定分析中，系統(tǒng)由線性化的微分方程組和代數(shù)方程組描寫，并用經(jīng)典的或現(xiàn) 代的線性系統(tǒng)理論來進行穩(wěn)定分析，分析可以在時域或頻域進行。當用計算機和現(xiàn)代線性 系統(tǒng)理論分析時，常把系統(tǒng)線性化的微分方程組和代數(shù)方程組消去代數(shù)變量，化為狀態(tài)方 程形式，并廣泛采用特征分析進行穩(wěn)定分析。電力系統(tǒng)是由不同類型的發(fā)電機組、多種電力負荷、不同電壓等級的電力網(wǎng)絡(luò)等組成 的十分龐大復(fù)雜的動力學(xué)系統(tǒng)。其暫態(tài)過渡過程不僅包括電磁方面的過渡過程，而且還有 機電方面的過渡過程。由此可見，電力系統(tǒng)的數(shù)學(xué)模型是一個強非線性的高維狀態(tài)方程組。 在動態(tài)穩(wěn)定仿真中使用簡單的電力系統(tǒng)模型，發(fā)電機用三階模型表

3、示。二、課程主要任務(wù)P. M. Anderson and A. A. Fouad編寫的Power System Control and WesternSystem Coordinated Council (WSCC)三機九節(jié)點系統(tǒng)模型。1.系統(tǒng)數(shù)據(jù)其中，節(jié)點數(shù)據(jù)如下:%節(jié)點數(shù)據(jù)% %節(jié)點電壓電壓發(fā)電機發(fā)電機負荷負荷節(jié)點%號幅值相角有功無功有功無功類型(1PQ 2PV 3平衡）衡）N=11.0400.71640.270500321.025 01.630.066500231.025 00.85-0.108600241000001510001.250.51610000.90.317100000181

4、00010.35，191000001;其中，支路數(shù)據(jù)如下:%線路數(shù)據(jù)%首端末端電阻電抗電納(1/2)變壓器非標準變比L=450.010.0850.0881460.0170.0920.0791570.0320.1610.1531690.0390.170.1791780.00850.0720.07451本次課程主要應(yīng)用Stability 一書中所引用的電力系統(tǒng)仿真作業(yè)于永生2890.01190.10080.105411400.0576012700.0625013900.058601;發(fā)電機數(shù)據(jù)如下：%發(fā)電機母線XdXdTd0XqXqTq0TjXfGe= 110.14600.06088.960.09

5、690.0969047.280.0576220.89580.11986.000.86450.19690.53512.800.0625331.31250.18138.591.25780.25000.6006.020.0585;系統(tǒng)電路結(jié)構(gòu)拓撲圖如下:I00MW圖1 WSCC 3機9節(jié)點系統(tǒng)（所有參數(shù)以100MVA為基準值的標幺值）2.潮流計算首先進行潮流計算，采用牛頓拉夫遜迭代法，電力系統(tǒng)潮流計算是電力系統(tǒng)運行和規(guī) 劃中最基本和最經(jīng)常的計算，其任務(wù)是在已知某些運行參數(shù)的情況下，計算出系統(tǒng)中全部 的運行參數(shù)，一般來說，各個母線j0.06250_0085+j0.072B/2-J0.0745_ 35

6、MVARjO.O586廠、230/13.830 0119壩.1008B/2=jOJO45ulcT笙0OTZM / 5MJoirzxQMor6mo125 MW50 MVARso.orsZ60OESO6s0rzs90 MW30 MVAR6OR2E于永生電力系統(tǒng)仿真作業(yè)3所供負荷的功率是已知的，各個節(jié)點電壓是未知的（平衡節(jié)點除外），可以根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)形成節(jié)點導(dǎo)納矩陣，然后由節(jié)點導(dǎo)納矩陣和網(wǎng)絡(luò)拓撲結(jié)構(gòu)列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計 算的問題就轉(zhuǎn)化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述 功率方程改寫成功率平衡方程，并對功率平衡方程求

7、偏導(dǎo)，得出對應(yīng)的雅可比矩陣，給未 知節(jié)點賦電壓初值，一般為額定電壓，將初值帶入功率平衡方程，得到功率不平衡量，這 樣由功率不平衡量、雅可比矩陣、節(jié)點電壓不平衡量（未知的）構(gòu)成了誤差方程，解誤差 方程，得到節(jié)點電壓不平衡量，節(jié)點電壓加上節(jié)點電壓不平衡量構(gòu)成新的節(jié)點電壓初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡 量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。牛頓拉夫遜算法修正方程W = -J &其中W是節(jié)點不平衡量向量，包括有功，無功，電壓；J是雅克比矩陣；AV是節(jié)點電壓修正量。令Vi =ejfi；Yij=GijjBij則直角坐標形式的功率

8、不平衡量方程為PQ節(jié)點：nn：R=Rs-/ *）廣）廣（GjfjBjej）=0j1jnnQ=Qis fi(GjejBijfj) ej(GjfjBj%) =0j1jwPV節(jié)點：nn. P =Ps-/(Gijej-Bijfj)-fj(GijfjBijej)=0jij 日日. v2=Vi：-Vi2=Vi： - -免免2fi2)極坐標形式的功率不平衡量方程nP =PsViVj(GijcosijBijsinij) =0jTnQ=Qis-V? Vj(Gjsin一Bjcos)=0j日雅可比矩陣J各元素的表達式H NJ =? ?N L,電力系統(tǒng)仿真作業(yè)于永生43.負荷等效和支路簡化然后求出支路電流，將發(fā)電機內(nèi)

9、電抗X加入系統(tǒng)導(dǎo)納矩陣，求出發(fā)電機內(nèi)電勢E。加入發(fā)電機內(nèi)節(jié)點后，系統(tǒng)導(dǎo)納矩陣變成12*12階的矩陣，并將負荷等效成阻抗。然后將支路導(dǎo)納矩陣分塊，如下：A E )Y = I9 DJ其中，A是3*3的方陣，E是3*9的矩陣，C是9*3的矩陣，D是9*9的方陣。Hij郵郵BgGUi:fjNj：ejMj-：fj=BjfiGjeLj= GijfiBje當j=i時:Hn：fi=Biie-Giifi -biNH半R=一Giie- Biifi -ai：eQ=Giie +Biifi-ai-fiLii其中,a八八(GT -Bjfj)bi 八八(Gf Bjj進行牛頓拉夫遜算法迭代后得到電壓幅值V和相角0。于永生電力

10、系統(tǒng)仿真作業(yè)5經(jīng)過網(wǎng)絡(luò)簡化得到故障前的3*3簡化導(dǎo)納矩陣電力系統(tǒng)仿真作業(yè)于永生6Ypre=A-E*(inv(D)*C(1)其中“inv ( D廣是MATLAB中D矩陣的求逆。故障中導(dǎo)納矩陣的第七行和第七列從矩陣中刪除，此時有fA E )Y =1e-7)P=zeros(nbus,1);%迭代初值P為零Q=zeros(nbus,1);%迭代初值Q為零for i=1:nbusfor k=1:nbusP(i)=P(i)+V(i)*V(k)*(G(i,k)*cos(del(i)-del(k)+B(i,k)*sin(del(i)-del(k);Q(i)=Q(i)+V(i)*V(k)*(G(i,k)*sin

11、(del(i)-del(k)-B(i,k)*cos(del(i)-del(k);endenddPa=Psp-P;%額定有功功率減去流出或流入的有功功率dQa=Qsp-Q;%額定無功功率減去流出或流入的無功功率k=1;dQ=zeros(npq,1);%pq節(jié)點的無功功率，pq乘1的矩陣,6乘1for i=1:nbusif type(i)=1%pq節(jié)點dQ(k,1)=dQa(i);k=k+1;endenddP=dPa(2:nbus);Z=dP;dQ;%下面計算雅克比矩陣J1=zeros(nbus-1,nbus-1);%有功功率對功角求偏微分for i=1:(nbus-1)m=i+1;for k=1

12、:(nbus-1)n=k+1;if n=mfor n=1:nbusJ1(i,k)=J1(i,k)+V(m)*V(n)*(-G(m,n)*sin(del(m)-del(n)+B(m,n)*cos(del(m)-del(n);endJ1(i,k)=J1(i,k)-V(m)A2*B(m,m);elseJ1(i,k)=V(m)*V(n)*(G(m,n)*sin(del(m)-del(n)-B(m,n)*cos(del(m)-del(n);endendendJ2=zeros(nbus-1,npq);%有功功率對電壓求偏微分for i=1:(nbus-1)于永生電力系統(tǒng)仿真作業(yè)11m=i+1;for k=

13、1:npqn=pq(k);if n=mfor n=1:nbusJ2(i,k)=J2(i,k)+V(n)*(G(m,n)*cos(del(m)-del(n)+B(m,n)*sin(del(m)-del(n);endJ2(i,k)=J2(i,k)+V(m)*G(m,m);elseJ2(i,k)=V(m)*(G(m,n)*cos(del(m)-del(n)+B(m,n)*sin(del(m)-del(n);endendendJ3=zeros(npq,nbus-1);%無功功率對相角求偏微分for i=1:npqm=pq(i);for k=1:(nbus-1)n=k+1;if n=mfor n=1:n

14、busJ3(i,k)=J3(i,k)+V(m)*V(n)*(G(m,n)*cos(del(m)-del(n)+B(m,n)*sin(del(m)-del(n);endJ3(i,k)=J3(i,k)-V(m)A2*G(m,m);elseJ3(i,k)=V(m)*V(n)*(-G(m,n)*cos(del(m)-del(n)-B(m,n)*sin(del(m)-del(n);endendendJ4=zeros(npq,npq);%無功功率對電壓求偏微分for i=1:npqm=pq(i)；for k=1:npqn=pq(k)；if n=mfor n=1:nbusJ4(i,k)=J4(i,k)+V(

15、n)*(G(m,n)*sin(del(m)-del(n)-B(m,n)*cos(del(m)-del(n);endJ4(i,k)=J4(i,k)-V(m)*B(m,m);elseJ4(i,k)=V(m)*(G(m,n)*sin(del(m)-del(n)-B(m,n)*cos(del(m)-del(n);endendendJ=J1 J2;J3 J4;電力系統(tǒng)仿真作業(yè)于永生12%X=inv(J)*Z;X=JZ;dTh=X(1:nbus-1);dV=X(nbus:end);del(2:nbus)=dTh+del(2:nbus);k=1;for i=2:nbusif type(i)=1V(i)=V(

16、i)+dV(i-3);k=k+1;endendIter=Iter+1;Tol=max(abs(Z);endVm = pol2rect(V ,del);%極坐標轉(zhuǎn)換成直角坐標Del = 180/pi*del;%弧度轉(zhuǎn)化成角度Iij = zeros(nb,nb);%支路電流for m = 1:nl%每一條支路都過一遍p = fb(m); q = tb(m);Iij(p,q) = -(Vm(p) - Vm(q)*Y(p,q); % Y(m,n) = -y(m,n). iij(q,p) = -iij(p,q);endIij = sparse(Iij);%線路數(shù)據(jù)%首端末端電阻電抗電納(:1/2)變壓器

17、非標準變比L=450.010.0850.0881460.0170.0920.0791570.0320.1610.1531690.0390.170.1791780.00850.0720.07451890.01190.10080.105411400.0576012700.0625013900.058601;Yzengjia=zeros(12,12);Zci=j*0.0608 j*0.1198 j*0.1813;Yci=1./Zci;for k=1:3于永生電力系統(tǒng)仿真作業(yè)13Yzengjia(k,k)=Yci(k);endfor k=1:9for m=1:9Yzengjia(3+k,3+m)=Y(

18、k,m);endendYzengjia(1,4)=-Yzengjia(1,1);Yzengjia(4,1)=Yzengjia(1,4);Yzengjia(2,5)=-Yzengjia(2,2);Yzengjia(5,2)=Yzengjia(2,5);Yzengjia(3,6)=-Yzengjia(3,3);Yzengjia(6,3)=Yzengjia(3,6);Yzengjia(4,4)=Yzengjia(4,4)+Yzengjia(1,1);Yzengjia(5,5)=Yzengjia(5,5)+Yzengjia(2,2);Yzengjia(6,6)=Yzengjia(6,6)+Yzengj

19、ia(3,3);%求負荷等效導(dǎo)納Y5=N(5,6)/(V(5)A2-j*(N(5,7)/(V(5)A2);Y6=N(6,6)/(V(6)A2-j*(N(6,7)/(V(5)A2);Y8=N(8,6)/(V(8)A2-j*(N(8,7)/(V(5)A2);Yzengjia(8,8)=Yzengjia(8,8)+Y5;Yzengjia(9,9)=Yzengjia(9,9)+Y6;Yzengjia(11,11)=Yzengjia(11,11)+Y8;A=zeros(3,3);for m=1:3for n=1:3A(m,n)=Yzengjia(m,n);endendE=zeros(3,9);for m

20、=1:3for n=4:12E(m,n-3)=Yzengjia(m,n);endendC=zeros(9,3);for m=4:12for n=1:3C(m-3,n)=Yzengjia(m,n);endendD=zeros(9,9);電力系統(tǒng)仿真作業(yè)于永生14for m=4:12for n=4:12D(m-3,n-3)=Yzengjia(m,n);endendYpre=A-E*(inv(D)*C;%得到故障前的3*3矩陣Yzengjiadur=Yzengjia;Yzengjiadur(:,10)=;Yzengjiadur(10,:)=;Adur=A;Edur=zeros(3,8);for m=

21、1:3for n=1:8Edur(m,n)=Yzengjiadur(m,n+3);%因為Edur的第四到第九列是零，而所以這里為了方便，直接取E的前六列endendCdur=zeros(8,3);for m=1:8for n=1:3Cdur(m,n)=Yzengjiadur(m+3,n);%因為Cdur的第四到第九行是零，而所以這里為了方便，直接取C的前六行endendDdur=zeros(8,8);for m=1:8for n=1:8Ddur(m,n)=Yzengjiadur(m+3,n+3);endendYdur=Adur-Edur*(inv(Ddur)*Cdur;%故障中的3*3矩陣Yz

22、engjiaaft=Yzengjia;temp=Yzengjiaaft(8,10);Yzengjiaaft(8,10)=0;Yzengjiaaft(8,8)=Yzengjiaaft(8,8)+temp;Yzengjiaaft(10,8)=0;Yzengjiaaft(10,10)=Yzengjiaaft(10,10)+temp;Aaft=zeros(3,3);for m=1:3for n=1:3Aaft(m,n)=Yzengjiaaft(m,n);endendEaft=zeros(3,9);for m=1:3for n=4:12于永生電力系統(tǒng)仿真作業(yè)15Eaft(m,n-3)=Yzengjiaa

23、ft(m,n);endendCaft=zeros(9,3);for m=4:12for n=1:3Caft(m-3,n)=Yzengjiaaft(m,n);endendDaft=zeros(9,9);for m=1:9for n=1:9Daft(m,n)=Yzengjiaaft(m+3,n+3);endendYaft=Aaft-Eaft*(inv(Daft)*Caft;%故障后的3*3矩陣%求電磁功率Enei=ones(3,1);%發(fā)電機內(nèi)電勢Pe=ones(3,1);Enei(2)=Vm(2)+Zci(2)*Iij(2,7);Enei(1)=Vm(1)+Zci(1)*(-Iij(4,1);E

24、nei(3)=Vm(3)+Zci(3)*Iij(3,9);polEnei1=rect2pol(Enei(1);%直角坐標轉(zhuǎn)換成極坐標后得到角度polE nei2=rect2pol(Enei(2);polE nei3=rect2pol(Enei(3);Pe(1)=real(Enei(1)*(-Iij(4,1);%故障前的電氣功率Pe(2)=real(Ene i(2) *(Iij(2,7);Pe(3)=real(Enei(3)*(Iij(3,9);delta=zeros(3,1);delta(1)=po lE nei1(2);delta(2)=po丘nei2(2);delta(3)=p ol En

25、ei3(2);Tj=47.28;12.8;6.02;Pm0=Pe;%機械功率=故障前的電磁功率X0=ones(3,1).*(2*pi*60),delta;Y_Gen=Ydur;E=abs(Enei);Time_s=0,0.083;%0.083s后切除故障t1,Xout1=ode45(t,X) Gen_fw(t,X,Y_Gen,E,Pm0,Tj),Time_s,X0);Time_s=0.083,2;Y_Gen=Yaft;X0=Xout1(end,:);t2,Xout2=ode45(t,X) Gen_fw(t,X,Y_Gen,E,Pm0,Tj),Time_s,X0);XOUT=Xout1;Xout

26、2;T=t1;t2;電力系統(tǒng)仿真作業(yè)于永生16XOUT=XOUT*180/pi;figure,plot(T,XOUT(:,4),r,T,XOUT(:,5),g,T,XOUT(:,6),y);figure,plot(T,XOUT(:,5)-XOUT(:,4),r,T,XOUT(:,6)-XOUT(:,4),Y);(2) .極坐標轉(zhuǎn)換成直角坐標函數(shù)pol2rect(V,del)function rect = pol2rect(rho,theta)rect = rho.*cos(theta) + j*rho.*sin(theta);(3) .直角坐標轉(zhuǎn)換成極坐標函數(shù)rect2pol(Z)functi

27、on pol=rect2pol(Z)t=real(Z);k=imag(Z);pol(1)=(tA2+kA2)A(1/2);%模pol(2)=atan(k/t);(4) .求解微分方程所用的得到微分量的函數(shù)Gen_fw(t,X,Y_Gen,E,Pm0,Tj)function Xp = Gen_fw(t,X,Y_Gen,E,Pm0,Tj)Wn=ones(3,1).*2*pi*60;D=zeros(3,1);dW=zeros(3,1);dDeta=zeros(3,1);W=X(1:3);Deta=X(4:6);Gen_U=E.*exp(1i*Deta);Pe=real(Gen_U.*conj(Y_Gen*Gen_U);dW=(Pm0-Pe-D.*W)./Tj.*Wn;dDeta=W-Wn;Xp=dW;dDeta;