电力系统仿真课程报告院系：自动化学院专业：电力系统及其自动化学号：************ *名：***时间：2015年1月6日1研究对象1.1三机九节点系统模型图1.1 WSCC-9系统模型图1.1是一个三机九节点的系统阻抗图，图中给出的阻抗参数都是以100MV A为基准的标幺值。

该图中包括三台发电机，三台双绕组变压器，九条母线（节点）和三个负荷。

本文将对该系统的动态过程进行相应的仿真分析。

1.2系统参数1.2.1节点参数按照节点类型，9个节点分为，给出已知参数如下表：上表中发电机有功、无功出力和负荷的有功无功功率均为以100MV A为基准时的标幺值。

1.2.2支路参数表1.2 支路参数上表中所有的参数均为标幺值，对于变压器支路。

最后三行表示三台变压器参数，已经计算出变压器的等效电抗并直接在表格中给出。

1.2.3发电机参数对于发电机，采用二阶经典模型，并对系统作出如下假设：(1)输入的机械功率保持恒定；(2)忽略阻尼效应；(3)负荷采用恒阻抗模型。

以上阻抗参数均以标幺值表示，额定转速下存储的能量（H）也转化为以100MV A 为基准的标幺值。

1.3需要求解的动态过程系统稳定运行时，0s时刻在线路5 – 7末端靠近母线7附近发生三相短路故障，故障持续5个周期（0.083s），仿真0s – 0.083s时的动态过程；然后切除5 – 7线路，仿真0.083s – 2s时的系统动态过程[1]。

2动态仿真过程2.1仿真总流程图2.1 仿真过程总流程2.2 潮流计算流程潮流计算采用牛顿——拉夫逊法，此方法的具有良好的收敛性，其解算过程理论依据如下。

极坐标形式的功率方程(cos sin ),(1,2,,n)(sin cos )i i j ij ij ij ij j ii i j ij ij ij j i P U U G B i Q U U G B θθθθ∈∈⎧=+⎪=⎨=-⎪⎩∑∑ (2.1)将其改写为增量的形式 有功功率：111111112222222211,11,1,1,1,1(cos sin )0(cos sin )0(cos sin )0s j j j j j j s j j j j j j n n s n j n j n j n j n j j n P P U U G B P P U U G B P P U U G B θθθθθθ∈∈-------∈-⎧∆=-+=⎪⎪∆=-+=⎪⎨⎪⎪∆=-+=⎪⎩∑∑∑ (2.2)无功功率：11,11,1,1,1,122,22,2,2,2,21n 1,1n r 1,n r 1,n r 1,n r 1,1(sin cos )0(sin cos )0(sin cos )0s j j j j j j s j j j j j j n r r s n r j j j j j j n r Q Q U U G B Q Q U U G B Q Q U U G B θθθθθθ∈∈--------------∈--⎧∆=--=⎪⎪∆=--=⎪⎨⎪⎪∆=--=⎪⎩∑∑∑ (2.3) 以上两式经泰勒级数展开可线性化为：111121,n 111121,1221222,121222,11,11,21,11,11,21,1111121,111121,1121222,12121n n n n n n n n n n n n n n n n r P H H H N N N P HH H N N N H H H N N N P M M M L L L Q MM M L Q Q ------------------∆⎡⎤⎢⎥∆⎢⎥⎢⎥⎢⎥∆⎢⎥⎢⎥=-⎢⎥∆⎢⎥⎢⎥∆⎢⎥⎢⎥⎢⎥∆⎣⎦1111222,1221,11,21,11,11,21,111///n n n r n r n r n r n r n r n r n r n r n r U U L L U U MM M L L L U U θθθ----------------------∆⎡⎤⎡⎤⎢⎥⎢⎥∆⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆⎢⎥⎢⎥⎢⎥⎢⎥∆⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆⎣⎦⎣⎦上式可简写为：⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ΔθΔP H N =-ΔU ΔQ M L U (2.4)进一步改写为： Δy =-J Δx (2.5)即-1Δx =-J Δy(2.6)其中J 为雅可比矩阵，在每次潮流计算结束后，通过式(2.6)得出修正量，根据修正量对x 的值进行修正，并再次计算潮流，直到精度达到给定的精度限制值。

