matlab三机九节点电力系统仿真(带程序)(共19页).doc

《matlab三机九节点电力系统仿真(带程序)(共19页).doc》由会员分享，可在线阅读，更多相关《matlab三机九节点电力系统仿真(带程序)(共19页).doc（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上三机九节点电力系统暂态仿真学院：自动化学院专业：电气工程专业学号：姓名：授课教师：江宁强专心-专注-专业一、摘要 电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。本文采用IEEE 3机9节点的经典多机模型，基于隐式梯形积分法对系统发生三相金属性短路故障进行仿真，分析系统在这种情况下的暂态稳定。发电机模型采用经典的二阶模型;负荷采用恒定阻抗负荷。在Matlab2010上编写程序进行调试和运行。电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞

2、大复杂的动力学系统。其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。在动态稳定仿真中使用简单的电力系统模型，通过仿真计算分析说明，此仿真方法可以进行简单的电力系统暂态分析，对提高电力系统暂态稳定具有重要意义。二、案例本次课程主要应用P. M. Anderson and A. A. Fouad编写的Power System Control and Stability一书中所引用的Western System Coordinated Council (WSCC)三机九节点系统模型。系统电路结构拓扑图如下：图2-1 3

3、机9节点系统系统数据其中，节点数据如下：节点号 有无负载 类型 电压 相角 有功负荷 无功负荷 有功出力 无功出力 电压基准 期望电压 N=1 0 3 1.0400 0.00 0.00 0.00 71.60 27.00 16.50 1.040 2 0 2 1.0250 0.00 0.00 0.00 163.00 6.70 18.00 1.025 3 0 2 1.0250 0.00 0.00 0.00 85.00 -10.90 13.80 1.025 4 0 0 1.0000 0.00 0.00 0.00 0.00 0.00 230.00 1.026 5 1 0 1.0000 0.00 125.

4、00 50.00 0.00 0.00 0.00 0.996 6 1 0 1.0000 0.00 90.00 30.00 0.00 0.00 0.00 1.013 7 0 0 1.0000 0.00 0.00 0.00 0.00 0.00 230.00 1.026 8 1 0 1.0000 0.00 100.00 35.00 0.00 0.00 0.00 1.016 9 0 0 1.0000 0.00 0.00 0.00 0.00 0.00 230.00 1.032;%支路数据% 从 到 电阻 电抗 容纳 类型 变比B=1 4 0.0 0.0576 0.0 1 1 2 7 0.0 0.0625

5、0.0 1 1 3 9 0.0 0.0586 0.0 1 1 4 5 0.010 0.085 0.176 0 0 4 6 0.017 0.092 0.158 0 0 5 7 0.032 0.161 0.306 0 0 6 9 0.039 0.170 0.358 0 0 7 8 0.0085 0.072 0.149 0 0 8 9 0.0119 0.1008 0.209 0 0;发电机数据如下：% 发电机 母线 Xd Xd Td0 Xq Xq Tq0 Tj XfGe= 1 1 0.1460 0.0608 8.96 0.0969 0.0969 0 47.28 0.0576 2 2 0.8958 0

6、.1198 6.00 0.8645 0.1969 0.535 12.80 0.0625 3 3 1.3125 0.1813 8.59 1.2578 0.2500 0.600 6.02 0.0585;三、仿真框图 在仿真之前，首先，应明确仿真的所要到达的结果，即仿真目标:本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。故障前，系统处于稳定状态，发电机的攻角、转速基本稳定。而当系统发生故障以及故障切除，系统结构拓扑发生变化，系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到

7、系统状态的随时间的变化值。此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的状态量为发电机攻角、发电机转速。其次，当明确仿真目标后，我们就得明确大体的仿真框架流程。仿真框架流程如下： 数据准备(支路、节点、发电机)后处理调用ode45计算发电攻角、 转速变化情况列写系统状态方程 (转子运动方程)计算故障前中后发电机 内节点的导纳矩阵 发电机初值计算 潮流计算图3-1仿真流程图四、 仿真模型在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。一、 潮流计算 由于本文以三

8、机九节点为模型，假定节点一为参考节点，这样就有2两个发电机的PV节点，6个PQ节点，未知量为8个节点(包括2个PV节点和6个PQ节点)的电压相角，还有6个节点(PQ节点)的电压幅值。可以先求出Y矩阵图4-1 Y矩阵 然后，我们列写方程，也就是利用各个节点的有功、无功功率的平衡关系，列写14个功率平衡方程。这样就能使用牛顿一拉夫逊算法来求解这14个非线性方程。其中的关键是要计算出雅克比矩阵图4-2雅克比矩阵然后计算出修正量。在设定精度和最大迭代次数的前提下进行迭代，直到满足要求。 电力网络的节点功率方程可用一般形式表示如下:牛顿拉夫逊算法修正方程W = -JV其中W是节点不平衡量向量，包括有功，

9、无功，电压；J是雅克比矩阵；V是节点电压修正量。令，则极坐标形式的功率不平衡量方程雅可比矩阵J各元素的表达式当ji时： 当j=i时： 其中，。进行牛顿拉夫逊算法得到潮流结果图4-3潮流结果二、故障前中后仅含发电机内节点的导纳矩阵图4-4故障前中后仅含发电机内节点的导纳矩阵三、 求解电磁功率得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有Pe=real(E*I) 如上中，E为发电机内电势，I为从发电机流出的电流。但在参考文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothar

10、i.MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式为：稳态情况下有，机械功率Pme=Pe四、求解运动方程发电机的运动方程可以写成常微分方程组：其中Pmi为第i个机组故障前稳态的电磁功率。在本次仿真中Dji为零，即阻尼为零。仿真开始，t=0时引入故障，0.083s后切除故障。求解运动方程后得到曲线如下：五、 结果分析上图分别显示了各台发电机的转子角与时间的关系曲线，显示了发电机转速差的曲线，和、的曲线，由图可以看到，最大角差为，出现在处，无论是还是第二个摇摆

