2020三机九节点电力系统建模与仿真

学院
专业
姓名
学号
指导教师
邮箱
提交日期
一、摘要
电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。

通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。

本文采用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.30
6 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 Xf
Ge=[ 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];
三、仿真框图
在仿真之前，首先，应明确仿真的所要到达的结果，即仿真目标:本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。

故障前，系统处于稳定状态，发电机的攻角、转速基本稳定。

而当系统发生故障以及故障切除，系统结构拓扑发生变化，系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。

此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的状态量为发电机攻角、发电机转速。

其次，当明确仿真目标后，我们就得明确大体的仿真框架流程。

仿真框架流程如下：
图3-1 仿真流程图
四、仿真模型
在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。

一、潮流计算
由于本文以三机九节点为模型，假定节点一为参考节点，这样就有2 两个发电机的PV 节点，6 个PQ 节点，未知量为8 个节点(包括2 个PV 节点和6 个PQ 节点)的电压相角，还有6 个节点(PQ 节点)的电压幅值。

可以先求出Y 矩阵
图4-1 Y 矩阵
然后，我们列写方程，也就是利用各个节点的有功、无功功率的平衡关系，列写14 个功率平衡方程。

这样就能使用牛顿一拉夫逊算法来求解这14 个非线性方程。

其中的关键是要计算出雅克比矩阵
⎨
L ⎬
图4-2 雅克比矩阵
然后计算出修正量。

在设定精度和最大迭代次数的前提下进行迭代，直到满足要求。

电力网络的节点功率方程可用一般形式表示如下:
牛顿拉夫逊算法修正方程
W = -JΔV
其中W 是节点不平衡量向量，包括有功，无功，电压；J 是雅克比矩阵；ΔV 是节点电压修正量。

令
V i =e
i
+jf
i
;Y
ij
=G
ij
+jB
ij ，
则极坐标形式的功率不平衡量方程
n
∆P
i =P
is
-V
i ∑V j (G ij cos δij +B ij sin δij ) = 0
j =1
n
∆Q
i =Q
is
-V
i ∑V j (G ij sin δij -B ij cosδij ) = 0
j =1
雅可比矩阵J 各元素的表达式
J =⎧H N ⎫ ⎩M ⎭
当j≠i 时：当j=i 时：
其中，
n
a i = ∑ (G ij e j - B ij f j )
i =1 n
b i = ∑ (G ij f j + B ij e j )
i =1。

进行牛顿拉夫逊算法得到潮流结果
图 4-3 潮流结果
二、故障前中后仅含发电机内节点的导纳矩阵
2
图 4-4 故障前中后仅含发电机内节点的导纳矩阵
三、求解电磁功率
得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率， 机械功率等于稳态的电磁功率中的有功分量。

所以可以有
P e=r eal (E *I)
如上中，E 为发电机内电势，I 为从发电机流出的电流。

但在参考文献 Ramnarayan Patel, T. S. Bhatti and D. P . Kothari.MA TLAB/Simulink-based transient stability analysis of a multimachine power system 中给出的电磁功率计算公式为：
P ei = E i
n
G ii + ∑
E i E j Y ij cos(θ ij - δ i + δ j ) j =1 j ≠i
稳态情况下有，机械功率 Pme=Pe
四、求解运动方程 发电机的运动方程可以写成常微分方程组：
o
其中 Pmi 为第 i 个机组故障前稳态的电磁功率。

在本次仿真中 D j ωi 为零，即阻尼为零。

t=0仿真开始，t=1 时引入故障，1.083s 后切除故障,t=3.083s 仿真结束。

求解运动方程后得到曲线如下：
五、结果分析
上图分别显示了各台发电机的转子角与时间的关系曲线，显示了发电机转速差
的曲线，和δ 21 = δ 2 - δ1 、δ31 = δ3 - δ1 的曲线，由图可以看到，最大角差δ21 为 85，出现在 t = 0.4s 处，无论是δ21 还是δ31 第二个摇摆都不大于第一个摇摆，可见系统是稳定的。

