故障切除時間與三相快速重合閘對電力系統(tǒng)穩(wěn)定性的影響仿真(含MATLAB程序)

1、影響電力系統(tǒng)暫態(tài)穩(wěn)定性的因素1 研究對象1.1 系統(tǒng)模型研究影響系統(tǒng)暫態(tài)穩(wěn)定性的因素，首先要確定一個系統(tǒng)模型，本例選取參數(shù)可靠的美國西部電網(wǎng)等值模型SWCC-9，該系統(tǒng)為三機(jī)九節(jié)點(diǎn)系統(tǒng)，如圖1.1所示。圖1.1 WSCC-9系統(tǒng)模型1.2 系統(tǒng)參數(shù)1.2.1 節(jié)點(diǎn)參數(shù)表1.1 節(jié)點(diǎn)已知參數(shù)節(jié)點(diǎn)類型電壓幅值電壓角度發(fā)電機(jī)有功發(fā)電機(jī)無功負(fù)荷有功負(fù)荷無功1V1.04000.71640.2705002PV1.0251.63000.0665003PV1.0250.8500-0.1086004PQ00005PQ001.25000.50006PQ000.90000.30007PQ00008PQ001.000

2、00.35009PQ0000上表中發(fā)電機(jī)有功、無功出力和負(fù)荷的有功無功功率均為以100MVA為基準(zhǔn)時的標(biāo)幺值。1.2.2 支路參數(shù)表1.2 支路參數(shù)首節(jié)點(diǎn)末節(jié)點(diǎn)電阻電抗電納一半450.01000.08500.0880460.01700.09200.0790570.03200.16100.1530690.03900.17000.1790780.00850.07200.0745890.01190.10080.1045140.00000.05760.0000270.00000.06250.0000390.00000.05860.0000上表中所有的參數(shù)均為標(biāo)幺值，對于變壓器支路。最后三行表示三臺變壓

3、器參數(shù)，已經(jīng)計(jì)算出變壓器的等效電抗并直接在表格中給出。1.2.3 發(fā)電機(jī)參數(shù)對于發(fā)電機(jī)，采用二階經(jīng)典模型，并對系統(tǒng)作如下假設(shè)：(1) 輸入的機(jī)械功率保持恒定；(2) 忽略阻尼效應(yīng)；(3) 負(fù)荷采用恒阻抗模型。表1.3 發(fā)電機(jī)參數(shù)發(fā)電機(jī)XdXd10.14600.060823.6420.89580.11986.4031.31250.18133.01以上阻抗參數(shù)均以標(biāo)幺值表示，額定轉(zhuǎn)速下存儲的能量（H）也轉(zhuǎn)化為以100MVA為基準(zhǔn)的標(biāo)幺值。2 待研究的影響因素及其仿真流程設(shè)計(jì)2.1 故障切除時間的影響由等面積法則可在理論上分析得出如下結(jié)論：系統(tǒng)發(fā)生短路故障后，故障切除的越快，越有利于提高系統(tǒng)的暫態(tài)穩(wěn)

4、定性，且存在臨界切除時間。假設(shè)系統(tǒng)0時刻以前處于穩(wěn)定運(yùn)行狀態(tài)，0時刻在線路5-7上靠近7的母線出口處發(fā)生三相短路，故障切除時間為tc，仿真三臺發(fā)電機(jī)的轉(zhuǎn)速與功角變化曲線，通過觀察曲線判斷系統(tǒng)是否失去穩(wěn)定。不斷改變故障時間，求出恰好能使系統(tǒng)維持穩(wěn)定的故障切除時間，該時間就是臨界切除時間。2.2 三相快速重合閘的影響由等面積法則分析可知，在系統(tǒng)發(fā)生短路故障后，在減速過程中自動重合閘動作將會增大加速面積，可能使原來無法維持暫態(tài)穩(wěn)定的系統(tǒng)保持穩(wěn)定。但是重合閘要求的時間比較苛刻，如果重合閘在系統(tǒng)已經(jīng)失去穩(wěn)定才動作則對系統(tǒng)的暫態(tài)穩(wěn)定起不到作用，因此要求時間通常很短。這時如果采用單相重合閘由于有潛供電流的影

