2020三機九節(jié)點電力系統(tǒng)建模與仿真

1、三機九節(jié)點電力系統(tǒng)建模與仿真報告學 院專 業(yè)姓 名學 號指導教師郵 箱提交日期一、摘要電力系統(tǒng)仿真計算己經成為電力系統(tǒng)設計、運行與控制中不可缺少的手段。通過設 置不同故障類型、不同故障地點運用仿真技術可以對電力系統(tǒng)的暫態(tài)穩(wěn)定進行分析。本文采用IEEE 3機9節(jié)點的經典多機模型，基于隱式梯形積分法對系統(tǒng)發(fā)生三相金屬性短路故障 進行仿真，分析系統(tǒng)在這種情況下的暫態(tài)穩(wěn)定。發(fā)電機模型采用經典的二階模型；負荷采用恒定阻抗負荷。在 Matlab2010上編寫程序進行調試和運行。電力系統(tǒng)是由不同類型的發(fā)電機組、多種電力負荷、 不同電壓等級的電力網絡等組成的十分龐大復雜的動力學系統(tǒng)。其暫態(tài)過渡過程不僅包括電磁

2、方面的過渡過程，而且還有機電方面的過渡過程。由此可見，電力系統(tǒng)的數(shù)學模型是一個強非線性的高維狀態(tài)方程 組。在動 態(tài)穩(wěn)定仿真中使用簡單的電力系統(tǒng)模型，通過仿真計算分析說明，此仿真方法可以進行簡單的電力系統(tǒng)暫態(tài)分析，對提高電力系統(tǒng)暫態(tài)穩(wěn)定具有重要意義。二、案例本次課程主要應用 P. M. Anderson and A. A. Fouad 編寫的Power System Controland Stability一書中所引用的 Western System Coordi nated Cou ncil (WSCC)三機九節(jié)點系統(tǒng)模型。系統(tǒng)電路結構拓撲圖如下：B荷負GsorI6.5XT圖2-1 3機9節(jié)點

3、系統(tǒng)系統(tǒng)數(shù)據其中，節(jié)點數(shù)據如下：節(jié)點號有無負載類型電壓相角有功負荷無功負荷有功出力無功出力電壓基準期望電壓N=1031.04000.000.000.0071.6027.0016.501.0402021.02500.000.000.00163.006.7018.001.0253021.02500.000.000.0085.00-10.9013.801.0254001.00000.000.000.000.000.00230.001.0265101.00000.00125.0050.000.000.000.000.9966101.00000.0090.0030.000.000.000.001.013

4、7001.00000.000.000.000.000.00230.001.0268101.00000.00100.0035.000.000.000.001.0169001.00000.000.000.000.000.00230.001.032;%支路數(shù)據% 從到電阻電抗容納類型變比B=140.00.05760.011270.00.06250.011390.00.05860.011450.0100.0850.17600460.0170.0920.15800570.0320.1610.30600690.0390.1700.35800780.00850.0720.14900890.01190.100

5、80.20900;發(fā)電機數(shù)據如下：% 發(fā)電機母線XdXdTd0XqXqTq0TjXfGe= 110.14600.06088.960.09690.0969047.280.0576220.89580.11986.000.86450.19690.53512.800.0625331.31250.18138.591.25780.25000.6006.020.0585;三、仿真框圖在仿真之前，首先，應明確仿真的所要到達的結果，即仿真目標：本此仿真的結果主要是得到發(fā)電機攻角、轉速隨時間變化的值，包括故障前、故障中、故障后。故障前，系統(tǒng)處 于穩(wěn)定狀態(tài)，發(fā)電機的攻角、轉速基本穩(wěn)定。而當系統(tǒng)發(fā)生故障以及故障切除，

