故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序剖析.docx
《故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序剖析.docx》由会员分享，可在线阅读，更多相关《故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序剖析.docx（21页珍藏版）》请在冰豆网上搜索。
故障切除时间与三相快速重合闸对电力系统稳定性的影响仿真含MATLAB程序剖析
影响电力系统暂态稳定性的因素
1研究对象
1.1系统模型
研究影响系统暂态稳定性的因素，首先要确定一个系统模型，本例选取参数可靠的美国西部电网等值模型SWCC-9，该系统为三机九节点系统，如图1.1所示。
图1.1WSCC-9系统模型
1.2系统参数
1.2.1节点参数
表1.1节点已知参数
节点
类型
电压幅值
电压角度
发电机有功
发电机无功
负荷有功
负荷无功
1
Vθ
1.040
0
0.7164
0.2705
0
0
2
PV
1.025
\
1.6300
0.0665
0
0
3
PV
1.025
\
0.8500
-0.1086
0
0
4
PQ
\
\
0
0
0
0
5
PQ
\
\
0
0
1.2500
0.5000
6
PQ
\
\
0
0
0.9000
0.3000
7
PQ
\
\
0
0
0
0
8
PQ
\
\
0
0
1.0000
0.3500
9
PQ
\
\
0
0
0
0
上表中发电机有功、无功出力和负荷的有功无功功率均为以100MVA为基准时的标幺值。
1.2.2支路参数
表1.2支路参数
首节点
末节点
电阻
电抗
电纳一半
4
5
0.0100
0.0850
0.0880
4
6
0.0170
0.0920
0.0790
5
7
0.0320
0.1610
0.1530
6
9
0.0390
0.1700
0.1790
7
8
0.0085
0.0720
0.0745
8
9
0.0119
0.1008
0.1045
1
4
0.0000
0.0576
0.0000
2
7
0.0000
0.0625
0.0000
3
9
0.0000
0.0586
0.0000
上表中所有的参数均为标幺值，对于变压器支路。
最后三行表示三台变压器参数，已经计算出变压器的等效电抗并直接在表格中给出。
1.2.3发电机参数
对于发电机，采用二阶经典模型，并对系统作如下假设：
（1）输入的机械功率保持恒定；
（2）忽略阻尼效应；
（3）负荷采用恒阻抗模型。
表1.3发电机参数
发电机
Xd
X’d
1
0.1460
0.0608
23.64
2
0.8958
0.1198
6.40
3
1.3125
0.1813
3.01
以上阻抗参数均以标幺值表示，额定转速下存储的能量（H）也转化为以100MVA为基准的标幺值。
2待研究的影响因素及其仿真流程设计
2.1故障切除时间的影响
由等面积法则可在理论上分析得出如下结论：
系统发生短路故障后，故障切除的越快，越有利于提高系统的暂态稳定性，且存在临界切除时间。
假设系统0时刻以前处于稳定运行状态，0时刻在线路5-7上靠近7的母线出口处发生三相短路，故障切除时间为tc，仿真三台发电机的转速与功角变化曲线，通过观察曲线判断系统是否失去稳定。
不断改变故障时间，求出恰好能使系统维持稳定的故障切除时间，该时间就是临界切除时间。
2.2三相快速重合闸的影响
由等面积法则分析可知，在系统发生短路故障后，在减速过程中自动重合闸动作将会增大加速面积，可能使原来无法维持暂态稳定的系统保持稳定。
但是重合闸要求的时间比较苛刻，如果重合闸在系统已经失去稳定才动作则对系统的暂态稳定起不到作用，因此要求时间通常很短。
这时如果采用单相重合闸由于有潜供电流的影响，单相重合闸的动作时间不能太快，因此本例只考虑三相重合闸。
在2.1的基础上，tar时刻自动重合闸动作，调整参数，做对比仿真实验，一次不投入重合闸，另一次投入重合闸，通过判断系统是否失去稳定来得出重合闸对系统暂态稳定性是否有影响。
然后，固定故障切除时间tc为某一个定值，改变tar的值并观察系统是否失去稳定来确定自动重合闸的最长动作延时。
2.3仿真流程设计
首先要求出系统的稳态运行参数，即系统潮流分布，这些参数作为暂态过程的初始值。
潮流计算采用Newton-Raphson迭代法，求出各个母线的电压和角度。
2.3.1发电机初态
在动态分析时，首先应将发电机和负荷用相应的模型等效。
本文将发电机等效为二阶经典模型，将负荷等效为恒阻抗负荷。
图2.1发电机等效模型
本系统中具有三台发电机，因此会引入三个内节点。
因此节点导纳矩阵Y将增广到12阶。
可表示为：
在
中：
在
和
中：
在
中：
将网络等效在发电机内节点构成的网络中，可以得到一个3阶的降阶节点导纳矩阵，该矩阵可由以下关系解出：
对于故障前和故障后的降阶节点导纳矩阵，利用式可以计算出对应的降阶节点导纳矩阵；对于故障中的节点导纳矩阵，在12阶的增广矩阵中去掉故障母线所在的那一行和那一列，利用下式计算降阶矩阵。
发电机初态中包含初始的电压幅值和功角，根据图2.1，可以求解发电机1、2、3的初始状态。
2.3.2列写发电机动态方程
发电机采用经典模型，其动态方程为：
三台发电机共有六个状态量，定义：
得到：
2.3.3求解发电机状态
利用龙格——库塔数值积分方法，对发电机状态方程积分，求解发电机的状态。
认为在故障过程中，发电机内电压的幅值是不变的，只有功角改变；认为发电机的机械功率是不变的，其值等于故障前的发电机电磁功率。
3仿真过程与结果
3.1潮流计算结果
表3.1潮流计算结果
节点
类型
电压幅值
电压角度
发电机有功
发电机无功
负荷有功
负荷无功
1
Vθ
1.0400
0.0000
0.7164
0.2705
0
0
2
PV
1.0250
9.2800
1.6300
0.0665
0
0
3
PV
1.0250
4.6648
0.8500
-0.1086
0
0
4
PQ
1.0258
-2.2168
0
0
0
0
5
PQ
0.9956
-3.9888
0
0
1.2500
0.5000
6
PQ
1.0127
-3.6874
0
0
0.9000
0.3000
7
PQ
1.0258
3.7197
0
0
0
0
8
PQ
1.0159
0.7275
0
0
1.0000
0.3500
9
PQ
1.0324
1.9667
0
0
0
0
3.2发电机初态
表3.2发电机初态
发电机
内电压
功角
1
1.0566
2.2716
2
1.0502
19.7316
3
1.0170
13.1664
3.3故障切除时间影响暂态稳定性仿真
假设0时刻发生故障，经过5工频周期后切除故障，即tc=0.083s，重合闸不投入。
得到三台发电机的转速和功角变化曲线。
图3.1发电机转速变化曲线
图3.2发电机功角变化曲线
图3.1和图3.2分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4s。
从曲线中可以看出系统没有失去暂态稳定。
现在改变故障切除时间tc为12个工频周期，即0.2s，三台发电机的动态曲线分别如下。
图3.3发电机转速变化曲线
图3.4发电机功角变化曲线
图3.3和图3.4分别为该过程中三台发电机的转速变化曲线和功角变化曲线，仿真总时间为4s。
从曲线中可以明显看出系统已经失去暂态稳定。
因此可以得出故障切除时间会影响系统的暂态稳定性，并且时间越短越有利于提高系统的暂态稳定性。
通过改变故障切除时间，可以求出临界故障切除时间tcc=0.16s，这时的发电机暂态曲线如下：
图3.5临界切除时间下的发电机转速变化曲线
图3.6临界切除时间下的发电机功角变化曲线
3.4
重合闸影响暂态稳定性仿真
图3.7tc=0.17s未投入重合闸时发电机转速变化曲线
图3.8tc=0.17s未投入重合闸时发电机功角变化曲线
取故障切除时间tc=0.17s，由3.3节可知，该时间大于临界故障切除时间，因此暂态稳定不能维持，发电机暂态曲线如图3.7和图3.8所示。
图3.9tc=0.17s投入重合闸时发电机转速变化曲线
图3.10tc=0.17s投入重合闸时发电机功角变化曲线
现将自动重合闸投入，重合闸动作延时为0.25s，即动作时刻为tar=0.42s，这时发电机的转速和功角变化曲线分别如图3.9和图3.10所示。
从图中课明显看出此时系统没有失去暂态稳定性，说明三相快速重合闸能够提高系统的暂态稳定性。
附录
MATLAB仿真程序代码：
%生成节点导纳矩阵
Ypofl=zeros（9）;
forn=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）,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）;
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）;
%计算潮流
whiletol>tolerance
%Delta_y
form=2:
9
sum=0;
forn=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
delta_y（m-1）=BUS（m,5）-BUS（m,7）-sum;
end
form=4:
9
sum=0;
forn=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
form=2:
9
forn=2:
9
ifm==n
sum=0;
fork=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
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
form=2:
9
forn=4:
9
ifm~=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
form=2:
7
sum=0;
forn=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,m）-sum;
end
%M
form=4:
9
forn=2:
9
ifm==n
sum=0;
fork=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,n）*sin（BUS（m,4）-BUS（n,4）））;
end
end
end
%L
form=4:
9
forn=4:
9
ifm==n
sum=0;
fork=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）*sin（BUS（m,4）-BUS（n,4））-B（m,n）*cos（BUS（m,4）-BUS（n,4）））;
end
end
end
%雅可比矩阵
JACOB=[HN;ML];
%修正量
delta_x=-JACOB\delta_y;
%修正电压幅值和相位
form=4:
9
BUS（m,3）=BUS（m,3）+delta_x（m+5）*BUS（m,3）;
end
form=2:
9
BUS（m,4）=BUS（m,4）+delta_x（m-1）;
end
%计算精度
tol=abs（max（delta_y））;
LoopCount=LoopCount+1;
end
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））;
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/（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_pf33Y_pf39;Y_pf93Y_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_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））-1j*BRANCH（3,5）;
Y_af99（5,7）=0;
Y_af99（7,5）=0;
Yraf=Y_af33-Y_af39/Y_af99*Y_af93;
disp（['Beforefault',date]）;
disp（Yrpf）;
disp（Yrdf）;
disp（Yraf）;
%故障前发电机状态
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））./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'）;
%故障中发电机状态
t_cut=0.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）.*180./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（:
4:
6）.*180./pi;
Omega_af=Xaf_out（:
1:
3）;
%重合闸后发电机状态
tspan=[t_ar,t_end];
x0_ar=Xaf_out（end,:
）;
[Tar_out,Xar_out]=ode45（@（t,x）Gen_Fun（t,x,GEN,Yrpf,E_pf）,tspan,x0_ar,options）;
Delta_ar=Xar_out（:
4:
6）.*180./pi;
Omega_ar=Xar_out（:
1:
3）;
%整合结果
T_out=[Tdf_out;Taf_out;Tar_out];
X_out=[Xdf_out;Xaf_out;Xar_out];
Delta_out=X_out（:
4:
6）.*180./pi;
Omega_out=X_out（:
1:
3）;
%绘制曲线
figure（'Name','deltacurve','NumberTitle','off'）;
plot（T_out,Delta_out（:
1）,'-',T_out,Delta_out（:
2）,'-.',T_out,Delta_out（:
3）,'--'）;
xlabel（'t（s）'）;
ylabel（'\delta_1,\delta_2,\delta_3（Deg）'）;
figure（'Name','omegacurve','NumberTitle','off'）;
plot（T_out,Omega_out）;
xl
升级会员 