故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序剖析Word下载.docx

1、781.00000.35009上表中发电机有功、无功出力和负荷的有功无功功率均为以100MVA为基准时的标幺值。1.2.2 支路参数表1.2 支路参数首节点末节点电阻电抗电纳一半0.01000.08500.08800.01700.09200.07900.03200.16100.15300.03900.17000.17900.00850.07200.07450.01190.10080.10450.00000.05760.06250.0586上表中所有的参数均为标幺值，对于变压器支路。最后三行表示三台变压器参数，已经计算出变压器的等效电抗并直接在表格中给出。1.2.3 发电机参数对于发电机，采用二

2、阶经典模型，并对系统作如下假设：（1） 输入的机械功率保持恒定；（2） 忽略阻尼效应；（3） 负荷采用恒阻抗模型。表1.3 发电机参数发电机XdXd0.14600.060823.640.89580.11986.401.31250.18133.01以上阻抗参数均以标幺值表示，额定转速下存储的能量（H）也转化为以100MVA为基准的标幺值。2 待研究的影响因素及其仿真流程设计2.1 故障切除时间的影响由等面积法则可在理论上分析得出如下结论：系统发生短路故障后，故障切除的越快，越有利于提高系统的暂态稳定性，且存在临界切除时间。假设系统0时刻以前处于稳定运行状态，0时刻在线路5-7上靠近7的母线出口处

3、发生三相短路，故障切除时间为tc，仿真三台发电机的转速与功角变化曲线，通过观察曲线判断系统是否失去稳定。不断改变故障时间，求出恰好能使系统维持稳定的故障切除时间，该时间就是临界切除时间。2.2 三相快速重合闸的影响由等面积法则分析可知，在系统发生短路故障后，在减速过程中自动重合闸动作将会增大加速面积，可能使原来无法维持暂态稳定的系统保持稳定。但是重合闸要求的时间比较苛刻，如果重合闸在系统已经失去稳定才动作则对系统的暂态稳定起不到作用，因此要求时间通常很短。这时如果采用单相重合闸由于有潜供电流的影响，单相重合闸的动作时间不能太快，因此本例只考虑三相重合闸。在2.1的基础上，tar时刻自动重合闸动

4、作，调整参数，做对比仿真实验，一次不投入重合闸，另一次投入重合闸，通过判断系统是否失去稳定来得出重合闸对系统暂态稳定性是否有影响。然后，固定故障切除时间tc为某一个定值，改变tar的值并观察系统是否失去稳定来确定自动重合闸的最长动作延时。2.3 仿真流程设计首先要求出系统的稳态运行参数，即系统潮流分布，这些参数作为暂态过程的初始值。潮流计算采用Newton-Raphson迭代法，求出各个母线的电压和角度。2.3.1 发电机初态在动态分析时，首先应将发电机和负荷用相应的模型等效。本文将发电机等效为二阶经典模型，将负荷等效为恒阻抗负荷。图2.1 发电机等效模型本系统中具有三台发电机，因此会引入三个

5、内节点。因此节点导纳矩阵Y将增广到12阶。可表示为：在中： 和将网络等效在发电机内节点构成的网络中，可以得到一个3阶的降阶节点导纳矩阵，该矩阵可由以下关系解出：对于故障前和故障后的降阶节点导纳矩阵，利用式可以计算出对应的降阶节点导纳矩阵；对于故障中的节点导纳矩阵，在12阶的增广矩阵中去掉故障母线所在的那一行和那一列，利用下式计算降阶矩阵。发电机初态中包含初始的电压幅值和功角，根据图2.1，可以求解发电机1、2、3的初始状态。2.3.2 列写发电机动态方程发电机采用经典模型，其动态方程为：三台发电机共有六个状态量，定义：得到：2.3.3 求解发电机状态利用龙格库塔数值积分方法，对发电机状态方程积