六、程序代码
主程序：
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.30
6 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 MV A xd'*10000 node xd xq xl xad xaq xf td0' rf
gen=[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)+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=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=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,YY0,YY2,YY3]=Ysim(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))+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);
for i=1:3
Pm(1,i)=N(i,8)./100;
end
options=odeset('RElTOL',1e-10); %设置精度
XX=[GEgj(1),GEgj(2),GEgj(3),2*pi*60*ones(1,3)];
tt=0;t0=1; tspan0=[tt,t0]; Yrun=Yc;
[T0,Y0]=ode45('fzz',tspan0,XX,options);
X0=Y0(end,:);t0=1; tc=1.083; tspan1=[t0,tc]; Yrun=Yb;
[T1,Y1]=ode45('fzz',tspan1,X0,options);
X1=Y1(end,:); tf=3.083; tspan2=[tc,tf]; Yrun=Ya;
[T2,Y2]=ode45('fzz',tspan2,X1,options);
delta1=Y0(:,1:3);
delta2=Y1(:,1:3);
delta3=Y2(:,1:3);
E1=repmat([1.05662336337847,1.05012618834491,1.01677591008947],45,1); E2=repmat([1.05662336337847,1.05012618834491,1.01677591008947],41,1); E3=repmat([1.05662336337847,1.05012618834491,1.01677591008947],369,1); e1=E1.*exp(delta1*1i);
e2=E2.*exp(delta2*1i);
e3=E3.*exp(delta3*1i);
V1=-(YY0(4:end,4:end)^-1)*(YY0(4:end,1:3))*e1';
V21=-(YY2(4:end,4:end)^-1)*YY2(4:end,1:3)*e2';
V2=[V21(1:6,:);zeros(1,length(V21));V21(7:8,:)];
V3=-(YY3(4:end,4:end)^-1)*YY3(4:end,1:3)*e3';
T=[T0;T1;T2]; Y12=[Y0;Y1;Y2];
V0=[V1,V2,V3]';
subplot(3,2,4);
plot(T,abs(V0));
subplot(3,2,5);
plot(T,unwrap(angle(V0))*180/pi);
subplot(3,2,1);
plot(T,Y12(:,1)*180/pi,T,Y12(:,2)*180/pi,T,Y12(:,3)*180/pi); %发电机功角
subplot(3,2,2);
plot(T,Y12(:,4),T,Y12(:,5),T,Y12(:,6)); %发电机转速
subplot(3,2,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)=-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
for 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-real(Y(i,j))*((N(i,4))^2);
end
end
end %N
for 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
end
end %K
for 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(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=[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;
end
for i=4:9
dQ(i-3,1)=(N(i,9)-N(i,7))/100;
end
for i=2:9
for j =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
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)));
end
end
故障前中后仅含发电机内节点的导纳矩阵：
function [Yp,Yd,Ya,YY0,YY2,YY3] = Ysim(gen,N,B,Y)
%故障前中后仅含发电机内节点的导纳矩阵：
% 此处显示详细说明
Yp=zeros(3);Yd=zeros(3);Ya=zeros(3); YY1=ones(12);YY2=zeros(11);YY3=zeros(12); YY1=[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; %负载等效导纳
YY1(i,i)= -(1/gen(i,3)*10000)*1i; YY1(i,i+3)=(1/gen(i,3)*10000)*1i; YY1(i+3,i)=YY1(i,i+3); end
for i=1:3
YY1(i+3,i+3)=YY1(i+3,i+3)-(1/gen(i,3)*10000)*1i; %发电机节点修改
end
YY1(3+5,3+5)=YY1(3+5,3+5)+Y5; YY1(3+6,3+6)=YY1(3+6,3+6)+Y6;
YY1(3+8,3+8)=YY1(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=YY1(1:3,1:3);Ynr=YY1(1:3,4:12);Yrn=YY1(4:12,1:3);Yrr=YY1(4:12,4:12); Yp=Ynn-Ynr*(Yrr^-1)*Yrn;
YY0=YY1;
YY3=YY1;
YY3(8,8)=YY3(8,8)-(1./(B(6,3)+B(6,4)*1i)+B(6,5)*1i/2);YY3(10,10)=YY3(10,10)-
(1./(B(6,3)+B(6,4)*1i)+B(6,5)*1i/2);YY3(8,10)=0;YY3(10,8)=0; %故障后增广
YY1(10,:)=[]; YY1(:,10)=[];
YY2=YY1; %故障中增广
Ynn=YY2(1:3,1:3);Ynr1=YY2(1:3,4:11);Yrn1=YY2(4:11,1:3);Yrr1=YY2(4:11,4:11); Yd=Ynn-Ynr1*(Yrr1^-1)*Yrn1;
Ynn=YY3(1:3,1:3);Ynr=YY3(1:3,4:12);Yrn=YY3(4:12,1:3);Yrr=YY3(4:12,4:12); Ya=Ynn-Ynr*(Yrr^-1)*Yrn;
end
求解运动方程：
function xtt=fzz(T,DX)
global Pm Yrun gen GenE
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];
七、总结
首先感谢任课老师江宁强老师对我仿真的指导与帮助，还要感谢在程序编写过程中给过我无私帮助的其他同学。

这次课程的仿真设计实验让我学到了很多，本文中切断时间在1.083s，最终的图形和老师的基本一致。

刚开始y 举证先是求的不准确,少了对地导纳，之后求解雅克比矩阵分块求解不是太准确，但最终结果终于达到了要求。

对于电磁功率如果按照文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothari.MA TLAB/Simulink-based transient stability analysis of a multimachine power system 中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，简化了电磁功率的计算，这样一来计算量就少了很多。