三机九节点电力系统仿真matlab.docx
《三机九节点电力系统仿真matlab.docx》由会员分享，可在线阅读，更多相关《三机九节点电力系统仿真matlab.docx（26页珍藏版）》请在冰豆网上搜索。
三机九节点电力系统仿真matlab
电力系统仿真作业------------
三机九节点电力系统
暂态仿真
）
学院：
能源与动力工程学院
专业：
电力系统及其自动化
学号：
姓名：
于永生
导师：
授课教师：
一、概述
在动态稳定分析中，系统由线性化的微分方程组和代数方程组描写，并用经典的或现代的线性系统理论来进行稳定分析，分析可以在时域或频域进行。
当用计算机和现代线性系统理论分析时，常把系统线性化的微分方程组和代数方程组消去代数变量，化为状态方程形式，并广泛采用特征分析进行稳定分析。
电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。
其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。
由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。
在动态稳定仿真中使用简单的电力系统模型，发电机用三阶模型表示。
二、课程主要任务
本次课程主要应用P.M.AndersonandA.A.Fouad编写的《PowerSystemControlandStability》一书中所引用的WesternSystemCoordinatedCouncil（WSCC）三机九节点系统模型。
1.系统数据
其中，节点数据如下：
%节点数据
%节点电压电压发电机发电机负荷负荷节点
%号幅值相角有功无功有功无功类型（1PQ2PV3平衡）
N=[10003
20002
30002
41000001
510001
610001
71000001
8100011
91000001];
其中，支路数据如下：
%线路数据
%首端末端电阻电抗电纳（1/2）变压器非标准变比
L=[451
461
571
691
781
891
14001
27001
39001];
发电机数据如下：
%发电机母线XdXd'Td0'XqXq'Tq0’TjXf
Ge=[110
22
33];
系统电路结构拓扑图如下：
图1WSCC3机9节点系统
（所有参数以100MVA为基准值的标幺值）
2.潮流计算
首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。
为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
牛顿拉夫逊算法修正方程
W=-JΔV
其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵；ΔV是节点电压修正量。
令
，
则直角坐标形式的功率不平衡量方程为
PQ节点：
PV节点：
极坐标形式的功率不平衡量方程
雅可比矩阵J各元素的表达式
当j≠i时：
当j=i时：
其中，
。
进行牛顿拉夫逊算法迭代后得到电压幅值V和相角θ。
3.负荷等效和支路简化
然后求出支路电流，将发电机内电抗X’加入系统导纳矩阵，求出发电机内电势E’。
加入发电机内节点后，系统导纳矩阵变成12*12阶的矩阵，并将负荷等效成阻抗。
然后将支路导纳矩阵分块，如下：
其中，A是3*3的方阵，E是3*9的矩阵，C是9*3的矩阵，D是9*9的方阵。
经过网络简化得到故障前的3*3简化导纳矩阵
（1）
其中“inv（D）”是MATLAB中D矩阵的求逆。
故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有
此时，A是3*3的方阵，E是3*8的矩阵，C是8*3的矩阵，D是8*8的方阵。
简化矩阵求法如同公式
（1）。
故障后的Y矩阵相对于故障前的Y矩阵只是第5个节点和第7个节点有变化，反映到12*12的矩阵中即为第（10,8），（8,10），（8,8），（10,10）位置的元素有变化，其中（10,8），（8,10）位置的元素变为零，（8,8），（10,10）节点在故障前的基础上加上（8,10）位置元素的值。
然后简化导纳矩阵的求法同式
（1）。
4.求解电磁功率
得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。
所以可以有
Pe=real（E*I）
（2）
公式
（2）中，E为发电机内电势，I为从发电机流出的电流。
但在参考文献RamnarayanPatel,T.S.BhattiandD.P.Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式为：
本人在此次仿真中用的是公式
（2）计算得到的电磁功率。
稳态情况下有，机械功率
Pme=Pe
5.求解运动方程
发电机的运动方程可以写成常微分方程组：
其中Pmi为第i个机组故障前稳态的电磁功率。
在本次仿真中Djωi为零，即阻尼为零。
仿真开始，t=0时引入故障，后切除故障。
求解运动方程后得到ω和δ曲线如下：
图1ω曲线
图2δ曲线
6.程序清单
以下是我编写的仿真程序，除主程序外还包含三个函数：
①极坐标转换成直角坐标函数pol2rect（V,del）,其中输入参数V为极坐标向量的幅值，del为极坐标向量的角度，函数返回值为一个复数，即为极坐标向量在直角坐标中的复数值；②直角坐标转换成极坐标函数rect2pol（Z）,其中输入参数Z为一个复数，函数返回值为一个二维向量，向量的第一个数为幅值，第二个数为相角；③求解微分方程所用的得到微分量的函数Gen_fw（t,X,Y_Gen,E,Pm0,Tj）,其中输入参数X为ω和δ的迭代初值，Y_Gen，为求解所用的导纳矩阵，这里是三阶的方阵，对应于故障前，故障中，和故障后的三个Y矩阵，E为发电机内电势，Pm0为机械功率，等于稳态时的有功功率，Tj为运动方程中的2H。
（1）.主程序：
clc;
clear;
%节点数据
%节点电压电压发电机发电机负荷负荷节点
%号幅值相角有功无功有功无功类型（1PQ2PV3平衡
%N=[10003
%20002
%30002
%41000001
%510001
%610001
%71000001
%8100011
%91000001];
%线路数据
%首端末端电阻电抗电纳（1/2）变压器非标准变比
L=[451
461
571
691
781
891
14001
27001
39001];
%形成节点导纳矩阵
fb=L（:
1）;%起始节点编号
tb=L（:
2）;%末节点编号
r=L（:
3）;%电阻
x=L（:
4）;%电抗
b=L（:
5）;%电纳
a=L（:
6）;%变压器非标准变比
z=r+i*x;%z矩阵，复数，9乘1的矩阵
y=1./z;%求倒数，9乘1的导纳矩阵
b=i*b;%虚数
nb=max（max（fb）,max（tb））;%nb=9
nl=length（fb）;%支路个数
Y=zeros（nb,nb）;%Y是9乘9的支路矩阵
fork=1:
nl
Y（fb（k）,tb（k））=Y（fb（k）,tb（k））-y（k）/a（k）;%Y方阵的元素，非对角元（4,5）（4，6）（5，7）（6,9）（7，8）（8,9）（1,4）（2,7）（3,9）
Y（tb（k）,fb（k））=Y（fb（k）,tb（k））;%Y方阵的元素，上一句的对称部分（5,4）（6,4）（7，5）（9，6）（8，7）（9,8）（4,1）（7,2）（9,3）
end
form=1:
nb
forn=1:
nl
iffb（n）==m%支路中的第n个数的首节点编号等于m
Y（m,m）=Y（m,m）+y（n）/（a（n）^2）+b（n）;%Y矩阵中的对角元
elseiftb（n）==m%支路中的第n个数的末节点编号等于m
Y（m,m）=Y（m,m）+y（n）+b（n）;%因为支路是首端或末端等于m，不能是首端编号等于末端编号
end
end
end
%节点数据
%节点电压电压发电机发电机负荷负荷节点
%号幅值相角有功无功有功无功类型（1PQ2PV3平衡
N=[10003
20002
30002
41000001
510001
610001
71000001
8100011
91000001];
nbus=9;%节点数
bus=N（:
1）;%节点编号
type=N（:
8）;%节点类型
V=N（:
2）;%额定电压
del=N（:
3）;%电压相角
Pg=N（:
4）;%发电机有功功率
Qg=N（:
5）;%发电机无功功率
Pl=N（:
6）;%负荷有功功率
Ql=N（:
7）;%负荷无功功率
P=Pg-Pl;%发出减消耗
Q=Qg-Ql;%发出减消耗
Psp=P;%额定有功功率
Qsp=Q;%额定无功功率
G=real（Y）;
B=imag（Y）;
pv=find（type==3|type==2）;%pv节点编号和平衡节点编号
pq=find（type==1）;%pq节点编号
npv=length（pv）;%npv=3
npq=length（pq）;%npq=6
Tol=1;
Iter=1;
while（Tol>1e-7）
P=zeros（nbus,1）;%迭代初值P为零
Q=zeros（nbus,1）;%迭代初值Q为零
fori=1:
nbus
fork=1:
nbus
P（i）=P（i）+V（i）*V（k）*（G（i,k）*cos（del（i）-del（k））+B（i,k）*sin（del（i）-del（k）））;
Q（i）=Q（i）+V（i）*V（k）*（G（i,k）*sin（del（i）-del（k））-B（i,k）*cos（del（i）-del（k）））;
end
end
dPa=Psp-P;%额定有功功率减去流出或流入的有功功率
dQa=Qsp-Q;%额定无功功率减去流出或流入的无功功率
k=1;
dQ=zeros（npq,1）;%pq节点的无功功率，pq乘1的矩阵,6乘1
fori=1:
nbus
iftype（i）==1%pq节点
dQ（k,1）=dQa（i）;
k=k+1;
end
end
dP=dPa（2:
nbus）;
Z=[dP;dQ];
%下面计算雅克比矩阵
J1=zeros（nbus-1,nbus-1）;%有功功率对功角求偏微分
fori=1:
（nbus-1）
m=i+1;
fork=1:
（nbus-1）
n=k+1;
ifn==m
forn=1:
nbus
J1（i,k）=J1（i,k）+V（m）*V（n）*（-G（m,n）*sin（del（m）-del（n））+B（m,n）*cos（del（m）-del（n）））;
end
J1（i,k）=J1（i,k）-V（m）^2*B（m,m）;
else
J1（i,k）=V（m）*V（n）*（G（m,n）*sin（del（m）-del（n））-B（m,n）*cos（del（m）-del（n）））;
end
end
end
J2=zeros（nbus-1,npq）;%有功功率对电压求偏微分
fori=1:
（nbus-1）
m=i+1;
fork=1:
npq
n=pq（k）;
ifn==m
forn=1:
nbus
J2（i,k）=J2（i,k）+V（n）*（G（m,n）*cos（del（m）-del（n））+B（m,n）*sin（del（m）-del（n）））;
end
J2（i,k）=J2（i,k）+V（m）*G（m,m）;
else
J2（i,k）=V（m）*（G（m,n）*cos（del（m）-del（n））+B（m,n）*sin（del（m）-del（n）））;
end
end
end
J3=zeros（npq,nbus-1）;%无功功率对相角求偏微分
fori=1:
npq
m=pq（i）;
fork=1:
（nbus-1）
n=k+1;
ifn==m
forn=1:
nbus
J3（i,k）=J3（i,k）+V（m）*V（n）*（G（m,n）*cos（del（m）-del（n））+B（m,n）*sin（del（m）-del（n）））;
end
J3（i,k）=J3（i,k）-V（m）^2*G（m,m）;
else
J3（i,k）=V（m）*V（n）*（-G（m,n）*cos（del（m）-del（n））-B（m,n）*sin（del（m）-del（n）））;
end
end
end
J4=zeros（npq,npq）;%无功功率对电压求偏微分
fori=1:
npq
m=pq（i）;
fork=1:
npq
n=pq（k）;
ifn==m
forn=1:
nbus
J4（i,k）=J4（i,k）+V（n）*（G（m,n）*sin（del（m）-del（n））-B（m,n）*cos（del（m）-del（n）））;
end
J4（i,k）=J4（i,k）-V（m）*B（m,m）;
else
J4（i,k）=V（m）*（G（m,n）*sin（del（m）-del（n））-B（m,n）*cos（del（m）-del（n）））;
end
end
end
J=[J1J2;J3J4];
%X=inv（J）*Z;
X=J\Z;
dTh=X（1:
nbus-1）;
dV=X（nbus:
end）;
del（2:
nbus）=dTh+del（2:
nbus）;
k=1;
fori=2:
nbus
iftype（i）==1
V（i）=V（i）+dV（i-3）;
k=k+1;
end
end
Iter=Iter+1;
Tol=max（abs（Z））;
end
Vm=pol2rect（V,del）;%极坐标转换成直角坐标
Del=180/pi*del;%弧度转化成角度
Iij=zeros（nb,nb）;%支路电流
form=1:
nl%每一条支路都过一遍
p=fb（m）;q=tb（m）;
Iij（p,q）=-（Vm（p）-Vm（q））*Y（p,q）;%Y（m,n）=-y（m,n）..
Iij（q,p）=-Iij（p,q）;
end
Iij=sparse（Iij）;
%线路数据
%首端末端电阻电抗电纳（1/2）变压器非标准变比
L=[451
461
571
691
781
891
14001
27001
39001];
Yzengjia=zeros（12,12）;
Zci=[j*j*j*];
Yci=1./Zci;
fork=1:
3
Yzengjia（k,k）=Yci（k）;
end
fork=1:
9
form=1:
9
Yzengjia（3+k,3+m）=Y（k,m）;
end
end
Yzengjia（1,4）=-Yzengjia（1,1）;
Yzengjia（4,1）=Yzengjia（1,4）;
Yzengjia（2,5）=-Yzengjia（2,2）;
Yzengjia（5,2）=Yzengjia（2,5）;
Yzengjia（3,6）=-Yzengjia（3,3）;
Yzengjia（6,3）=Yzengjia（3,6）;
Yzengjia（4,4）=Yzengjia（4,4）+Yzengjia（1,1）;
Yzengjia（5,5）=Yzengjia（5,5）+Yzengjia（2,2）;
Yzengjia（6,6）=Yzengjia（6,6）+Yzengjia（3,3）;
%求负荷等效导纳
Y5=N（5,6）/（V（5））^2-j*（N（5,7）/（V（5））^2）;
Y6=N（6,6）/（V（6））^2-j*（N（6,7）/（V（5））^2）;
Y8=N（8,6）/（V（8））^2-j*（N（8,7）/（V（5））^2）;
Yzengjia（8,8）=Yzengjia（8,8）+Y5;
Yzengjia（9,9）=Yzengjia（9,9）+Y6;
Yzengjia（11,11）=Yzengjia（11,11）+Y8;
A=zeros（3,3）;
form=1:
3
forn=1:
3
A（m,n）=Yzengjia（m,n）;
end
end
E=zeros（3,9）;
form=1:
3
forn=4:
12
E（m,n-3）=Yzengjia（m,n）;
end
end
C=zeros（9,3）;
form=4:
12
forn=1:
3
C（m-3,n）=Yzengjia（m,n）;
end
end
D=zeros（9,9）;
form=4:
12
forn=4:
12
D（m-3,n-3）=Yzengjia（m,n）;
end
end
Ypre=A-E*（inv（D））*C;%得到故障前的3*3矩阵
Yzengjiadur=Yzengjia;
Yzengjiadur（:
[10]）=[];
Yzengjiadur（[10],:
）=[];
Adur=A;
Edur=zeros（3,8）;
form=1:
3
forn=1:
8
Edur（m,n）=Yzengjiadur（m,n+3）;%因为Edur的第四到第九列是零，而所以这里为了方便，直接取E的前六列
end
end
Cdur=zeros（8,3）;
form=1:
8
forn=1:
3
Cdur（m,n）=Yzengjiadur（m+3,n）;%因为Cdur的第四到第九行是零，而所以这里为了方便，直接取C的前六行
end
end
Ddur=zeros（8,8）;
form=1:
8
forn=1:
8
Ddur（m,n）=Yzengjiadur（m+3,n+3）;
end
end
Ydur=Adur-Edur*（inv（Ddur））*Cdur;%故障中的3*3矩阵
Yzengjiaaft=Yzengjia;
temp=Yzengjiaaft（8,10）;
Yzengjiaaft（8,10）=0;
Yzengjiaaft（8,8）=Yzengjiaaft（8,8）+temp;
Yzengjiaaft（10,8）=0;
Yzengjiaaft（10,10）=Yzengjiaaft（10,10）+temp;
Aaft=zeros（3,3）;
form=1:
3
forn=1:
3
Aaft（m,n）=Yzengjiaaft（m,n）;
end
end
Eaft=zeros（3,9）;
form=1:
3
forn=4:
12
Eaft（m,n-3）=Yzengjiaaft（m,n）;
end
end
Caft=zeros（9,3）;
form=4:
12
forn=1:
3
Caft（m-3,n）=Yzengjiaaft（m,n）;
end
end
Daft=zeros（9,9）;
form=1:
9
forn=1:
9
Daft（m,n）=Yzengjiaaft（m+3,n+3）;
end
end
Yaft=Aaft-Eaft*（inv（Daft））*Caft;%故障后的3*3矩阵
%求电磁功率
Enei=ones（3,1）;%发电机内电势
Pe=ones（3,1）;
Enei
（2）=Vm
（2）+Zci
（2）*Iij（2,7）;
Enei
（1）=Vm
（1）+Zci
（1）*（-Iij（4,1））;
Enei（3）=Vm（3）+Zci（3）*Iij（3,9）;
polEnei1=rect2pol（Enei
（1））;%直角坐标转换成极坐标后得到角度
polEnei2=rect2pol（Enei
（2））;
polEnei3=rect2pol（Enei（3））;
Pe
（1）=real（Enei
（1）*（（-Iij（4,1））'））;%故障前的电气功率
Pe
（2）=real（Enei
（2）*（（Iij（2,7））'））;
Pe（3）=real（Enei（3）*（Iij（3,9）））;
delta=zeros（3,1）;
delta
（1）=polEnei1
（2）;
delta
（2）=polEnei2
（2）;
delta（3）=polEnei3
（2）;
Tj=[;;];
Pm0=Pe;%机械功率=故障前的电磁功率
X0=[ones（3,1）.*（2*pi*60）,delta];
Y_Gen=Ydur;E=abs（Enei）;
Time_s=[0,];%后切除故障
[t1,Xout1]=ode45（@（t,X）Gen_fw（t,X,Y_Gen,E,Pm0,Tj）,Time_s,X0）;
Time_s=[,2];
Y_Gen=Yaft;X0=Xout1（end,:
）;
[t2,Xout2]=ode45（@（t,X）Gen_fw（t,X,Y_Gen,E,Pm0,Tj）,Time_s,X0）;
升级会员 