matlab三机九节点电力系统仿真带程序Word下载.docx
《matlab三机九节点电力系统仿真带程序Word下载.docx》由会员分享,可在线阅读,更多相关《matlab三机九节点电力系统仿真带程序Word下载.docx(17页珍藏版)》请在冰豆网上搜索。
7001.00000.000.000.000.000.00230.001.026
8101.00000.00100.0035.000.000.000.001.016
9001.00000.000.000.000.000.00230.001.032];
%支路数据
%从到电阻电抗容纳类型变比
B=[140.00.05760.011
270.00.06250.011
390.00.05860.011
450.0100.0850.17600
460.0170.0920.15800
570.0320.1610.30600
690.0390.1700.35800
780.00850.0720.14900
890.01190.10080.20900];
发电机数据如下:
%发电机母线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];
三、仿真框图
在仿真之前,首先,应明确仿真的所要到达的结果,即仿真目标:
本此仿真的结果主要是得到发电机攻角、转速随时间变化的值,包括故障前、故障中、故障后。
故障前,系统处于稳定状态,发电机的攻角、转速基本稳定。
而当系统发生故障以及故障切除,系统结构拓扑发生变化,系统的状态也将随时间发生变化,为了求取系统状态的变化,我们通过对系统进行简化建立数学模型,得到相关的代数一微分方程组,进行数值计算,从而得到系统状态的随时间的变化值。
此次仿真的系统以发电机二阶经典模型来进行系统是数学建模,系统的
状态量为发电机攻角、发电机转速。
其次,当明确仿真目标后,我们就得明确大体的仿真框架流程。
仿真框架流程如下:
图3-1仿真流程图
四、仿真模型
在电力系统的机电暂态仿真中,常根据实际要求的不同,采用不同时间尺度的仿真模型,而仿真算法和采用的模型有直接的关系,下面就本次仿真实例机电暂态过程的仿真模型及其仿真算法。
一、潮流计算
由于本文以三机九节点为模型,假定节点一为参考节点,这样就有2两个发电机的PV节点,6个PQ节点,未知量为8个节点(包括2个PV节点和6个PQ节点)的电压相角,还有6个节点(PQ节点)的电压幅值。
可以先求出Y矩阵
图4-1Y矩阵
然后,我们列写方程,也就是利用各个节点的有功、无功功率的平衡关系,列写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为从发电机流出的电流。
但在参考文献RamnarayanPatel,T.S.BhattiandD.P.Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式为:
稳态情况下有,机械功率Pme=Pe
四、求解运动方程
发电机的运动方程可以写成常微分方程组:
其中Pmi为第i个机组故障前稳态的电磁功率。
在本次仿真中Djωi为零,即阻尼为零。
仿真开始,t=0时引入故障,0.083s后切除故障。
求解运动方程后得到曲线如下:
五、结果分析
上图分别显示了各台发电机的转子角与时间的关系曲线,显示了发电机转速差的曲线,和
、
的曲线,由图可以看到,最大角差
为
,出现在
处,无论是
还是
第二个摇摆都不大于第一个摇摆,可见系统是稳定的。
六、程序代码
主程序:
globalPmYrungenGenE
%节点数据
%节点号有无负载类型电压相角有功负荷无功负荷有功负荷无功负荷电压基准期望电压
2021.02500.000.000.00163.006.7018.001.025
3021.02500.000.000.0085.00-10.9013.801.025
%发电机数据
%HMVAxd'
*10000nodexdxqxlxadxaqxftd0'
rf
gen=[2364247.560810.14600.09690.03360.11240.06330.14838.960.0000439
640192.0119820.89580.86450.05210.84370.81240.91736.000.0004054
301128.0181331.31251.25780.07421.23831.28361.35555.890.0006105];
Y=zeros(9,9);
%导纳矩阵
fori=1:
9
a=B(i,1);
b=B(i,2);
ifB(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%导纳矩阵
forT=1:
100
[dP,dQ]=Caoliu(N,Y);
%潮流
J=Ykb(N,Y);
%雅克比矩阵
U=zeros(6);
fori=4:
U(i-3,i-3)=N(i,4);
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
[Yc,Yb,Ya]=Ynew(gen,N,B,Y);
GEgj=zeros(1,3);
GenE=zeros(1,3);
fori=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))))));
Yrun=zeros(3);
Pm=zeros(1,3);
Pm(1,i)=N(i,8)./100;
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('
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);
5)-Y12(:
4),T,Y12(:
6)-Y12(:
4));
%发电机转速差
subplot(3,1,3);
2)*180/pi-Y12(:
3)*180/pi-Y12(:
1)*180/pi)%发电机攻角差
雅克比矩阵:
functionJ=Ykb(N,Y)
H=zeros(8);
N1=zeros(8,6);
K=zeros(6,8);
L=zeros(6);
fori=2:
forj=2:
ifi~=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)));
H(i-1,j-1)=(N(i,9)-N(i,7))/100+imag(Y(i,j))*((N(i,4))^2);
end%H
forj=4:
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)));
N1(i-1,j-3)=-(N(i,8)-N(i,6))/100-real(Y(i,j))*((N(i,4))^2);
end%N
fori=4:
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)));
K(i-3,j-1)=-(N(i,8)-N(i,6))/100+real(Y(i,j))*((N(i,4))^2);
end%K
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)));
L(i-3,j-3)=-(N(i,9)-N(i,7))/100+imag(Y(i,j))*((N(i,4))^2);
end
end%L
J=[HN1;
KL];
潮流计算
function[dP,dQ]=Caoliu(N,Y)
dP=zeros(8,1);
dQ=zeros(6,1);
dP(i-1,1)=(N(i,8)-N(i,6))/100;
dQ(i-3,1)=(N(i,9)-N(i,7))/100;
forj=1:
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
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)));
故障前中后仅含发电机内节点的导纳矩阵:
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;
%负载等效导纳
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);
Y1(i+3,i+3)=Y1(i+3,i+3)-(1/gen(i,3)*10000)*1i;
%发电机节点修改
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:
Yrr=Y1(4:
12,4:
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:
Ynr1=Y2(1:
11);
Yrn1=Y2(4:
11,1:
Yrr1=Y2(4:
11,4:
Yd=Ynn-Ynr1*(Yrr1^-1)*Yrn1;
Ynn=Y3(1:
Ynr=Y3(1:
Yrn=Y3(4:
Yrr=Y3(4:
Ya=Ynn-Ynr*(Yrr^-1)*Yrn;
求解运动方程:
functionxtt=fouad3(T,DX)
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举证先是求的不准确,少了对地导纳。
之后求解雅克比矩阵分块求解不是太准确。
对于电磁功率如果按照文献RamnarayanPatel,T.S.BhattiandD.P.Kothari.MATLAB/Simulink-basedtransientstabilityanalysisofamultimachinepowersystem中给出的电磁功率计算公式来计算会变得很繁琐,但经过其他同学的帮助与指点,简化了电磁功率的计算,这样一来计算量就少了很多。