5、響，單相重合閘的動作時間不能太快，因此本例只考慮三相重合閘。在2.1的基礎(chǔ)上，tar時刻自動重合閘動作，調(diào)整參數(shù)，做對比仿真實(shí)驗(yàn)，一次不投入重合閘，另一次投入重合閘，通過判斷系統(tǒng)是否失去穩(wěn)定來得出重合閘對系統(tǒng)暫態(tài)穩(wěn)定性是否有影響。然后，固定故障切除時間tc為某一個定值，改變tar的值并觀察系統(tǒng)是否失去穩(wěn)定來確定自動重合閘的最長動作延時。2.3 仿真流程設(shè)計(jì)首先要求出系統(tǒng)的穩(wěn)態(tài)運(yùn)行參數(shù)，即系統(tǒng)潮流分布，這些參數(shù)作為暫態(tài)過程的初始值。潮流計(jì)算采用Newton-Raphson迭代法，求出各個母線的電壓和角度。2.3.1 發(fā)電機(jī)初態(tài)在動態(tài)分析時，首先應(yīng)將發(fā)電機(jī)和負(fù)荷用相應(yīng)的模型等效。本文將發(fā)電機(jī)等效為

6、二階經(jīng)典模型，將負(fù)荷等效為恒阻抗負(fù)荷。圖2.1 發(fā)電機(jī)等效模型本系統(tǒng)中具有三臺發(fā)電機(jī)，因此會引入三個內(nèi)節(jié)點(diǎn)。因此節(jié)點(diǎn)導(dǎo)納矩陣Y將增廣到12階。可表示為：在中： 在和中： 在中： 將網(wǎng)絡(luò)等效在發(fā)電機(jī)內(nèi)節(jié)點(diǎn)構(gòu)成的網(wǎng)絡(luò)中，可以得到一個3階的降階節(jié)點(diǎn)導(dǎo)納矩陣，該矩陣可由以下關(guān)系解出： 對于故障前和故障后的降階節(jié)點(diǎn)導(dǎo)納矩陣，利用式可以計(jì)算出對應(yīng)的降階節(jié)點(diǎn)導(dǎo)納矩陣；對于故障中的節(jié)點(diǎn)導(dǎo)納矩陣，在12階的增廣矩陣中去掉故障母線所在的那一行和那一列，利用下式計(jì)算降階矩陣。 發(fā)電機(jī)初態(tài)中包含初始的電壓幅值和功角，根據(jù)圖2.1，可以求解發(fā)電機(jī)1、2、3的初始狀態(tài)。 2.3.2 列寫發(fā)電機(jī)動態(tài)方程發(fā)電機(jī)采用經(jīng)典模型，

7、其動態(tài)方程為： 三臺發(fā)電機(jī)共有六個狀態(tài)量，定義： 得到： 2.3.3 求解發(fā)電機(jī)狀態(tài)利用龍格庫塔數(shù)值積分方法，對發(fā)電機(jī)狀態(tài)方程積分，求解發(fā)電機(jī)的狀態(tài)。認(rèn)為在故障過程中，發(fā)電機(jī)內(nèi)電壓的幅值是不變的，只有功角改變；認(rèn)為發(fā)電機(jī)的機(jī)械功率是不變的，其值等于故障前的發(fā)電機(jī)電磁功率。3 仿真過程與結(jié)果3.1 潮流計(jì)算結(jié)果表3.1 潮流計(jì)算結(jié)果節(jié)點(diǎn)類型電壓幅值電壓角度發(fā)電機(jī)有功發(fā)電機(jī)無功負(fù)荷有功負(fù)荷無功1V1.04000.00000.71640.2705002PV1.02509.28001.63000.0665003PV1.02504.66480.8500-0.1086004PQ1.0258-2.21680