6、分，求解发电机的状态。认为在故障过程中，发电机内电压的幅值是不变的，只有功角改变；认为发电机的机械功率是不变的，其值等于故障前的发电机电磁功率。3 仿真过程与结果3.1 潮流计算结果表3.1 潮流计算结果1.04001.02509.28004.66481.0258-2.21680.9956-3.98881.0127-3.68743.71971.01590.72751.03241.96673.2 发电机初态表3.2 发电机初态内电压功角1.05662.27161.050219.73161.017013.16643.3 故障切除时间影响暂态稳定性仿真假设0时刻发生故障，经过5工频周期后切除故障，即

7、tc = 0.083 s，重合闸不投入。得到三台发电机的转速和功角变化曲线。图3.1 发电机转速变化曲线图3.2 发电机功角变化曲线图3.1和图3.2分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4 s。从曲线中可以看出系统没有失去暂态稳定。现在改变故障切除时间tc为12个工频周期，即0.2 s，三台发电机的动态曲线分别如下。图3.3发电机转速变化曲线图3.4发电机功角变化曲线图3.3和图3.4分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4 s。从曲线中可以明显看出系统已经失去暂态稳定。因此可以得出故障切除时间会影响系统的暂态稳定性，并且时间越短越有利

8、于提高系统的暂态稳定性。通过改变故障切除时间，可以求出临界故障切除时间tcc = 0.16 s，这时的发电机暂态曲线如下：图3.5 临界切除时间下的发电机转速变化曲线图3.6 临界切除时间下的发电机功角变化曲线3.4 重合闸影响暂态稳定性仿真图3.7 tc = 0.17 s未投入重合闸时发电机转速变化曲线图3.8 tc = 0.17 s未投入重合闸时发电机功角变化曲线取故障切除时间tc = 0.17 s，由3.3节可知，该时间大于临界故障切除时间，因此暂态稳定不能维持，发电机暂态曲线如图3.7和图3.8所示。图3.9 tc = 0.17 s投入重合闸时发电机转速变化曲线图3.10 tc = 0

9、.17 s投入重合闸时发电机功角变化曲线现将自动重合闸投入，重合闸动作延时为0.25 s，即动作时刻为tar = 0.42 s，这时发电机的转速和功角变化曲线分别如图3.9和图3.10所示。从图中课明显看出此时系统没有失去暂态稳定性，说明三相快速重合闸能够提高系统的暂态稳定性。附 录MATLAB仿真程序代码：%生成节点导纳矩阵 Ypofl = zeros（9）; for n = 1: Ypofl（BRANCH（n,1）,BRANCH（n,1） = Ypofl（BRANCH（n,1）,BRANCH（n,1） + 1j*BRANCH（n,5） + 1/（BRANCH（n,3） + 1j*BRANC

10、H（n,4）; Ypofl（BRANCH（n,2）,BRANCH（n,2） = Ypofl（BRANCH（n,2）,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） +

11、1j*BRANCH（n,4）; end G = real（Ypofl）; B = imag（Ypofl）; tolerance = 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）; %计算潮流 while tol tolerance % Delta_y for m = 2: sum = 0; sum = sum + BUS（m,3） * BUS（n,3） * （G（m,n） * cos（BUS（m,4） - BUS（n

12、,4） + B（m,n） * sin（BUS（m,4） - BUS（n,4）; delta_y（m-1） = BUS（m,5） - BUS（m,7） - sum; for m = 4: 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）; delta_y（m+5） = BUS（m,6） - BUS（m,8） - sum; % H for n = 2: if m = n for k = 1: sum = sum + BUS（m,3） * BUS

13、（k,3） * （G（m,k） * sin（BUS（m,4） - BUS（k,4） - B（m,k） * cos（BUS（m,4） - BUS（k,4）; 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）; % N for n = 4: if m = n %行列不等时 N（m-1,n-3） = - BUS（m,3） * BUS（n,3） * （G（

