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

(完整版)matlab三机九节点电力系统仿真(带程序)三机九节点电力系统暂态仿真学院：自动化学院专业：电气工程专业学号：姓名：授课教师：江宁强一、摘要 电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。本文采用IEEE 3机9节点的经典多机模型，基于隐式梯形积分法对系统发生三相金属性短路故障进行仿真，分析系统在这种情况下的暂态稳定。发电机模型采用经典的二阶模型;负荷采用恒定阻抗负荷。在Matlab2010上编写程序进行调试和运行。电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统.其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组.在动态稳定仿真中使用简单的电力系统模型,通过仿真计算分析说明，此仿真方法可以进行简单的电力系统暂态分析，对提高电力系统暂态稳定具有重要意义。二、案例本次课程主要应用P. M。 Anderson and A. A。 Fouad编写的Power System Control and Stability一书中所引用的Western System Coordinated Council （WSCC）三机九节点系统模型。系统电路结构拓扑图如下：图2-1 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。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 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。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；三、仿真框图 在仿真之前，首先，应明确仿真的所要到达的结果,即仿真目标：本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。故障前，系统处于稳定状态，发电机的攻角、转速基本稳定.而当系统发生故障以及故障切除，系统结构拓扑发生变化,系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的状态量为发电机攻角、发电机转速。其次，当明确仿真目标后，我们就得明确大体的仿真框架流程.仿真框架流程如下： 数据准备(支路、节点、发电机)后处理调用ode45计算发电攻角、 转速变化情况列写系统状态方程 （转子运动方程）计算故障前中后发电机 内节点的导纳矩阵 发电机初值计算 潮流计算图31仿真流程图四、 仿真模型在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。一、 潮流计算 由于本文以三机九节点为模型，假定节点一为参考节点，这样就有2两个发电机的PV节点，6个PQ节点，未知量为8个节点（包括2个PV节点和6个PQ节点)的电压相角，还有6个节点（PQ节点）的电压幅值。可以先求出Y矩阵图41 Y矩阵 然后,我们列写方程,也就是利用各个节点的有功、无功功率的平衡关系，列写14个功率平衡方程。这样就能使用牛顿一拉夫逊算法来求解这14个非线性方程.其中的关键是要计算出雅克比矩阵图4-2雅克比矩阵然后计算出修正量.在设定精度和最大迭代次数的前提下进行迭代,直到满足要求。 电力网络的节点功率方程可用一般形式表示如下：牛顿拉夫逊算法修正方程W = -JV其中W是节点不平衡量向量，包括有功,无功,电压；J是雅克比矩阵；V是节点电压修正量.令，则极坐标形式的功率不平衡量方程雅可比矩阵J各元素的表达式当ji时: 当j=i时： 其中，。进行牛顿拉夫逊算法得到潮流结果图43潮流结果二、故障前中后仅含发电机内节点的导纳矩阵图44故障前中后仅含发电机内节点的导纳矩阵三、 求解电磁功率得到故障前，故障中，故障后三个不同的导纳矩阵后,就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有Pe=real(EI） 如上中，E为发电机内电势，I为从发电机流出的电流.但在参考文献Ramnarayan Patel， T. S。 Bhatti and D. P。 Kothari。MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式为：稳态情况下有，机械功率Pme=Pe四、求解运动方程发电机的运动方程可以写成常微分方程组：其中Pmi为第i个机组故障前稳态的电磁功率.在本次仿真中Dji为零，即阻尼为零.仿真开始，t=0时引入故障,0。083s后切除故障.求解运动方程后得到曲线如下：五、 结果分析上图分别显示了各台发电机的转子角与时间的关系曲线，显示了发电机转速差的曲线，和、的曲线,由图可以看到，最大角差为，出现在处，无论是还是第二个摇摆都不大于第一个摇摆,可见系统是稳定的。六、 程序代码主程序：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.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=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 608 1 0.1460 0。0969 0。0336 0。1124 0.0633 0.1483 8.96 0。0000439 640 192。0 1198 2 0。8958 0.8645 0.0521 0.8437 0。8124 0.9173 6.00 0.0004054 301 128。0 1813 3 1。3125 1。2578 0.0742 1.2383 1。2836 1。3555 5.89 0。0006105;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(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)+1iB(i,4)）B(i，6); Y（b,a）=1。/（(B（i,3)+1iB(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）jB（i，6)2）+B（i,5）1i./2； end end 导纳矩阵for T=1：100 dP,dQ=Caoliu(N,Y);潮流 J=Ykb（N，Y）;雅克比矩阵 U=zeros(6）; for i=4:9 U（i-3，i3）=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(1，3）；for i=1:3 GenE(1,i)=abs（N（i,4)*exp（1i*N(i,5）+1igen(i,3）/10000conj(（N(i，8）/100+1iN（i,9)/100)/(N（i，4）*exp（1i*N（i,5）)）)）; GEgj(1，i)=angle(N（i，4)exp（1iN(i，5)）+1i*gen（i,3)/10000*conj(（N（i,8)/100+1i*N(i，9）/100)/（N（i，4)*exp（1iN(i，5）)）；endYrun=zeros（3）; Pm=zeros(1,3);for i=1:3 Pm（1，i)=N（i，8)。/100;endoptions=odeset（RElTOL，1e-10)； 设置精度X0=GEgj（1),GEgj(2），GEgj(3），2*pi60*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) 发电机攻角差 雅克比矩阵：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,j1）=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，j1)=（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（i1,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，j3)=（N(i，8）-N（i，6）/100real(Y(i，j)）*（(N（i，4）2）; end endend %Nfor i=4:9 for j=2：9 if i=j K（i3,j1)=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(i3，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,j3）=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 LJ=H N1;K L；潮流计算function dP，dQ=Caoliu(N，Y）dP=zeros（8，1）;dQ=zeros（6，1)；for i=2:9 dP（i1，1)=(N（i，8)N(i,6)/100;endfor i=4：9 dQ（i3,1）=(N（i，9)-N(i，7)/100;endfor i=2:9 for j =1:9 dP(i-1，1)=dP(i1，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(i3，1）=dQ(i3，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)；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; 负载等效导纳for i=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）;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（1:3，1:3）;Ynr=Y1(1：3，4:12)；Yrn=Y1(4：12，1：3)；Yrr=Y1(4：12,4:12)；Yp=YnnYnr*（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*(Yrr11)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=YnnYnr(Yrr-1）*Yrn;求解运动方程:function xtt=fouad3(T,DX) global Pm Yrun gen GenExtt=DX（4）2pi*60;DX（5)2*pi*60;DX(6）2pi*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*2pi60; （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*2pi*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)1002*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中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，简化了电磁功率的计算，这样一来计算量就少了很多。