图2.2 潮流子程序流程图通过潮流计算可以得到每一个节点的电压幅值和角度，得到每个节点的有功功率和无功功率，这些数据是进一步分析动态过程的基础。

2.3 计算发电机初态在动态分析时，首先应将发电机和负荷用相应的模型等效。

本文将发电机等效为二阶经典模型，将负荷等效为恒阻抗负荷。

图2.3 发电机等效模型本系统中具有三台发电机，因此会引入三个内节点。

因此节点导纳矩阵Y 将增广到12阶[2]。

可表示为：333912129399⨯⨯⨯⨯⨯⎡⎤=⎢⎥⎣⎦Y Y Y Y Y 在33⨯Y 中：'1,1,2,3j ii diY i x == (2.7)在39⨯Y 和93⨯Y 中：'1,1,2,3j ii diY i x =-= (2.8)在99⨯Y 中：'1,1,2,3j ii ii diY Y i x =+= (2.9)22j ,5,6,8i iii ii i i P Q Y Y i V V =+-= (2.10)将网络等效在发电机内节点构成的网络中，可以得到一个3阶的降阶节点导纳矩阵，该矩阵可由以下关系解出：133399993-⨯⨯⨯⨯=-Y Y Y Y Y(2.11)i iE δ∠I对于故障前和故障后的降阶节点导纳矩阵，利用式(2.11)可以计算出对应的降阶节点导纳矩阵；对于故障中的节点导纳矩阵，在12阶的增广矩阵中去掉故障母线所在的那一行和那一列，利用下式计算降阶矩阵。

133388883-⨯⨯⨯⨯=-Y Y Y Y Y(2.12)发电机初态中包含初始的电压幅值和功角，根据图 2.3，可以求解发电机1、2、3的初始状态。

*'jQ j(),1,2,3i ii i di iP E V x i V +=+=(2.13)i i i E E δ=∠(2.14)2.4 列写发电机动态方程发电机采用经典模型，其动态方程为：()2,1,2,3m e ri i i rP P d dt H i d dtωωδωω-⋅⎧=⎪⎪=⎨⎪=-⎪⎩ (2.15)三台发电机共有六个状态量，定义：123123x ωωωδδδ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦得到：(),x F t x =(2.16)2.5 求解发电机状态利用龙格——库塔数值积分方法，求解对发电机状态方程积分，求解发电机的状态。

认为在故障过程中，发电机内电压的幅值是不变的，只有功角改变；认为发电机的机械功率是不变的，其值等于故障前的发电机电磁功率。

3仿真程序运行结果与分析3.1潮流计算结果图3.1 潮流计算结果3.2降阶矩阵图3.2 降阶矩阵3.3发电机初态图3.3 发电机初态3.4动态过程3.4.10.083秒时刻的发电机状态图3.4 0.083s时的发电机状态3.4.22秒时刻的发电机状态图3.5 2s时刻的发电机状态3.4.3功角和转速的变化曲线图3.6 功角变化曲线图3.7 功角差变化曲线图3.8 转速变化曲线3.5结果分析对于三台发电机、九母线的系统，在0 – 2s的时间内，计算机给出了摇摆方程的数值积分。

图3.6给出了功角变化曲线，图3.7给出了功角差的变化曲线。

从中可以看出，这个系统是稳定的。

最大功角差为85°，由于第二次摇摆以及以后的摇摆都不可能大于第一次，那么第一次摇摆的特性就完全可以断定系统是否稳定。

如果功角差达到最大值然后就减小，则系统是稳定的。

如果功角差持续地增大，则系统是不稳定的，至少会有一台发电机失去同步。