8、0005PQ0.9956-3.9888001.25000.50006PQ1.0127-3.6874000.90000.30007PQ1.02583.719700008PQ1.01590.7275001.00000.35009PQ1.03241.966700003.2 發(fā)電機(jī)初態(tài)表3.2 發(fā)電機(jī)初態(tài)發(fā)電機(jī)內(nèi)電壓功角11.05662.271621.050219.731631.017013.16643.3 故障切除時間影響暫態(tài)穩(wěn)定性仿真假設(shè)0時刻發(fā)生故障，經(jīng)過5工頻周期后切除故障，即tc = 0.083 s，重合閘不投入。得到三臺發(fā)電機(jī)的轉(zhuǎn)速和功角變化曲線。圖3.1 發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.2 發(fā)電

9、機(jī)功角變化曲線圖3.1和圖3.2分別為該過程中三臺發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真總時間為4 s。從曲線中可以看出系統(tǒng)沒有失去暫態(tài)穩(wěn)定?，F(xiàn)在改變故障切除時間tc為12個工頻周期，即0.2 s，三臺發(fā)電機(jī)的動態(tài)曲線分別如下。圖3.3發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.4發(fā)電機(jī)功角變化曲線圖3.3和圖3.4分別為該過程中三臺發(fā)電機(jī)的轉(zhuǎn)速變化曲線和功角變化曲線，仿真總時間為4 s。從曲線中可以明顯看出系統(tǒng)已經(jīng)失去暫態(tài)穩(wěn)定。因此可以得出故障切除時間會影響系統(tǒng)的暫態(tài)穩(wěn)定性，并且時間越短越有利于提高系統(tǒng)的暫態(tài)穩(wěn)定性。通過改變故障切除時間，可以求出臨界故障切除時間tcc = 0.16 s，這時的發(fā)電機(jī)暫態(tài)曲線如

10、下：圖3.5 臨界切除時間下的發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.6 臨界切除時間下的發(fā)電機(jī)功角變化曲線3.4 重合閘影響暫態(tài)穩(wěn)定性仿真圖3.7 tc = 0.17 s未投入重合閘時發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.8 tc = 0.17 s未投入重合閘時發(fā)電機(jī)功角變化曲線取故障切除時間tc = 0.17 s，由3.3節(jié)可知，該時間大于臨界故障切除時間，因此暫態(tài)穩(wěn)定不能維持，發(fā)電機(jī)暫態(tài)曲線如圖3.7和圖3.8所示。圖3.9 tc = 0.17 s投入重合閘時發(fā)電機(jī)轉(zhuǎn)速變化曲線圖3.10 tc = 0.17 s投入重合閘時發(fā)電機(jī)功角變化曲線現(xiàn)將自動重合閘投入，重合閘動作延時為0.25 s，即動作時刻為tar = 0.

