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

4 0 0 1.0000 0.00 0.00 0.00 0.000.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
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));
三机九节点电力系统
暂态仿真
学院：
专业：
学号：
姓名：
授课教师：
一、摘要
电力系统仿真计算己经成为电力系统设计、运行与控制中不可缺少的手段。通过设置不同故障类型、不同故障地点运用仿真技术可以对电力系统的暂态稳定进行分析。本文采用IEEE 3机9节点的经典多机模型，基于隐式梯形积分法对系统发生三相金属性短路故障进行仿真，分析系统在这种情况下的暂态稳定。发电机模型采用经典的二阶模型;负荷采用恒定阻抗负荷。在Matlab2010上编写程序进行调试和运行。
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;
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 01.0000 0.00 100.00 35.00 0.00 0.00 0.00 1.016
电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。在动态稳定仿真中使用简单的电力系统模型，通过仿真计算分析说明，此仿真方法可以进行简单的电力系统暂态分析，对提高电力系统暂态稳定具有重要意义。
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);
图2-1 3机9节点系统
系统数据
节点号有无负载类型电压相角有功负荷无功负荷有功出力无功出力电压基准期望电压
N=[1 0 3 1.0400 0.00 0.00 0.0071.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-2雅克比矩阵
然后计算出修正量。在设定精度和最大迭代次数的前提下进行迭代，直到满足要求。电力网络的节点功率方程可用一般形式表示如下:
牛顿拉夫逊算法修正方程
W = -JΔV
其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵；ΔV是节点电压修正量。
令
，
则极坐标形式的功率不平衡量方程
一、潮流计算
由于本文以三机九节点为模型，假定节点一为参考节点，这样就有2两个发电机的PV节点，6个PQ节点，未知量为8个节点(包括2个PV节点和6个PQ节点)的电压相角，还有6个节点(PQ节点)的电压幅值。
可以先求出Y矩阵
图4-1 Y矩阵
然后，我们列写方程，也就是利用各个节点的有功、无功功率的平衡关系，列写14个功率平衡方程。这样就能使用牛顿一拉夫逊算法来求解这14个非线性方程。
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.01 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
雅可比矩阵J各元素的表达式
当j≠i时：当j=i时：
其中，
。
进行牛顿拉夫逊算法得到潮流结果
图4-3潮流结果
二、故障前中后仅含发电机内节点的导纳矩阵
图4-4故障前中后仅含发电机内节点的导纳矩阵
三、求解电磁功率
得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有
2 2 0.8958 0.1198 6.00 0.8645 0.1969 0.53512.80 0.0625
33 1.3125 0.1813 8.59 1.2578 0.2500 0.6006.02 0.0585];
三、仿真框图
在仿真之前，首先，应明确仿真的所要到达的结果，即仿真目标:本此仿真的结果主要是得到发电机攻角、转速随时间变化的值，包括故障前、故障中、故障后。故障前，系统处于稳定状态，发电机的攻角、转速基本稳定。而当系统发生故障以及故障切除，系统结构拓扑发生变化，系统的状态也将随时间发生变化，为了求取系统状态的变化，我们通过对系统进行简化建立数学模型，得到相关的代数一微分方程组，进行数值计算，从而得到系统状态的随时间的变化值。此次仿真的系统以发电机二阶经典模型来进行系统是数学建模，系统的
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
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))))));
状态量为发电机攻角、发电机转速。
其次，当明确仿真目标后，我们就得明确大体的仿真框架流程。
仿真框架流程如下：
图3-1仿真流程图
四、仿真模型
在电力系统的机电暂态仿真中，常根据实际要求的不同，采用不同时间尺度的仿真模型，而仿真算法和采用的模型有直接的关系，下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。
%支路数据
%从到电阻电抗容纳类型变比
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.0170.092 0.158 0 0
5 7 0.032 0.161 0.306 0 0
4 6 0.017 0.092 0.158 0 0
57 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'*10000node xd xq xl xad xaq xf td0' rf
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
Pe=real(E*I)
如上中，E为发电机内电势，I为从发电机流出的电流。
但在参考文献Ramnarayan Patel, T. S. Bhatti and D.P. Kothari.MATLAB/Simulink-based transient stabilityanalysis of a multimachine power system中给出的电磁功率计算公式为：
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 00];
发电机数据如下：
%发电机母线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
二、案例
本次课程主要应用P. M. Anderson and A. A. Fouad编写的《Power System Control and Stability》一书中所引用的Western System Coordinated Council (WSCC)三机九节点系统模型。
系统电路结构拓扑图如下：
7 0 0 1.0000 0.00 0.00 0.00 0.000.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];
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 01.0000 0.00 125.00 50.00 0.00 0.00 0.00 0.996
稳态情况下有，机械功率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 21.0250 0.00 0.0来自百度文库 0.00 163.00 6.70 18.00 1.025
[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;
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