11、都不大于第一个摇摆，可见系统是稳定的。六、 程序代码主程序：global Pm Yrun gen GenE%节点数据%节点号 有无负载 类型 电压 相角 有功负荷 无功负荷 有功负荷 无功负荷 电压基准 期望电压 N=1 0 3 1.0400 0.00 0.00 0.00 71.60 27.00 16.50 1.040 2 0 2 1.0250 0.00 0.00 0.00 163.00 6.70 18.00 1.025 3 0 2 1.0250 0.00 0.00 0.00 85.00 -10.90 13.80 1.025 4 0 0 1.0000 0.00 0.00 0.00 0.00 0

12、.00 230.00 1.026 5 1 0 1.0000 0.00 125.00 50.00 0.00 0.00 0.00 0.996 6 1 0 1.0000 0.00 90.00 30.00 0.00 0.00 0.00 1.013 7 0 0 1.0000 0.00 0.00 0.00 0.00 0.00 230.00 1.026 8 1 0 1.0000 0.00 100.00 35.00 0.00 0.00 0.00 1.016 9 0 0 1.0000 0.00 0.00 0.00 0.00 0.00 230.00 1.032;%支路数据% 从 到 电阻 电抗 容纳 类型 变比B

13、=1 4 0.0 0.0576 0.0 1 1 2 7 0.0 0.0625 0.0 1 1 3 9 0.0 0.0586 0.0 1 1 4 5 0.010 0.085 0.176 0 0 4 6 0.017 0.092 0.158 0 0 5 7 0.032 0.161 0.306 0 0 6 9 0.039 0.170 0.358 0 0 7 8 0.0085 0.072 0.149 0 0 8 9 0.0119 0.1008 0.209 0 0;%发电机数据% H MVA xd*10000 node xd xq xl xad xaq xf td0 rfgen=2364 247.5 60

14、8 1 0.1460 0.0969 0.0336 0.1124 0.0633 0.1483 8.96 0. 640 192.0 1198 2 0.8958 0.8645 0.0521 0.8437 0.8124 0.9173 6.00 0. 301 128.0 1813 3 1.3125 1.2578 0.0742 1.2383 1.2836 1.3555 5.89 0.;Y=zeros(9,9);%导纳矩阵 for i=1:9 a=B(i,1);b=B(i,2); if B(i,6)=0 Y(a,b)=-1./(B(i,3)+B(i,4)*1i); Y(b,a)=-1./(B(i,3)+B(

15、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 %导纳矩阵for T=

16、1:100 dP,dQ=Caoliu(N,Y);%潮流 J=Ykb(N,Y);%雅克比矩阵 U=zeros(6); for i=4:9 U(i-3,i-3)=N(i,4); end dAngU=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) break endendYc,Yb,Ya=Ynew(gen,N,B,Y);GEgj=zeros(1,3);GenE=zeros

17、(1,3);for i=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);endYrun=zeros(3); Pm=zeros(1,3);for i=1:3 Pm(1,i)=N(i,8)./100;endopti

18、ons=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(

19、:,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) %发电机攻角差 雅克比矩阵：function J=Ykb(N,Y)H=zeros(8);N1=zeros(8,6);K=zeros(6,8);L=zeros(6);for i=2:9 for j=2:9 if i=j H(i-1,j-1)=-

20、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 endend %Hfor i=2:9 for j=4:9 if i=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-

21、real(Y(i,j)*(N(i,4)2); end endend %Nfor i=4:9 for j=2:9 if i=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 endend %Kfor i=4:9 for j=4:9 if i=j L(i-3,j-3)=-N(i,4)*N(j,4)*(real(Y(i,j)*sin(N(i,5)-N

22、(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 %LJ=H N1;K L;潮流计算function dP,dQ=Caoliu(N,Y)dP=zeros(8,1);dQ=zeros(6,1);for i=2:9 dP(i-1,1)=(N(i,8)-N(i,6)/100;endfor i=4:9 dQ(i-3,1)=(N(i,9)-N(i,7)/100;endfor i=2:9 for j =1:9 dP(i-1,1)=dP(i-1

23、,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); endend for i=4:9 for j =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); endend 故障前中后仅含发电机内节点的导纳矩阵：function Yp,Yd,Ya=Ysim(gen,N,B,Y)Yp=zeros(3);Yd=zeros(3);Ya=zeros(3);Y

24、1=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; %负载等效导纳for i=1:3 Y1(i,i)= -(1/gen(i,3)*10000)*1i; Y1(i,i+3)=(1/gen(i,3)

25、*10000)*1i; Y1(i+3,i)=Y1(i,i+3);endfor i=1:3 Y1(i+3,i+3)=Y1(i+3,i+3)-(1/gen(i,3)*10000)*1i; %发电机节点修改endY1(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

26、(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=

27、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;求解运动方程：function xtt=fouad3(T,DX) global Pm Yrun gen GenExtt=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

28、(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(Yru

29、n(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

30、)*100*2*pi*60;七、 总结首先感谢任课老师江宁强老师对我仿真的指导与帮助，还要感谢在程序编写过程中给过我无私帮助的其他同学。这次课程的仿真设计实验让我学到了很多，本文中切断时间在0.083s，最终的图形和老师的有一定的误差，这是我今后要完善的地方。刚开始y举证先是求的不准确,少了对地导纳。之后求解雅克比矩阵分块求解不是太准确。对于电磁功率如果按照文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，简化了电磁功率的计算，这样一来计算量就少了很多。