6、系統(tǒng)結構拓 撲發(fā)生變化，系統(tǒng)的狀態(tài)也將隨時間發(fā)生變化，為了求取系統(tǒng)狀態(tài)的變化，我們通過對系統(tǒng)進行簡化建立數(shù)學模型，得到相關的代數(shù)一微分方程組，進行數(shù)值計算，從而得到系統(tǒng)狀態(tài)的隨時間的變化值。 此次仿真的系統(tǒng)以發(fā)電機二階經典模型來進行系統(tǒng)是數(shù)學建模，系統(tǒng)的狀態(tài)量為發(fā)電機攻角、發(fā)電機轉速。其次，當明確仿真目標后，我們就得明確大體的仿真框架流程。仿真框架流程如下：圖3-1仿真流程圖四、仿真模型在電力系統(tǒng)的機電暫態(tài)仿真中， 常根據實際要求的不同，采用不同時間尺度的仿真模型,而仿真算法和采用的模型有直接的關系，下面就本次仿真實例機電暫態(tài)過程的仿真模型及其仿真算法。一、潮流計算由于本文以三機九節(jié)點為模型，

7、假定節(jié)點一為參考節(jié)點，這樣就有2兩個發(fā)電機的 PV節(jié)點，6個PQ節(jié)點，未知量為8個節(jié)點（包括2個PV節(jié)點和6個PQ節(jié)點）的電壓相角，還有6個節(jié)點（PQ節(jié)點）的電壓幅值?？梢韵惹蟪鯵矩陣0 -LL3611iaD0 +IT.361110000oaoooi00d00 -i-|gB 00001000 -IF.,00000 +LL36LLia03. SOM -33.30851-1. 3452 +11- 00411-1.5422 +L0. &10?i00ai：-1.3052 +IB.604112. 5529 -17.3382i0-1.187S 十 5.75Li000-1.9422 +Hflfc5L07i0

8、3.2242 -15. MDSi000 +16. OOOOi00+ 6.9753102,8047 -35B4456i0a00C0-1.0171 +13.658010aD +1?.00-1.2820 + 5. 5BB2i0Cglursns S through 910a0a00+ 7.064810a0a0-1.2820+ 5.5SS2i-1.6171 +L3Pfi960i02.7722 -25.3032i-1. 155+ 9.784513551 + Si. m3i2,71-32.1539L圖4-1 Y矩陣然后，我們列寫方程，也就是利用各個節(jié)點的有功、無功功率的平衡關系，列寫14個功率平衡方程。這樣

9、就能使用牛頓一拉夫遜算法來求解這14個非線性方程。其中的關鍵是要計算出雅克比矩陣-16. 74300D0ID1& 7435000000-IS-QJ780aQQ爍QJW00000-4L 2SLLSSSS10.弼銅D00-14B02L02TI1-門開00naao?i5-If.68710001. iT6DL-1.280500n0-I66444D05,陰MJ, isro0-2-16.743500氏 2099014b34260a0.383400000ID14. ir-24.399410.2320000IS-0374QQE 94E7*DLO.2850-34.269D00,75100工頌2-L027I-L?

10、3a000-41. 3SSKII.SSSSio. g 制400-1.76013.7805a-2. 0200011.BOSS-16.6S71000-2.蜩004.Q0-,3Q9310,Q-5,9444-1. 630000-EL 383402.9511-0.9377C0g. 2099000000-2,4?r5- 610-L4330000-Q. S5QQQ-0-T5SD6 9S922-59?35-957ColunnE 02 thTouh. 14-1.63000000-fl. BWQD0D2. 020500009BLL0.9377a2. 427&-1.80101.U3000.9893-2, W?300

