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 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］；
三、仿真框图
在仿真之前，首先，应明确仿真的所要到达的结果,即仿真目标：本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。故障前，系统处于稳定状态，发电机的攻角、转速基本稳定.而当系统发生故障以及故障切除，系统结构拓扑发生变化,系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的
状态量为发电机攻角、发电机转速。
其次，当明确仿真目标后，我们就得明确大体的仿真框架流程.
仿真框架流程如下：
数据准备
(支路、节点、发电机)

后处理
调用ode45计算发电攻角、
转速变化情况

列写系统状态方程
（转子运动方程）

计算故障前中后发电机
内节点的导纳矩阵

发电机初值计算
潮流计算

图3—1仿真流程图
四、 仿真模型
在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。
一、 潮流计算
由于本文以三机九节点为模型，假定节点一为参考节点，这样就有2两个发电机的PV节点，6个PQ节点，未知量为8个节点（包括2个PV节点和6个PQ节点)的电压相角，还有6个节点（PQ节点）的电压幅值。
可以先求出Y矩阵

图4—1 Y矩阵
然后,我们列写方程,也就是利用各个节点的有功、无功功率的平衡关系，列写14个功率平衡方程。这样就能使用牛顿一拉夫逊算法来求解这14个非线性方程.
其中的关键是要计算出雅克比矩阵

图4-2雅克比矩阵
然后计算出修正量.在设定精度和最大迭代次数的前提下进行迭代,直到满足要求。 电力网络的节点功率方程可用一般形式表示如下：
牛顿拉夫逊算法修正方程
W = -JΔV
其中W是节点不平衡量向量，包括有功,无功,电压；J是雅克比矩阵；ΔV是节点电压修正量.
令
，
则极坐标形式的功率不平衡量方程

雅可比矩阵J各元素的表达式

当j≠i时: 当j=i时：

其中，
。
进行牛顿拉夫逊算法得到潮流结果

图4—3潮流结果
二、故障前中后仅含发电机内节点的导纳矩阵

图4—4故障前中后仅含发电机内节点的导纳矩阵
三、 求解电磁功率
得到故障前，故障中，故障后三个不同的导纳矩阵后,就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有
Pe=real(E＊I）
如上中，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个机组故障前稳态的电磁功率.在本次仿真中Djωi为零，即阻尼为零.仿真开始，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' 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]=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））+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)； ％设置精度

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(:，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）=—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]=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）;
end

for i=1:3
Y1（i+3,i+3）=Y1(i+3，i+3）-(1/gen（i,3)*10000）＊1i; %发电机节点修改
end

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