11、42 s，這時發(fā)電機(jī)的轉(zhuǎn)速和功角變化曲線分別如圖3.9和圖3.10所示。從圖中課明顯看出此時系統(tǒng)沒有失去暫態(tài)穩(wěn)定性，說明三相快速重合閘能夠提高系統(tǒng)的暫態(tài)穩(wěn)定性。附 錄MATLAB仿真程序代碼：%生成節(jié)點(diǎn)導(dǎo)納矩陣 Ypofl = zeros(9); for n = 1:9 Ypofl(BRANCH(n,1),BRANCH(n,1) = Ypofl(BRANCH(n,1),BRANCH(n,1) + 1j*BRANCH(n,5) + 1/(BRANCH(n,3) + 1j*BRANCH(n,4); Ypofl(BRANCH(n,2),BRANCH(n,2) = Ypofl(BRANCH(n,2),

12、BRANCH(n,2) + 1j*BRANCH(n,5) + 1/(BRANCH(n,3) + 1j*BRANCH(n,4); Ypofl(BRANCH(n,1),BRANCH(n,2) = Ypofl(BRANCH(n,1),BRANCH(n,2) - 1/(BRANCH(n,3) + 1j*BRANCH(n,4); Ypofl(BRANCH(n,2),BRANCH(n,1) = Ypofl(BRANCH(n,2),BRANCH(n,1) - 1/(BRANCH(n,3) + 1j*BRANCH(n,4); end G = real(Ypofl); B = imag(Ypofl); tole

13、rance = 1e-10; tol = 1; LoopCount = 0; delta_y = zeros(14,1); H = zeros(8); N = zeros(8,6); M = zeros(6,8); L = zeros(6); %計(jì)算潮流 while tol > tolerance % Delta_y for m = 2:9 sum = 0; for n = 1:9 sum = sum + BUS(m,3) * BUS(n,3) * (G(m,n) * cos(BUS(m,4) - BUS(n,4) + B(m,n) * sin(BUS(m,4) - BUS(n,4);

14、end delta_y(m-1) = BUS(m,5) - BUS(m,7) - sum; end for m = 4:9 sum = 0; for n = 1:9 sum = sum + BUS(m,3) * BUS(n,3) * (G(m,n) * sin(BUS(m,4) - BUS(n,4) - B(m,n) * cos(BUS(m,4) - BUS(n,4); end delta_y(m+5) = BUS(m,6) - BUS(m,8) - sum; end % H for m = 2:9 for n = 2:9 if m = n sum = 0; for k = 1:9 sum =

15、 sum + BUS(m,3) * BUS(k,3) * (G(m,k) * sin(BUS(m,4) - BUS(k,4) - B(m,k) * cos(BUS(m,4) - BUS(k,4); end H(m-1,m-1) = BUS(m,3)2 * B(m,m) + sum; else H(m-1,n-1) = - BUS(m,3) * BUS(n,3) * (G(m,n) * sin(BUS(m,4) - BUS(n,4) - B(m,n) * cos(BUS(m,4) - BUS(n,4); end end end % N for m = 2:9 for n = 4:9 if m =

16、 n %行列不等時 N(m-1,n-3) = - BUS(m,3) * BUS(n,3) * (G(m,n) * cos(BUS(m,4) - BUS(n,4) + B(m,n) * sin(BUS(m,4) - BUS(n,4); end end end for m = 2:7 sum = 0; for n = 1:9 sum = sum + BUS(m,3) * BUS(n,3) * (G(m,n) * cos(BUS(m,4) - BUS(n,4) + B(m,n) * sin(BUS(m,4) - BUS(n,4); end N(m-1,m-1) = - BUS(m,3)2 * G(m

17、,m) - sum; end % M for m = 4:9 for n = 2:9 if m = n sum = 0; for k = 1:9 sum = sum + BUS(m,3) * BUS(k,3) * (G(m,k) * cos(BUS(m,4) - BUS(k,4) + B(m,k) * sin(BUS(m,4) - BUS(k,4); end M(m-3,m-3) = BUS(m,3)2 * G(m,m) - sum; else M(m-3,n-1) = BUS(m,3) * BUS(n,3) * (G(m,n) * cos(BUS(m,4) - BUS(n,4) + B(m,

18、n) * sin(BUS(m,4) - BUS(n,4); end end end % L for m = 4:9 for n = 4:9 if m = n sum = 0; for k = 1:9 sum = sum + BUS(m,3) * BUS(k,3) * (G(m,k) * sin(BUS(m,4) - BUS(k,4) - B(m,k) * cos(BUS(m,4) - BUS(k,4); end L(m-3,m-3) = BUS(m,3)2 * B(m,m) - sum; else L(m-3,n-3) = - BUS(m,3) * BUS(n,3) * (G(m,n) * s

19、in(BUS(m,4) - BUS(n,4) - B(m,n) * cos(BUS(m,4) - BUS(n,4); end end end % 雅可比矩陣 JACOB = H N;M L; % 修正量 delta_x = - JACOB delta_y; % 修正電壓幅值和相位 for m = 4:9 BUS(m,3) = BUS(m,3) + delta_x(m+5) * BUS(m,3); end for m = 2:9 BUS(m,4) = BUS(m,4) + delta_x(m-1); end % 計(jì)算精度 tol = abs(max(delta_y); LoopCount = L

20、oopCount + 1; end BUS(:,4) = BUS(:,4) * 180/pi; disp('POWERFLOW ',date); disp(BUS); % 故障分析 % 簡化的Y陣 % 故障前 % 計(jì)算增廣Y陣 Y_pf33 = zeros(3); Y_pf39 = zeros(3,9); Y_pf99 = Ypofl; Y_pf33(1,1) = 1/(1j*GEN(1,3); Y_pf33(2,2) = 1/(1j*GEN(2,3); Y_pf33(3,3) = 1/(1j*GEN(3,3); Y_pf39(1,1) = -1/(1j*GEN(1,3); Y

21、_pf39(2,2) = -1/(1j*GEN(2,3); Y_pf39(3,3) = -1/(1j*GEN(3,3); Y_pf93 = transpose(Y_pf39); Y_pf99(1,1) = Y_pf99(1,1) + 1/(1j*GEN(1,3); Y_pf99(2,2) = Y_pf99(2,2) + 1/(1j*GEN(2,3); Y_pf99(3,3) = Y_pf99(3,3) + 1/(1j*GEN(3,3); Y_pf99(5,5) = Y_pf99(5,5) + BUS(5,7)/BUS(5,3)2 - 1j*BUS(5,8)/BUS(5,3)2; Y_pf99(

22、6,6) = Y_pf99(6,6) + BUS(6,7)/BUS(6,3)2 - 1j*BUS(6,8)/BUS(6,3)2; Y_pf99(8,8) = Y_pf99(8,8) + BUS(8,7)/BUS(8,3)2 - 1j*BUS(8,8)/BUS(8,3)2; Yex = Y_pf33 Y_pf39; Y_pf93 Y_pf99; Yrpf = Y_pf33 - Y_pf39 / Y_pf99 * Y_pf93; % 故障中 Y_df33 = Y_pf33; Y_df38 = Y_pf39; Y_df38(:,7) = ; Y_df83 = transpose(Y_df38); Y

23、_df88 = Y_pf99; Y_df88(7,:) = ; Y_df88(:,7) = ; Yrdf = Y_df33 - Y_df38 / Y_df88 * Y_df83; % 故障后 Y_af33 = Y_pf33; Y_af39 = Y_pf39; Y_af93 = Y_pf93; Y_af99 = Y_pf99; Y_af99(5,5) = Y_af99(5,5) - 1/(BRANCH(3,3) + 1j*BRANCH(3,4) - 1j*BRANCH(3,5); Y_af99(7,7) = Y_af99(7,7) - 1/(BRANCH(3,3) + 1j*BRANCH(3,4

24、) - 1j*BRANCH(3,5); Y_af99(5,7) = 0; Y_af99(7,5) = 0; Yraf = Y_af33 - Y_af39 / Y_af99 * Y_af93; disp('Before fault ',date); disp(Yrpf); disp(Yrdf); disp(Yraf); % 故障前發(fā)電機(jī)狀態(tài) tempAng = BUS(1:3,4)*pi/180; tempV1 = BUS(1:3,3).*exp(1j.*tempAng); tempV2 = tempV1 + conj(BUS(1:3,5) + 1j.*BUS(1:3,6) ./

25、 tempV1) .* 1j.*GEN(1:3,3); E_pf = sqrt(real(tempV2).2 + imag(tempV2).2); Delta_pf = atan(imag(tempV2)./real(tempV2).*180./pi; GEN(:,5) = real(tempV2.*conj(Yrpf*tempV2); disp('GenInit ',date); disp('E'''); disp(E_pf'); disp('delta' ); disp(Delta_pf'); % 故障中發(fā)電機(jī)

26、狀態(tài) t_cut = 0.17; %故障切除時間 t_end = 4; %仿真結(jié)束時間 t_ar = 0.42; %重合閘時間 odetol = 1e-10; options = odeset('RelTol',1e-10); tspan = 0 , t_cut; Delta_init = Delta_pf.*pi./180; x0_df = ones(3,1).*2.*pi.*60;Delta_init; Tdf_out , Xdf_out = ode45(t,x)Gen_Fun(t,x,GEN,Yrdf,E_pf),tspan,x0_df,options); Delta_df = Xdf_out(:,4:6).*180./pi; Omega_df = Xdf_out(:,1:3); % 故障后發(fā)電機(jī)狀態(tài) tspan = t_cut , t_ar; x0_af = Xdf_out(end,:); Taf_out , Xaf_out = ode45(t,x)Gen_Fun(t,x,GEN,Yra