11、05. 8&4500005, SI.637u296014.加期&14. 1667-23.99410.緲50IQ. 2350-34,2CTI圖4-2雅克比矩陣電力然后計算出修正量。在設定精度和最大迭代次數(shù)的前提下進行迭代，直到滿足要求。網絡的節(jié)點功率方程可用一般形式表示如下牛頓拉夫遜算法修正方程W = -J AVJ是雅克比矩陣；AV是節(jié)點電其中W是節(jié)點不平衡量向量，包括有功，無功，電壓; 壓修正量。令V eijfi ;YjGj jBj則極坐標形式的功率不平衡量方程npRs ViVj (Gj cos ijBij sin j )0j inQ Qis ViVj (Gij sin ij Bij cos

12、ij ) 0j 1雅可比矩陣J各元素的表達式譽5二 B.e.-Giifi - 勺N廣簽5B虛N a 二 Sei:-Gf尼叫:Sf.GeH i十 Biifi _ aiLy =警 3：=B .e.-+ b.其中，an(Gj ejBij f j )i 1bin(Gij f jBij ej )i 1o進行牛頓拉夫遜算法得到潮流結果1- 000003.0000 l.MOOc00 7LM 27.000016.50C0：,04002, DOOO02.00D：i1.02500.1000 1610000 8.? 00018.0000I.QC503.000002, DQD：i 1.02500JS1400 GQJO

13、OO -10,30001X80001.02504,0000001.0258-0.038700 0 0230.0000M5. ml.DCOj0i.9956125,0000別MM00n0,93606.0C0D1.0C0301J12?-0.064450/100030.003000ij：,0l30?JC0D00 L J258OM00 0 0230.0000L02S0& 0C031.0C0301.01590J12?100.000035.C00000uL01509. 00030:1.3324匚 034300 0 0230.00001.0320圖4-3潮流結果、故障前中后僅含發(fā)電機內節(jié)點的導納矩陣 YpVp

14、 =D.65680,0701 +0. 630613455 - 2.9383i2871 + 1.512910. 2871 + 1.512910.20S6 + L 2256i0.2133 + L087910.2770 - 2. 366li0.4200 -0.2133 +2. 723911.087912096 +1. 22561 Yd5,48551C,07010.630610.1740 -2.?959iYa 二1.1386 - 2.296610.1290 + 0. 7O63i0. 1824 + 1. 063710. 1290 + 0.706310. 3744 - 2.016110.1921 + 1.

15、206710.1324 + 1.0637i0.1921 + L 206710. 2691 - 2. 3516i圖4-4故障前中后僅含發(fā)電機內節(jié)點的導納矩陣三、求解電磁功率得到故障前，故障中，故障后三個不同的導納矩陣后，就開始計算電磁功率和機械功率，機械功率等于穩(wěn)態(tài)的電磁功率中的有功分量。所以可以有Pe=real(E*l)如上中，E為發(fā)電機內電勢，I為從發(fā)電機流出的電流。但在參考文獻 Ramnarayan Patel, T. S. Bhatti and D. P . Kothari. MATLAB/Simulink-based transient stability analysis of a

16、multimachine power system 中給出的電磁功率計算公式為：n2 FeiEi GiiEi EjYj COS( ij i j )j 1j i穩(wěn)態(tài)情況下有，機械功率Pme=Pe四、求解運動方程發(fā)電機的運動方程可以寫成常微分方程組：2Hi dq r dt= Pmi E； G祈 +d Yg CO- $ + 巧)d8石a其中Pmi為第i個機組故障前穩(wěn)態(tài)的電磁功率。在本次仿真中Dj wj為零，即阻尼為零。t=0仿真開始，t=1時引入故障，1.083s后切除故障，t=3.083s仿真結束。求解運動方程后得到曲 線如下：五、結果分析上圖分別顯示了各臺發(fā)電機的轉子角與時間的關系曲線，顯示了發(fā)

17、電機轉速差的曲線，和2121、 313 1的曲線，由圖可以看到，最大角差21 為85,出現(xiàn)在t 0.4s處，無論是21還是31第二個搖擺都不大于第一個搖擺, 可見系統(tǒng)是穩(wěn)定的。六、程序代碼主程序：global Pm Yrun gen GenE%節(jié)點數(shù)據%節(jié)點號 有無負載 類型 電壓 相角 有功負荷 無功負荷 有功負荷 無功負荷 電壓基 準 期望電壓N=1031.04000.000.002021.02500.000.003021.02500.000.004001.00000.000.005101.00000.00125.006101.00000.0090.007001.00000.000.008

