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

三机九节点电力系统暂态仿真学院：自动化学院专业：电气工程专业学号：姓名：授课教师：江宁强一、摘要电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。本文采用IEEE3机9节点的经典多机模型，基于隐式梯形积分法对系统发生三相金属性短路故障进行仿真，分析系统在这种情况下的暂态稳定。发电机模型采用经典的二阶模型;负荷采用恒定阻抗负荷。在Matlab2010上编写程序进行调试和运行。电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。在动态稳定仿真中使用简单的电力系统模型，通过仿真计算分析说明，此仿真方法可以进行简单的电力系统暂态分析，对提高电力系统暂态稳定具有重要意义。二、案例本次课程主要应用P.M.AndersonandA.A.Fouad编写的《PowerSystemControlandStability》一书中所引用的WesternSystemCoordinatedCouncil(WSCC)三机九节点系统模型。系统电路结构拓扑图如下：0.0625j0.0586j0.0576j18/230230/13.816.5/23018KV230KV230KV230KV13.8KV16.5KV112233456789A负荷负荷B负荷C0.00850.072j/20.0745Bj/20.1045Bj/20.153Bj/20.179Bj/20.088Bj/20.079Bj0.01190.1008j0.0320.161j0.0100.085j0.0390.170j0.0170.092j图2-13机9节点系统系统数据其中，节点数据如下：节点号有无负载类型电压相角有功负荷无功负荷有功出力无功出力电压基准期望电压N=[1031.04000.000.000.0071.6027.0016.501.0402021.02500.000.000.00163.006.7018.001.0253021.02500.000.000.0085.00-10.9013.801.0254001.00000.000.000.000.000.00230.001.0265101.00000.00125.0050.000.000.000.000.9966101.00000.0090.0030.000.000.000.001.0137001.00000.000.000.000.000.00230.001.0268101.00000.00100.0035.000.000.000.001.0169001.00000.000.000.000.000.00230.001.032];%支路数据%从到电阻电抗容纳类型变比B=[140.00.05760.011270.00.06250.011390.00.05860.011450.0100.0850.17600460.0170.0920.15800570.0320.1610.30600690.0390.1700.35800780.00850.0720.14900890.01190.10080.20900];发电机数据如下：%发电机母线XdXd'Td0'XqXq'Tq0’TjXfGe=[110.14600.06088.960.09690.0969047.280.0576220.89580.11986.000.86450.19690.53512.800.0625331.31250.18138.591.25780.25000.6006.020.0585];三、仿真框图在仿真之前，首先，应明确仿真的所要到达的结果，即仿真目标:本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。故障前，系统处于稳定状态，发电机的攻角、转速基本稳定。而当系统发生故障以及故障切除，系统结构拓扑发生变化，系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的状态量为发电机攻角、发电机转速。其次，当明确仿真目标后，我们就得明确大体的仿真框架流程。仿真框架流程如下：数据准备(支路、节点、发电机)潮流计算计算故障前中后发电机内节点的导纳矩阵发电机初值计算列写系统状态方程(转子运动方程)调用ode45计算发电攻角、转速变化情况后处理图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是节点电压修正量。令ijijijiiijBGYjfeV;，则极坐标形式的功率不平衡量方程0)sincos(1njijijijijjiisiBGVVPP0)cossin(1njijijijijjiisiBGVVQQ雅可比矩阵J各元素的表达式LMNHJ当j≠i时：当j=i时：iijijiijijPHBeGffiiiiiiiiiiiPHBeGfbfiijijiijijPNGeBfeiiiiiiiiiiiPNGeBfaeiijijiijijQMBfGefiiiiiiiiiiiQMGeBfafiijijiijijQLGfBeeiiiiiiiiiiiQLBeGfbe其中，nijijjijinijijjijieBfGbfBeGa11)()(。进行牛顿拉夫逊算法得到潮流结果图4-3潮流结果二、故障前中后仅含发电机内节点的导纳矩阵图4-4故障前中后仅含发电机内节点的导纳矩阵三、求解电磁功率得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有Pe=real(E*I)如上中，E为发电机内电势，I为从发电机流出的电流。但在参考文献RamnarayanPatel,T.S.BhattiandD.P.Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式为：nijjjiijijjiiiieiYEEGEP12)cos(稳态情况下有，机械功率Pme=Pe四、求解运动方程发电机的运动方程可以写成常微分方程组：RinijjjiijijjiiiimijjiRidtdYEEGEPDdtdH12)cos(2其中Pmi为第i个机组故障前稳态的电磁功率。在本次仿真中Djωi为零，即阻尼为零。仿真开始，t=0时引入故障，0.083s后切除故障。求解运动方程后得到曲线如下：五、结果分析上图分别显示了各台发电机的转子角与时间的关系曲线，显示了发电机转速差的曲线，和2121、3131的曲线，由图可以看到，最大角差21为85o，出现在0.4ts处，无论是21还是31第二个摇摆都不大于第一个摇摆，可见系统是稳定的。六、程序代码主程序：globalPmYrungenGenE%节点数据%节点号有无负载类型电压相角有功负荷无功负荷有功负荷无功负荷电压基准期望电压N=[1031.04000.000.000.0071.6027.0016.501.0402021.02500.000.000.00163.006.7018.001.0253021.02500.000.000.0085.00-10.9013.801.0254001.00000.000.000.000.000.00230.001.0265101.00000.00125.0050.000.000.000.000.9966101.00000.0090.0030.000.000.000.001.0137001.00000.000.000.000.000.00230.001.0268101.00000.00100.0035.000.000.000.001.0169001.00000.000.000.000.000.00230.001.032];%支路数据%从到电阻电抗容纳类型变比B=[140.00.05760.011270.00.06250.011390.00.05860.011450.0100.0850.17600460.0170.0920.15800570.0320.1610.30600690.0390.1700.35800780.00850.0720.14900890.01190.10080.20900];%发电机数据%HMVAxd'*10000nodexdxqxlxadxaqxftd0'rfgen=[2364247.560810.14600.09690.03360.11240.06330.14838.960.0000439640192.0119820.89580.86450.05210.84370.81240.91736.000.0004054301128.0181331.31251.25780.07421.23831.28361.35555.890.0006105];Y=zeros(9,9);%导纳矩阵fori=1:9a=B(i,1);b=B(i,2);ifB(i,6)==0Y(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;elseY(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;endend%导纳矩阵forT=1:100[dP,dQ]=Caoliu(N,Y);%潮流J=Ykb(N,Y);%雅克比矩阵U=zeros(6);fori=4:9U(i-3,i-3)=N(i,4);enddAngU=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