参考文献[1] 福阿德, 安德森. 电力系统控制与稳定[M]. 电子工业出版社, 2012.[2] R. Patel, T. S. Bhatti. MATLAB/Simulink-based transient stability analysis of a multimachine power system[J]. International Journal of Electrical Engineering Education,4(39):320-336.附录clear;clc;%节点参数% 节点节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功BUS =[13 1.04000.71600.27050022 1.0250 1.63000.06650032 1.02500.8500-0.10860041 1.0000000051 1.000000 1.25000.500061 1.0000000.90000.300071 1.0000000081 1.000000 1.00000.350091 1.00000000];%支路参数% 首节点末节点电阻（p.u.）电抗（p.u.）电纳（p.u.）/2 变压器非标准变比BRANCH =[450.01000.08500.08801460.01700.09200.07901570.03200.16100.15301690.03900.17000.17901780.00850.07200.07451890.01190.10080.10451140.00000.05760.00001270.00000.06250.00001390.00000.05860.00001];%发电机参数% 节点 Xd X'd HGEN =[10.14600.060823.6420.89580.1198 6.403 1.31250.1813 3.01];%生成节点导纳矩阵Ypofl = zeros(9);for n =1:9Ypofl(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));endG = 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);%计算潮流while tol > tolerance% Delta_yfor m =2:9sum =0;for n =1:9sum = 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)));enddelta_y(m-1)= BUS(m,5)- BUS(m,7)- sum;endfor m =4:9sum =0;for n =1:9sum =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)));enddelta_y(m+5)= BUS(m,6)- BUS(m,8)- sum;end% Hfor m =2:9for n =2:9if m == nsum =0;for k =1:9sum =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)));endH(m-1,m-1)= BUS(m,3)^2* B(m,m)+ sum;elseH(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)));endendend% Nfor m =2:9for n =4:9if m ~= 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)));endendendfor m =2:7sum =0;for n =1:9sum = 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)));endN(m-1,m-1)=- BUS(m,3)^2* G(m,m)- sum;end% Mfor m =4:9for n =2:9if m == nsum =0;for k =1:9sum = 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)));endM(m-3,m-3)= BUS(m,3)^2* G(m,m)- sum;elseM(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)));endendend% Lfor m =4:9for n =4:9if m == nsum =0;for k =1:9sum = 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)));endL(m-3,m-3)= BUS(m,3)^2* B(m,m)- sum;elseL(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)));endendend% 雅可比矩阵JACOB =[H N;M L];% 修正量delta_x =- JACOB \ delta_y;% 修正电压幅值和相位for m =4:9BUS(m,3)= BUS(m,3)+ delta_x(m+5)* BUS(m,3);endfor m =2:9BUS(m,4)= BUS(m,4)+ delta_x(m-1);end% 计算精度tol = abs(max(delta_y));LoopCount = LoopCount +1;endBUS(:,4)= BUS(:,4)*180/pi;disp(['潮流计算结果 ',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_pf33 Y_pf39; Y_pf93 Y_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(['故障前、故障中和故障后的降阶矩阵 ',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(['发电机初态 ',date]);disp('内电压');disp(E_pf');disp('功角');disp(Delta_pf');% 故障中发电机状态t_cut =0.083;%故障切除时间t_end =2;%仿真结束时间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);disp('0.083s时的功角');disp(Delta_df(end,:));disp('0.083s时的转速');disp(Omega_df(end,:));% 故障后发电机状态tspan =[t_cut , t_end];x0_af = transpose(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);disp('2s时的功角');disp(Delta_af(end,:));disp('2s时的转速');disp(Omega_af(end,:));% 整合结果T_out =[Tdf_out ; Taf_out];X_out =[Xdf_out ; Xaf_out];Delta_out = X_out(:,4:6).*180./pi;Delta21_out = Delta_out(:,2)- Delta_out(:,1);Delta31_out = Delta_out(:,3)- Delta_out(:,1);Omega_out =[Omega_df ; Omega_af];% 绘制曲线figure('Name','功角变化曲线','NumberTitle','off');plot(T_out,Delta_out(:,1),'-',T_out,Delta_out(:,2),'-.',T_out,Delta_out(:,3 ),'--');。