18、101.00000.00100.009001.00000.000.000.0071.6027.0016.501.0400.00163.006.7018.001.0250.0085.00-10.9013.801.0250.000.000.00230.00 1.02650.000.000.000.000.99630.000.000.000.001.0130.000.000.00230.001.02635.000.000.000.001.0160.000.000.00230.001.032;%支路數(shù)據% 從到電阻電抗容納類型 變比B=140.00.05760.011270.00.06250.0113

19、90.00.05860.011450.0100.0850.17600460.0170.0920.15800570.0320.1610.30600690.0390.1700.35800780.00850.0720.14900890.01190.10080.20900;%發(fā)電機數(shù)據% H MVA xd*10000 node xd xq xl xad xaq xf td0 rfgen=2364 247.5 608 1 0.1460 0.0969 0.0336 0.1124 0.0633 0.1483 8.96 0.0000439640 192.0 1198 2 0.8958 0.8645 0.052

20、1 0.8437 0.8124 0.9173 6.00 0.0004054301 128.0 1813 3 1.3125 1.2578 0.0742 1.2383 1.2836 1.3555 5.89 0.0006105;Y=zeros(9,9);% 導納矩陣for i=1:9a=B(i,1);b=B(i,2);if B(i,6)=0Y(a,b)=-1./(B(i,3)+B(i,4)*1i);Y(b,a)=-1./(B(i,3)+B(i,4)*1i);Y(a,a)=Y(a,a)+1./(B(i,3)+B(i,4)*j)+B(i,5)*1i./2;Y(b,b)=Y(b,b)+1./(B(i,3)

21、+B(i,4)*j)+B(i,5)*1i./2;elseY(a,b)=-1./(B(i,3)+1i*B(i,4)*B(i,6);Y(b,a)=-1./(B(i,3)+1i*B(i,4)*B(i,6);Y(a,a)=Y(a,a)+1./(B(i,3)+B(i,4)*1i);Y(b,b)=Y(b,b)+1./(B(i,3)+B(i,4)*j*B(i,6)A2)+B(i,5)*1i./2;endend %導納矩陣for T=1:100dP,dQ=Caoliu(N,Y);% 潮流J=Ykb(N,Y);% 雅克比矩陣U=zeros(6);for i=4:9U(i-3,i-3)=N(i,4);enddAn

22、gU=JdP;dQ;dAng=dAngU(1:8,1);dU=U*(dAngU(9:14,1);N(4:9,4)=N(4:9,4)-dU;N(2:9,5)=N(2:9,5)-dAng;if(max(abs(dU)0.00001)&(max(abs(dAng)0.00001)breakendEndYc,Yb,Ya,YY0,YY2,YY3=Ysim(gen,N,B,Y);GEgj=zeros(1,3);GenE=zeros(1,3);for i=1:3GenE(1,i)=abs(N(i,4)*exp(1i*N(i,5)+1i*gen(i,3)/10000*conj(N(i,8)/100+1i*N(

23、i,9)/100)/(N(i,4)* exp(1i*N(i,5);GEgj(1,i)=angle(N(i,4)*exp(1i*N(i,5)+1i*gen(i,3)/10000*conj(N(i,8)/100+1i*N(i,9)/100)/(N(i,4) *exp(1i*N(i,5);endYrun=zeros(3);Pm=zeros(1,3);for i=1:3Pm(1,i)=N(i,8)./100;endoptions=odeset(RElTOL,1e-10); % 設置精度XX=GEgj(1),GEgj(2),GEgj(3),2*pi*60*ones(1,3);tt=0;t0=1; tsp

24、an0=tt,t0; Yrun=Yc;T0,Y0=ode45(fzz,tspan0,XX,options);X0=Y0(end,:);t0=1; tc=1.083; tspan1=t0,tc; Yrun=Yb;T1,Y1=ode45(fzz,tspan1,X0,options);X1=Y1(end,:); tf=3.083; tspan2=tc,tf; Yrun=Ya;T2,Y2=ode45(fzz,tspan2,X1,options);delta1=Y0(:,1:3);delta2=Y1(:,1:3);delta3=Y2(:,1:3);E1=repmat(1.05662336337847,1

25、.05012618834491,1.01677591008947,45,1);E2=repmat(1.05662336337847,1.05012618834491,1.01677591008947,41,1);E3=repmat(1.05662336337847,1.05012618834491,1.01677591008947,369,1); e1=E1.*exp(delta1*1i);e2=E2.*exp(delta2*1i);e3=E3.*exp(delta3*1i);V1=-(YY0(4:e nd,4:e nd)A-1)*(YY0(4:e nd,1:3)*e1：V2 仁-(YY2(4

26、:e nd,4:e nd)A-1)*YY2(4:e nd,1:3)*e2：V2=V21(1:6,:);zeros(1,length(V21);V21(7:8,:);V3=-(YY3(4:e nd,4:e nd)A-1)*YY3(4:e nd,1:3)*e3：T=T0;T1;T2; Y12=Y0;Y1;Y2;V0=V1,V2,V3;subplot(3,2,4);plot(T,abs(V0);subplot(3,2,5);plot(T,unwrap(angle(V0)*180/pi);subplot(3,2,1);plot(T,Y12(:,1)*180/pi,T,Y12(:,2)*180/pi,T

27、,Y12(:,3)*180/pi); % 發(fā)電機功角subplot(3,2,2);plot(T,Y12(:,4),T,Y12(:,5),T,Y12(:,6); % 發(fā)電機轉速subplot(3,2,3);發(fā)電機攻角差plot(T,Y12(:,2)*180/pi-Y12(:,1)*180/pi,T,Y12(:,3)*180/pi-Y12(:,1)*180/pi) %雅克比矩陣： function J=Ykb(N,Y)H=zeros(8);N1=zeros(8,6);K=zeros(6,8);L=zeros(6);for i=2:9for j=2:9if i=jH(i-1,j-1)=-N(i,4)

28、*N(j,4)*(real(Y(i,j)*sin(N(i,5)-N(j,5)-imag(Y(i,j)*cos(N(i,5)-N(j,5);elseH(i-1,j-1)=(N(i,9)-N(i,7)/100+imag(Y(i,j)*(N(i,4)F2);endendend %Hfor i=2:9for j=4:9if i=jN1(i-1,j-3)=-N(i,4)*N(j,4)*(real(Y(i,j)*cos(N(i,5)-N(j,5)+imag(Y(i,j)*sin(N(i,5)-N(j,5); elseN1(i-1,j-3) =-(N(i,8)-N(i,6)/100-real(Y(i,j)*

29、(N(i,4)A2);endendend %Nfor i=4:9for j=2:9if i=jK(i-3,j-1)=N(i,4)*N(j,4)*(real(Y(i,j)*cos(N(i,5)-N(j,5)+imag(Y(i,j)*sin(N(i,5)-N(j,5); elseK(i-3,j-1)=-(N(i,8)-N(i,6)/100+real(Y(i,j)*(N(i,4)F2);endendend %Kfor i=4:9for j=4:9if i=jL(i-3,j-3)=-N(i,4)*N(j,4)*(real(Y(i,j)*sin(N(i,5)-N(j,5)-imag(Y(i,j)*cos

30、(N(i,5)-N(j,5);elseL(i-3,j-3)=-(N(i,9)-N(i,7)/100+imag(Y(i,j)*(N(i,4)A2);endendend %LJ=H N1;K L;潮流計算function dP,dQ=Caoliu(N,Y)dP=zeros(8,1);dQ=zeros(6,1);for i=2:9dP(i-1,1)=(N(i,8)-N(i,6)/100;endfor i=4:9dQ(i-3,1)=(N(i,9)-N(i,7)/100;endfor i=2:9for j =1:9dP(i-1,1)=dP(i-1,1)-N(i,4)*N(j,4)*(real(Y(i,j

31、)*cos(N(i,5)-N(j,5)+imag(Y(i,j) *sin(N(i,5)-N(j,5);endendfor i=4:9for j =1:9dQ(i-3,1)=dQ(i-3,1)-N(i,4)*N(j,4)*(real(Y(i,j)*sin(N(i,5)-N(j,5)-imag(Y(i,j) *cos(N(i,5)-N(j,5);endend故障前中后僅含發(fā)電機內節(jié)點的導納矩陣：function Yp,Yd,Ya,YY0,YY2,YY3 = Ysim(gen,N,B,Y) %故障前中后僅含發(fā)電機內節(jié)點的導納矩陣：% 此處顯示詳細說明Yp=zeros(3);Yd=zeros(3);Y

32、a=zeros(3); YY1=ones(12);YY2=zeros(11);YY3=zeros(12);YY1=zeros(3),zeros(3,9);zeros(9,3),Y; % 故障前增廣Y5=N(5,6)/100/(N(5,4)A2)-(N(5,7)/100/(N(5,4F2)*1i;Y6=N(6,6)/100/(N(6,4)A2)-(N(6,7)/100/(N(6,4F2)*1i;Y8=N(8,6)/100/(N(8,4)A2)-(N(8,7)/100/(N(8,4)A2)*1i; %負載等效導納for i=1:3YY1(i,i)= -(1/gen(i,3)*10000)*1i;

33、YY1(i,i+3)=(1/gen(i,3)*10000)*1i; YY1(i+3,i)=YY1(i,i+3); endfor i=1:3YY1(i+3,i+3)=YY1(i+3,i+3)-(1/gen(i,3)*10000)*1i; % 發(fā)電機節(jié)點修改endYY1(3+5,3+5)=YY1(3+5,3+5)+Y5; YY1(3+6,3+6)=YY1(3+6,3+6)+Y6;YY1(3+8,3+8)=YY1(3+8,3+8)+Y8; % 負載節(jié)點修改Ynn=zeros(3);Ynr=zeros(3,9);Ynr1=zeros(3,8);Yrn=zeros(9,3);Yrn1=zeros(8,3

34、); Yrr=zeros(9);Yrr1=zeros(8);Ynn=YY1(1:3,1:3);Ynr=YY1(1:3,4:12);Yrn=YY1(4:12,1:3);Yrr=YY1(4:12,4:12); Yp=Ynn-Yn r*(Yrd)*Yrn;YY0=YY1;YY3=YY1;YY3(8,8)=YY3(8,8)-(1./(B(6,3)+B(6,4)*1i)+B(6,5)*1i/2);YY3(10,10)=YY3(10,10)- (1./(B(6,3)+B(6,4)*1i)+B(6,5)*1i/2);YY3(8,10)=0;YY3(10,8)=0; % 故障后增廣YY1(10,:)=; YY1(:,10)=;YY2=YY1; % 故障中增廣Ynn=YY2(1:3,1:3);Ynr1=YY2(1:3,4:11);Yrn1=YY2(4:11,1:3);Yrr1=YY2(4:11,4:11); Yd=Ynn-Yn r1*(Yrr1A-1)*Yrn1;Ynn=YY3(1:3,1:3);Ynr=YY3(1:3,4:12);Yrn=YY3(4:12,1:3);Yrr=YY3(4:12,4:12); Ya=Ynn- Ynr*(YrrA-1)*