14、m,n） * cos（BUS（m,4） - BUS（n,4） + B（m,n） * sin（BUS（m,4） - BUS（n,4）; N（m-1,m-1） = - BUS（m,3）2 * G（m,m） - sum; % M 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）; M（m-3,m-3） = BUS（m,3）2 * G（m,m） - sum; M（m-3,n-1） = BUS（m,3） * BUS（n,3） * （G（m,n） *

15、cos（BUS（m,4） - BUS（n,4） + B（m,n） * sin（BUS（m,4） - BUS（n,4）; % L L（m-3,m-3） = BUS（m,3）2 * B（m,m） - sum; L（m-3,n-3） = - 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）; % 雅可比矩阵 JACOB = H N;M L; % 修正量 delta_x = - JACOB delta_y; % 修正电压幅值和相位 BUS（m,3） = BUS（m,3）

16、+ delta_x（m+5） * BUS（m,3）; BUS（m,4） = BUS（m,4） + delta_x（m-1）; % 计算精度 tol = abs（max（delta_y）; LoopCount = LoopCount + 1; BUS（:,4） = BUS（:,4） * 180/pi; disp（POWERFLOW ,date）; disp（BUS）; % 故障分析 % 简化的Y阵 % 故障前 % 计算增广Y阵 Y_pf33 = zeros（3）; Y_pf39 = zeros（3,9）; Y_pf99 = Ypofl; Y_pf33（1,1） = 1/（1j*GEN（1,3）;

17、 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_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/

18、（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（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 *

19、Y_pf93; % 故障中 Y_df33 = Y_pf33; Y_df38 = Y_pf39; Y_df38（:,7） = ; Y_df83 = transpose（Y_df38）; Y_df88 = Y_pf99; Y_df88（7,:） = ; Y_df88（: 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*BR

20、ANCH（3,4） - 1j*BRANCH（3,5）; Y_af99（7,7） = Y_af99（7,7） - 1/（BRANCH（3,3） + 1j*BRANCH（3,4） - 1j*BRANCH（3,5）; Y_af99（5,7） = 0; Y_af99（7,5） = 0; Yraf = Y_af33 - Y_af39 / Y_af99 * Y_af93;Before fault disp（Yrpf）; disp（Yrdf）; disp（Yraf）; % 故障前发电机状态 tempAng = BUS（1:3,4）*pi/180; tempV1 = BUS（1:3,3）.*exp（1j.*t

21、empAng）; tempV2 = tempV1 + conj（BUS（1:3,5） + 1j.*BUS（1:3,6） ./ 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）;GenInit disp（E）; disp（E_pfdelta ）; disp（Delta_pf % 故障中发电机状态 t_cut = 0

22、.17; %故障切除时间 t_end = 4; %仿真结束时间 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）.*18

23、0./pi; Omega_df = Xdf_out（:,1:3）; % 故障后发电机状态 tspan = t_cut , t_ar; x0_af = Xdf_out（end,: Taf_out , Xaf_out = ode45（t,x）Gen_Fun（t,x,GEN,Yraf,E_pf）,tspan,x0_af,options）; Delta_af = Xaf_out（: Omega_af = Xaf_out（: % 重合闸后发电机状态 tspan = t_ar , t_end; x0_ar = Xaf_out（end,: Tar_out , Xar_out = ode45（t,x）Gen_

24、Fun（t,x,GEN,Yrpf,E_pf）,tspan,x0_ar,options）; Delta_ar = Xar_out（: Omega_ar = Xar_out（: % 整合结果 T_out = Tdf_out ; Taf_out ; Tar_out; X_out = Xdf_out ; Xaf_out ; Xar_out; Delta_out = X_out（: Omega_out = X_out（: % 绘制曲线 figure（Name,delta curveNumberTitleoffplot（T_out,Delta_out（:,1）,-,T_out,Delta_out（:,2）,-.,3）,- xlabel（t （s） ylabel（delta_1,delta_2,delta_3 （Deg）omega curve plot（T_out,Omega_out）; xl