matlab三机九节点电力系统仿真带程序.docx
《matlab三机九节点电力系统仿真带程序.docx》由会员分享，可在线阅读，更多相关《matlab三机九节点电力系统仿真带程序.docx（17页珍藏版）》请在冰豆网上搜索。
matlab三机九节点电力系统仿真带程序
三机九节点电力系统
暂态仿真
学院：
自动化学院
专业：
电气工程专业
学号：
姓名：
授课教师：
江宁强
一、摘要
电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。
通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。
本文采用IEEE3机9节点的经典多机模型，基于隐式梯形积分法对系统发生三相金属性短路故障进行仿真，分析系统在这种情况下的暂态稳定。
发电机模型采用经典的二阶模型;负荷采用恒定阻抗负荷。
在Matlab2010上编写程序进行调试和运行。
电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。
其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。
由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。
在动态稳定仿真中使用简单的电力系统模型，通过仿真计算分析说明，此仿真方法可以进行简单的电力系统暂态分析，对提高电力系统暂态稳定具有重要意义。
二、案例
本次课程主要应用P.M.AndersonandA.A.Fouad编写的《PowerSystemControlandStability》一书中所引用的WesternSystemCoordinatedCouncil（WSCC）三机九节点系统模型。
系统电路结构拓扑图如下：
图2-13机9节点系统
系统数据其中，节点数据如下：
节点号有无负载类型电压相角有功负荷无功负荷有功出力无功出力电压基准期望电压
N=[1031.04000.000.000.0071.6027.0016.501.040
2021.02500.000.000.00163.006.7018.001.0253021.02500.000.000.0085.00-10.9013.801.025
4001.00000.000.000.000.000.00230.001.026
5101.00000.00125.0050.000.000.000.000.996
6101.00000.0090.0030.000.000.000.001.013
7001.00000.000.000.000.000.00230.001.026
8101.00000.00100.0035.000.000.000.001.016
9001.00000.000.000.000.000.00230.001.032];
%支路数据
%从到电阻电抗容纳类型变比
B=[140.00.05760.011
270.00.06250.011
390.00.05860.011
450.0100.0850.17600
460.0170.0920.15800
570.0320.1610.30600
690.0390.1700.35800
780.00850.0720.14900
890.01190.10080.20900];
发电机数据如下：
%发电机母线XdXd'Td0'XqXq'Tq0’TjXf
Ge=[110.14600.06088.960.09690.0969047.280.0576
220.89580.11986.000.86450.19690.53512.800.0625
331.31250.18138.591.25780.25000.6006.020.0585];
三、仿真框图
在仿真之前，首先，应明确仿真的所要到达的结果，即仿真目标:
本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。
故障前，系统处于稳定状态，发电机的攻角、转速基本稳定。
而当系统发生故障以及故障切除，系统结构拓扑发生变化，系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。
此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的
状态量为发电机攻角、发电机转速。
其次，当明确仿真目标后，我们就得明确大体的仿真框架流程。
仿真框架流程如下：
图3-1仿真流程图
四、仿真模型
在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。
一、潮流计算
由于本文以三机九节点为模型，假定节点一为参考节点，这样就有2两个发电机的PV节点，6个PQ节点，未知量为8个节点（包括2个PV节点和6个PQ节点）的电压相角，还有6个节点（PQ节点）的电压幅值。
可以先求出Y矩阵
图4-1Y矩阵
然后，我们列写方程，也就是利用各个节点的有功、无功功率的平衡关系，列写14个功率平衡方程。
这样就能使用牛顿一拉夫逊算法来求解这14个非线性方程。
其中的关键是要计算出雅克比矩阵
图4-2雅克比矩阵
然后计算出修正量。
在设定精度和最大迭代次数的前提下进行迭代，直到满足要求。
电力网络的节点功率方程可用一般形式表示如下:
牛顿拉夫逊算法修正方程
W=-JΔV
其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵；ΔV是节点电压修正量。
令
，
则极坐标形式的功率不平衡量方程
雅可比矩阵J各元素的表达式
当j≠i时：
当j=i时：
其中，
。
进行牛顿拉夫逊算法得到潮流结果
图4-3潮流结果
二、故障前中后仅含发电机内节点的导纳矩阵
图4-4故障前中后仅含发电机内节点的导纳矩阵
三、求解电磁功率
得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。
所以可以有
Pe=real（E*I）
如上中，E为发电机内电势，I为从发电机流出的电流。
但在参考文献RamnarayanPatel,T.S.BhattiandD.P.Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式为：
稳态情况下有，机械功率Pme=Pe
四、求解运动方程
发电机的运动方程可以写成常微分方程组：
其中Pmi为第i个机组故障前稳态的电磁功率。
在本次仿真中Djωi为零，即阻尼为零。
仿真开始，t=0时引入故障，0.083s后切除故障。
求解运动方程后得到曲线如下：
五、结果分析
上图分别显示了各台发电机的转子角与时间的关系曲线，显示了发电机转速差的曲线，和
、
的曲线，由图可以看到，最大角差
为
，出现在
处，无论是
还是
第二个摇摆都不大于第一个摇摆，可见系统是稳定的。
六、程序代码
主程序：
globalPmYrungenGenE
%节点数据
%节点号有无负载类型电压相角有功负荷无功负荷有功负荷无功负荷电压基准期望电压
N=[1031.04000.000.000.0071.6027.0016.501.040
2021.02500.000.000.00163.006.7018.001.025
3021.02500.000.000.0085.00-10.9013.801.025
4001.00000.000.000.000.000.00230.001.026
5101.00000.00125.0050.000.000.000.000.996
6101.00000.0090.0030.000.000.000.001.013
7001.00000.000.000.000.000.00230.001.026
8101.00000.00100.0035.000.000.000.001.016
9001.00000.000.000.000.000.00230.001.032];
%支路数据
%从到电阻电抗容纳类型变比
B=[140.00.05760.011
270.00.06250.011
390.00.05860.011
450.0100.0850.17600
460.0170.0920.15800
570.0320.1610.30600
690.0390.1700.35800
780.00850.0720.14900
890.01190.10080.20900];
%发电机数据
%HMVAxd'*10000nodexdxqxlxadxaqxftd0'rf
gen=[2364247.560810.14600.09690.03360.11240.06330.14838.960.0000439
640192.0119820.89580.86450.05210.84370.81240.91736.000.0004054
301128.0181331.31251.25780.07421.23831.28361.35555.890.0006105];
Y=zeros（9,9）;%导纳矩阵
fori=1:
9
a=B（i,1）;b=B（i,2）;
ifB（i,6）==0
Y（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）+B（i,4）*j）+B（i,5）*1i./2;
else
Y（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）^2）+B（i,5）*1i./2;
end
end%导纳矩阵
forT=1:
100
[dP,dQ]=Caoliu（N,Y）;%潮流
J=Ykb（N,Y）;%雅克比矩阵
U=zeros（6）;
fori=4:
9
U（i-3,i-3）=N（i,4）;
end
dAngU=J\[dP;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）
break
end
end
[Yc,Yb,Ya]=Ynew（gen,N,B,Y）;
GEgj=zeros（1,3）;
GenE=zeros（1,3）;
fori=1:
3
GenE（1,i）=abs（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））））））;
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））））））;
end
Yrun=zeros（3）;
Pm=zeros（1,3）;
fori=1:
3
Pm（1,i）=N（i,8）./100;
end
options=odeset（'RElTOL',1e-10）;%设置精度
X0=[GEgj
（1）,GEgj
（2）,GEgj（3）,2*pi*60*ones（1,3）];
t0=0;
tc=0.083;
tspan1=[t0,tc];
Yrun=Yb;
[T1,Y1]=ode45（'fzz',tspan1,X0,options）;
X1=Y1（end,:
）;
tf=2;
tspan2=[tc,tf];
Yrun=Ya;
[T2,Y2]=ode45（'fzz',tspan2,X1,options）;
T=[T1;T2];
Y12=[Y1;Y2];
subplot（3,1,1）;
plot（T,Y12（:
1）*180/pi,T,Y12（:
2）*180/pi,T,Y12（:
3）*180/pi）;%发电机功角
subplot（3,1,2）;
plot（T,Y12（:
5）-Y12（:
4）,T,Y12（:
6）-Y12（:
4））;%发电机转速差
subplot（3,1,3）;
plot（T,Y12（:
2）*180/pi-Y12（:
1）*180/pi,T,Y12（:
3）*180/pi-Y12（:
1）*180/pi）%发电机攻角差
雅克比矩阵：
functionJ=Ykb（N,Y）
H=zeros（8）;N1=zeros（8,6）;K=zeros（6,8）;L=zeros（6）;
fori=2:
9
forj=2:
9
ifi~=j
H（i-1,j-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）））;
else
H（i-1,j-1）=（N（i,9）-N（i,7））/100+imag（Y（i,j））*（（N（i,4））^2）;
end
end
end%H
fori=2:
9
forj=4:
9
ifi~=j
N1（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）））;
else
N1（i-1,j-3）=-（N（i,8）-N（i,6））/100-real（Y（i,j））*（（N（i,4））^2）;
end
end
end%N
fori=4:
9
forj=2:
9
ifi~=j
K（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）））;
else
K（i-3,j-1）=-（N（i,8）-N（i,6））/100+real（Y（i,j））*（（N（i,4））^2）;
end
end
end%K
fori=4:
9
forj=4:
9
ifi~=j
L（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（N（i,5）-N（j,5）））;
else
L（i-3,j-3）=-（N（i,9）-N（i,7））/100+imag（Y（i,j））*（（N（i,4））^2）;
end
end
end%L
J=[HN1;KL];
潮流计算
function[dP,dQ]=Caoliu（N,Y）
dP=zeros（8,1）;
dQ=zeros（6,1）;
fori=2:
9
dP（i-1,1）=（N（i,8）-N（i,6））/100;
end
fori=4:
9
dQ（i-3,1）=（N（i,9）-N（i,7））/100;
end
fori=2:
9
forj=1:
9
dP（i-1,1）=dP（i-1,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）））;
end
end
fori=4:
9
forj=1:
9
dQ（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）））;
end
end
故障前中后仅含发电机内节点的导纳矩阵：
function[Yp,Yd,Ya]=Ysim（gen,N,B,Y）
Yp=zeros（3）;Yd=zeros（3）;Ya=zeros（3）;
Y1=zeros（12）;Y2=zeros（11）;Y3=zeros（12）;
Y1=[zeros（3）,zeros（3,9）;zeros（9,3）,Y];%故障前增广
Y5=N（5,6）/100/（（N（5,4）^2））-（N（5,7）/100/（（N（5,4）^2）））*1i;
Y6=N（6,6）/100/（（N（6,4）^2））-（N（6,7）/100/（（N（6,4）^2）））*1i;
Y8=N（8,6）/100/（（N（8,4）^2））-（N（8,7）/100/（（N（8,4）^2）））*1i;%负载等效导纳
fori=1:
3
Y1（i,i）=-（1/gen（i,3）*10000）*1i;
Y1（i,i+3）=（1/gen（i,3）*10000）*1i;
Y1（i+3,i）=Y1（i,i+3）;
end
fori=1:
3
Y1（i+3,i+3）=Y1（i+3,i+3）-（1/gen（i,3）*10000）*1i;%发电机节点修改
end
Y1（3+5,3+5）=Y1（3+5,3+5）+Y5;
Y1（3+6,3+6）=Y1（3+6,3+6）+Y6;
Y1（3+8,3+8）=Y1（3+8,3+8）+Y8;%负载节点修改
Ynn=zeros（3）;Ynr=zeros（3,9）;Ynr1=zeros（3,8）;Yrn=zeros（9,3）;Yrn1=zeros（8,3）;
Yrr=zeros（9）;Yrr1=zeros（8）;
Ynn=Y1（1:
3,1:
3）;Ynr=Y1（1:
3,4:
12）;Yrn=Y1（4:
12,1:
3）;Yrr=Y1（4:
12,4:
12）;
Yp=Ynn-Ynr*（Yrr^-1）*Yrn;
Y3=Y1;
Y3（8,8）=Y3（8,8）-（1./（B（6,3）+B（6,4）*1i）+B（6,5）*1i/2）;Y3（10,10）=Y3（10,10）-（1./（B（6,3）+B（6,4）*1i）+B（6,5）*1i/2）;Y3（8,10）=0;Y3（10,8）=0;%故障后增广
Y1（10,:
）=[];
Y1（:
10）=[];
Y2=Y1;%故障中增广
Ynn=Y2（1:
3,1:
3）;Ynr1=Y2（1:
3,4:
11）;Yrn1=Y2（4:
11,1:
3）;Yrr1=Y2（4:
11,4:
11）;
Yd=Ynn-Ynr1*（Yrr1^-1）*Yrn1;
Ynn=Y3（1:
3,1:
3）;Ynr=Y3（1:
3,4:
12）;Yrn=Y3（4:
12,1:
3）;Yrr=Y3（4:
12,4:
12）;
Ya=Ynn-Ynr*（Yrr^-1）*Yrn;
求解运动方程：
functionxtt=fouad3（T,DX）
globalPmYrungenGenE
xtt=[DX（4）-2*pi*60;DX（5）-2*pi*60;DX（6）-2*pi*60;
（Pm
（1）-GenE
（1）*（real（Yrun（1,1））*GenE
（1）+real（Yrun（1,2））*GenE
（2）*cos（DX
（1）-DX
（2））+real（Yrun（1,3））*GenE（3）*cos（DX
（1）-DX（3））...
+imag（Yrun（1,2））*GenE
（2）*sin（DX
（1）-DX
（2））+imag（Yrun（1,3））*GenE（3）*sin（DX
（1）-DX（3））））/2/gen（1,1）*100*2*pi*60;
（Pm
（2）-GenE
（2）*（real（Yrun（2,2））*GenE
（2）+real（Yrun（2,1））*GenE
（1）*cos（DX
（2）-DX
（1））+real（Yrun（2,3））*GenE（3）*cos（DX
（2）-DX（3））...
+imag（Yrun（2,1））*GenE
（1）*sin（DX
（2）-DX
（1））+imag（Yrun（2,3））*GenE（3）*sin（DX
（2）-DX（3））））/2/gen（2,1）*100*2*pi*60;
（Pm（3）-GenE（3）*（real（Yrun（3,3））*GenE（3）+real（Yrun（3,1））*GenE
（1）*cos（DX（3）-DX
（1））+real（Yrun（3,2））*GenE
（2）*cos（DX（3）-DX
（2））...
+imag（Yrun（3,1））*GenE
（1）*sin（DX（3）-DX
（1））+imag（Yrun（3,2））*GenE
（2）*sin（DX（3）-DX
（2））））/2/gen（3,1）*100*2*pi*60];
七、总结
首先感谢任课老师江宁强老师对我仿真的指导与帮助，还要感谢在程序编写过程中给过我无私帮助的其他同学。
这次课程的仿真设计实验让我学到了很多，本文中切断时间在0.083s，最终的图形和老师的有一定的误差，这是我今后要完善的地方。
刚开始y举证先是求的不准确,少了对地导纳。
之后求解雅克比矩阵分块求解不是太准确。
对于电磁功率如果按照文献RamnarayanPatel,T.S.BhattiandD.P.Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，简化了电磁功率的计算，这样一来计算量就少了很多。
升级会员 