三机九节点电力系统仿真matlab.docx
《三机九节点电力系统仿真matlab.docx》由会员分享,可在线阅读,更多相关《三机九节点电力系统仿真matlab.docx(27页珍藏版)》请在冰点文库上搜索。
三机九节点电力系统仿真matlab
电力系统仿真作业-—---—---———
三机九节点电力系统
暂态仿真
学院:
能源与动力工程学院
专业:
电力系统及其自动化
学号:
姓名:
于永生
导师:
授课教师:
一、概述
在动态稳定分析中,系统由线性化的微分方程组和代数方程组描写,并用经典的或现代的线性系统理论来进行稳定分析,分析可以在时域或频域进行。
当用计算机和现代线性系统理论分析时,常把系统线性化的微分方程组和代数方程组消去代数变量,化为状态方程形式,并广泛采用特征分析进行稳定分析.
电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。
其暂态过渡过程不仅包括电磁方面的过渡过程,而且还有机电方面的过渡过程。
由此可见,电力系统的数学模型是一个强非线性的高维状态方程组。
在动态稳定仿真中使用简单的电力系统模型,发电机用三阶模型表示。
二、课程主要任务
本次课程主要应用P。
M。
AndersonandA。
A.Fouad编写的《PowerSystemControlandStability》一书中所引用的WesternSystemCoordinatedCouncil(WSCC)三机九节点系统模型。
1.系统数据
其中,节点数据如下:
%节点数据
%节点电压电压发电机发电机负荷负荷节点
%号幅值相角有功无功有功无功类型(1PQ2PV3平衡)
N=[11。
0400.71640。
2705003
21。
02501.630。
0665002
31。
02500。
85—0.1086002
41000001
510001.250。
51
610000。
90。
31
71000001
8100010。
351
91000001];
其中,支路数据如下:
%线路数据
%首端末端电阻电抗电纳(1/2)变压器非标准变比
L=[450。
010。
0850.0881
460。
0170.0920。
0791
570。
0320。
1610.1531
690。
0390。
170。
1791
780.00850。
0720。
07451
890.01190.10080.10541
1400.057601
2700。
062501
3900。
058601];
发电机数据如下:
%发电机母线XdXd'Td0'XqXq’Tq0’TjXf
Ge=[110.14600.06088。
960.09690.0969047.280.0576
220.89580。
11986。
000.86450。
19690.53512.800.0625
331。
31250。
18138.591.25780.25000。
6006。
020.0585];
系统电路结构拓扑图如下:
图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。
Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式为:
本人在此次仿真中用的是公式
(2)计算得到的电磁功率。
稳态情况下有,机械功率
Pme=Pe
5.求解运动方程
发电机的运动方程可以写成常微分方程组:
其中Pmi为第i个机组故障前稳态的电磁功率。
在本次仿真中Djωi为零,即阻尼为零.仿真开始,t=0时引入故障,0.083s后切除故障。
求解运动方程后得到ω和δ曲线如下:
图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=[11。
0400。
71640.2705003
%21.02501.630.0665002
%31.02500.85-0.1086002
%41000001
%510001。
250。
51
%610000.90。
31
%71000001
%8100010.351
%91000001];
%线路数据
%首端末端电阻电抗电纳(1/2)变压器非标准变比
L=[450。
010。
0850.0881
460。
0170.0920。
0791
570。
0320.1610.1531
690。
0390。
170.1791
780。
00850。
0720.07451
890。
01190。
10080。
10541
1400。
057601
2700.062501
3900。
058601];
%形成节点导纳矩阵
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=[11.0400。
71640。
2705003
21.02501。
630.0665002
31。
02500.85—0。
1086002
41000001
510001.250。
51
610000。
90.31
71000001
8100010.351
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=[450。
010.0850.0881
460.0170。
0920。
0791
570.0320.1610.1531
690.0390.170。
1791
780.00850。
0720。
07451
890.01190.10080。
10541
1400.057601
2700。
062501
3900。
058601];
Yzengjia=zeros(12,12);
Zci=[j*0。
0608j*0。
1198j*0.1813];
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-Ea