暂态稳定性分析大程
作者 : 李志
学号 : 21410076
Tel : 13732254786
目录
[TOC]
1. 引言
对于地区性的电力系统来说，一般失去暂态稳定过程发展很快，通常分析系统遭受干扰后第一摇摆周期（1~1.5s）的机电暂态过程就可以判断系统是否能够维持稳定运行。这种情况下的暂态稳定分析中，由于调速系统的惯性，使得在短时间内原动机的功率不会发生很大变化，因此可以忽略调速系统的作用而假定原动机的功率保持不变；此外由于发电机励磁绕组的时间常数较大，这样在短时间内其磁链也不会发生显著变化，而对励磁调节系统的作用，可以用发电机暂态电势$E_q{’}$或$E{’}$保持恒定来近似模拟，即认为在第一摇摆周期内，励磁绕组中自由电流分量的变化由励磁调节系统的作用所补偿，从而使励磁绕组的磁链$Ψ_f$在这段时间内保持不变。相应地，阻尼绕组的影响将略去不计。
简单模型下的暂态稳定分析程序在电力系统运行与规划中获得了广泛的应用，用它可以验证电力系统接线方式和运行方式的合理性，计算输电线路的最大输送功率，确定系统故障切除的临界时间，以及研究某些提高电力系统稳定措施的效果，等等。
2. 仿真原理
对发电机、负荷及电力网络采用不同的数学模型可以构成各种简化的暂态稳定分析程序，采用何种组合要根据所研究问题的性质而定。本文采用的数学模型和计算方法如下所示。
| 内容 | 方法 |
|---|---|
| 发电机 | 发电机暂态电势Eq’保持恒定 |
| 负荷 | 较小的负荷用恒定阻抗 |
| 电力网络 | 导纳矩阵 |
| 微分方程 | 改进欧拉法 |
| 网络方程 | 直接法 |
3. 流程
流程图\label{flow}
3.1 初值计算
主要求取进行稳态分析前确定微分方程求解所需的初值。包括当前的虚拟电势$E_{q0}$、原动机功率$P_{m0}$、发电机转子角度初值$\delta_{0}$、发电机转子角速度初值$\omega_0$、负荷的等值导纳。
3.2 用直接法求解网络方程
此步骤中主要是求解网络潮流，发电机以不计阻尼绕组的模型并入电力网络，同时负荷或者故障也以等值导纳形式并入网络。利用发电机此时的$\delta$值确定参数$b、g、G_{xi}、B_{yi}、G_{yi}、$。其中$G_{xi}、B_{yi}、G_{yi}、$是发嗲你就并入网络的导纳，$b、g$参数用于计算发电机节点的虚拟注入电流。在网络矩阵生成、发电机虚拟注入电流得出后，利用$YV=I$，直接求出网络方程，得出节点电压后算出发电机节点注入电流。
3.3 改进欧拉法求解微分方程
转子运动方程模型如下：
$$\begin{cases}
\frac{d\delta_i}{dt}=\omega_s(\omega_i-1)\cr
\frac{d\omega_i}{dt}=\frac 1{T_{Ji}}(P_{mi}-P_{ei})
\end{cases}$$
根据求取初值的方法用t时刻各发电机的$\delta_i(t)$求取系统所有节点的电压$V_x(t)$、$V_y(t)$和发电机节点的注入电流$I_{xi}(t)$、$I_{yi}(t)$。根据转子运动方程求得$\delta_i$、$\omega_i$在t时刻的导数值，从而求得$t+\Delta{t}$时刻的状态估计量。根据状态估计量重新按求取初值的方法计算出所有系统的节点电压$V_x(t+\Delta{t})$、$V_y(t+\Delta{t})$和发电机注入电流$I_{xi}(t+\Delta{t})$、$I_{yi}(t+\Delta{t})$。用同样的方法根据转子运动方程计算求得$t+\Delta{t}$时刻$\delta_{i}$、$\omega_{i}$的导数值。最后，将$t$时刻的导数值和$t+\Delta{t}$时刻的导数值的平均值（即相加除以2）与$\Delta{t}$的乘积和$t$时刻的$\delta_{i}$、$\omega_{i}$相加即得到了$t+\Delta{t}$时刻$\delta_{i}$、$\omega_{i}$值。
4. 程序结果
相对摇摆角和时间的关系曲线\label{flow}
5. 源码分析
PSTS选择界面主要进行了程序的可视化输出，让用户选择需要进行的计算模式。具体模式分析如下：
考虑相同切除时间，是否计及突极效应，在模式1下计算了在$ct=0.08333s$时候考虑和不考虑突极效应的情况。
考虑同样情况，不同切除时间，在模式2下计算了不计突极效应$ct=0.163s$和$ct=0.162s$的两种情况。
考虑同样情况，不同切除时间，在模式3计算了计及突极效应$ct=0.086s$和$ct=0.085s$的两种情况。
源码如下：
disp('选择模式：\n');
disp('1:相对摇摆角与时间的关系曲线\n');
disp('2:临界切除时间周围的相对摇摆角与时间的关系曲线（不计凸极效应）\n');
disp('3:临界切除时间周围的相对摇摆角与时间的关系曲线（计及凸极效应）\n');
prompt='请选择：\n';
a = input(prompt);
if a== 1
disp('相对摇摆角与时间的关系曲线');
doPlot(1);
elseif a==2
disp('临界切除时间周围的相对摇摆角与时间的关系曲线（不计凸极效应）');
doPlot(2);
elseif a==3
disp('临界切除时间周围的相对摇摆角与时间的关系曲线（计及凸极效应）');
doPlot(3);
else
disp('输入错误,请重新输入！\n');
end
主函数主要进行了数据分析和绘图。
源码如下：
%% doPlot: function description
function [] = doPlot(arg)
if arg==2
% fClear=0.08333;%故障切除时刻
% fClear=0.08600;%计及凸极效应临界故障切除时刻
fClear=[0.16300 0.16200];%不计凸极效应临界故障切除时刻
temp = fClear;
spEffect=0;%0为不计凸极效应
elseif arg==3
fClear=[0.086 0.085];%不计凸极效应临界故障切除时刻
temp = fClear;
spEffect=1;%1为计凸极效应
elseif arg==1
fClear=0.08333;
spEffect=0;
temp = spEffect;
end
dt=0.0005;%步长
t_total=2;%总时间
global freq
freq=60;%系统频率
%%已知数据
%发电机数据
global Tj Ra Xd X_d Xq X_q T_d0 D P_mech
Tj=[47.28 12.80 6.02]';
Ra=[0 0 0]';
Xd=[0.146 0.8958 1.3125]';
X_d=[0.0608 0.1198 0.1813]';
Xq=[0.0969 0.8645 1.2578]';
X_q=[0.0969 0.1969 0.25]';
T_d0=[8.96 6.00 5.89]';
D=[0 0 0]';
%正常运行情况下的系统潮流
V_amp=[1.04 1.025 1.025 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324]';%节点电压幅值
V_phase=deg2rad([0 9.28 4.6648 -2.2168 -3.9888 -3.6874 3.7197 0.7275 1.9667])';%节点电压相角，弧度值
P_gen=[0.7164 1.63 0.85 0 0 0 0 0 0]';%发电机注入有功功率
Q_gen=[0.2705 0.0665 -0.1086 0 0 0 0 0 0]';%发电机注入无功功率
P_load=[0 0 0 0 1.25 0.9 0 1 0]';%节点负荷有功功率
Q_load=[0 0 0 0 0.5 0.3 0 0.35 0]';%节点负荷无功功率
gen_index=find(P_gen);%标示发电机节点
load_index=find(P_load);%标示负荷节点
deltaDiff=zeros(ceil(t_total/dt),2);%相对摇摆角定义
for temp=0:1%0为不计凸极效应，1为计及凸极效应
w=[1 1 1]'; %角速度初值
V=V_amp.*exp(j*V_phase);%节点电压
S_gen=P_gen(gen_index)+j*Q_gen(gen_index);%发电机注入复功率
S_load=P_load(load_index)+j*Q_load(load_index);%负荷复功率
V_gen=V(gen_index);%发电机电压
V_load=V(load_index);%负荷电压
V_amp_gen=V_amp(gen_index);%发电机电压幅值
V_amp_load=V_amp(load_index);%负荷电压幅值
%%初值计算
%发电机数据
I_gen=conj(S_gen)./conj(V_gen);%
P_mech=P_gen(gen_index)+(real(I_gen).^2+imag(I_gen).^2).*Ra;
if spEffect==1%计及凸极效应
Xq=[0.0969 0.8645 1.2578]';
elseif spEffect==0%不计凸极效应
Xq=X_d;
end
EQ=V_gen+(Ra+j*Xq).*I_gen;
delta=angle(EQ);%初始功角
Vq=cos(delta).*real(V_gen)+sin(delta).*imag(V_gen);
Iq=cos(delta).*real(I_gen)+sin(delta).*imag(I_gen);
Id=sin(delta).*real(I_gen)-cos(delta).*imag(I_gen);
E_q=Vq+Ra.*Iq+X_d.*Id;%初始暂态电势
%负荷等值并联导纳
Y_load=conj(S_load)./(V_amp_load.^2);
%%故障系统与故障后系统描述
%线路原始导纳矩阵
Yn=zeros(9,9);
Yn(5,4)=-1/(0.010+j*0.085);Yn(4,5)=Yn(5,4);
Yn(6,4)=-1/(0.017+j*0.092);Yn(4,6)=Yn(6,4);
Yn(7,5)=-1/(0.032+j*0.161);Yn(5,7)=Yn(7,5);
Yn(9,6)=-1/(0.039+j*0.170);Yn(6,9)=Yn(9,6);
Yn(8,7)=-1/(0.0085+j*0.072);Yn(7,8)=Yn(8,7);
Yn(9,8)=-1/(0.0119+j*0.1008);Yn(8,9)=Yn(9,8);
Yn(4,1)=-1/(j*0.0576);Yn(1,4)=Yn(4,1);
Yn(7,2)=-1/(j*0.0625);Yn(2,7)=Yn(7,2);
Yn(9,3)=-1/(j*0.0586);Yn(3,9)=Yn(9,3);
Yn(1,1)=-Yn(1,4);
Yn(2,2)=-Yn(2,7);
Yn(3,3)=-Yn(3,9);
Yn(4,4)=-Yn(4,5)-Yn(4,6)-Yn(1,4)+j*0.088+j*0.079;
Yn(5,5)=-Yn(5,4)-Yn(5,7)+j*0.088+j*0.153;
Yn(6,6)=-Yn(6,4)-Yn(6,9)+j*0.079+j*0.179;
Yn(7,7)=-Yn(7,8)-Yn(2,7)-Yn(5,7)+j*0.153+j*0.0745;
Yn(8,8)=-Yn(7,8)-Yn(8,9)+j*0.0745+j*0.1045;
Yn(9,9)=-Yn(6,9)-Yn(8,9)-Yn(3,9)+j*0.179+j*0.1045;
%加入负荷之后的节点导纳矩阵
Yn(load_index(1),load_index(1))=Yn(load_index(1),load_index(1))+Y_load(1);
Yn(load_index(2),load_index(2))=Yn(load_index(2),load_index(2))+Y_load(2);
Yn(load_index(3),load_index(3))=Yn(load_index(3),load_index(3))+Y_load(3);
%将故障并入网络
Yn_fault=Yn;
Yn_fault(7,7)=10^10;
Yn_real_fault=yConvert(Yn_fault);%得到实数形式节点导纳矩阵
%将发电机导纳并入矩阵
Yn_real_fault_gen=genMerge(Yn_real_fault,delta);
[VxVy,IxIy]=getPF(Yn_real_fault_gen,E_q,delta);
%故障后系统求解
n_total=ceil(t_total/dt);
n_fault=ceil(fClear(fcm+1)/dt);
deltaOutput = zeros(n_total,3);
for li=1:n_fault
[delta1,w1]=doEuler(VxVy,IxIy,Yn_real_fault,delta,w,dt,E_q);
Yn_real_fault_gen=genMerge(Yn_real_fault,delta1);
[VxVy,IxIy]=getPF(Yn_real_fault_gen,E_q,delta1);
deltaOutput(li,:)=rad2deg(delta1');
delta=delta1;
w=w1;
end
%故障切除后系统求解
%生成切除故障后的节点导纳矩阵
Yn57=zeros(9,9);
Yn57(5,5)=1/(0.032+j*0.161)+j*0.153;
Yn57(7,7)=1/(0.032+j*0.161)+j*0.153;
Yn57(5,7)=-1/(0.032+j*0.161);
Yn57(7,5)=-1/(0.032+j*0.161);
Yn_faultclear=Yn-Yn57;%9阶矩阵
Yn_real_faultclear=yConvert(Yn_faultclear);
%将发电机导纳并入矩阵
Yn_real_faultclear_gen=genMerge(Yn_real_faultclear,delta);
[VxVy,IxIy]=getPF(Yn_real_faultclear_gen,E_q,delta);
for li=(n_fault+1):n_total
[delta1,w1]=doEuler(VxVy,IxIy,Yn_real_faultclear,delta,w,dt,E_q);
Yn_real_faultclear_gen=genMerge(Yn_real_faultclear,delta1);
[VxVy,IxIy]=getPF(Yn_real_faultclear_gen,E_q,delta1);
deltaOutput(li,:)=rad2deg(delta1');
delta=delta1;
w=w1;
end
%相对摇摆角计算
deltaDiff(:,fcm+1)=deltaOutput(:,2)-deltaOutput(:,1);
end
t=dt:dt:n_total*dt;
plot(t,deltaDiff(:,1),'--',t,deltaDiff(:,2));
%ylim([0 300]);
set(gca,'XTick',0:0.1:t_total);
grid on;
xlabel('时间t/s');
ylabel('相对摇摆角/°');
if arg==2
title('临界切除时间周围的相对摇摆角与时间的关系曲线（不计凸极效应）');
elseif arg==3
title('临界切除时间周围的相对摇摆角与时间的关系曲线（计及凸极效应）');
elseif arg==1
title('相对摇摆角与时间的关系曲线');
end
主函数中调用了一些功能函数，主要实现以下功能：
- 直接法求网络方程
源码如下：
function [VxVy,IxIy]=getPF(Yn_real,E_q,delta)
[m,n]=size(Yn_real);
I=zeros(m,1);
for li=1:3
[bg ,GB]=getGen(li,delta(li));
I(li*2-1,1)=E_q(li)*bg(1);
I(li*2,1)=E_q(li)*bg(2);
end
VxVy=Yn_real\I;
IxIy=zeros(m,1);
for li=1:3
[bg ,GB]=getGen(li,delta(li));
IxIy(li*2-1,1)=E_q(li)*bg(1)-GB(1,1)*VxVy(li*2-1)-GB(1,2)*VxVy(li*2);
IxIy(li*2,1)=E_q(li)*bg(2)-GB(2,1)*VxVy(li*2-1)-GB(2,2)*VxVy(li*2);
end
end
- 发电机并入导纳矩阵参数和虚拟注入电流参数计算
源码如下：
function [ bg , GB ] = getGen( number , delta )
global Ra X_d Xq %发电机参数
GB=zeros(2,2);
GB(1,1)=(Ra(number)-(X_d(number)-Xq(number))*sin(delta)*cos(delta))/(Ra(number)^2+X_d(number)*Xq(number));
GB(1,2)=(X_d(number)*cos(delta)^2+Xq(number)*sin(delta)^2)/(Ra(number)^2+X_d(number)*Xq(number));
GB(2,1)=(-X_d(number)*sin(delta)^2-Xq(number)*cos(delta)^2)/(Ra(number)^2+X_d(number)*Xq(number));
GB(2,2)=(Ra(number)+(X_d(number)-Xq(number))*sin(delta)*cos(delta))/(Ra(number)^2+X_d(number)*Xq(number));
bg=zeros(1,2);
bg(1)=(Ra(number)*cos(delta)+Xq(number)*sin(delta))/(Ra(number)^2+X_d(number)*Xq(number));
bg(2)=(Ra(number)*sin(delta)-Xq(number)*cos(delta))/(Ra(number)^2+X_d(number)*Xq(number));
end
- 导纳矩阵转化程序
源码如下：
%% function to convert the Yn(9*9) to Yn(18*18)
function [Yn18] = yConvert(Yn)
Yn18=zeros(18,18);
for li=1:1:9
for lj=1:9
ii=li*2-1;jj=lj*2-1;
Yn18(ii,jj)=real(Yn(li,lj));
Yn18(ii,jj+1)=-imag(Yn(li,lj));
Yn18(ii+1,jj)=imag(Yn(li,lj));
Yn18(ii+1,jj+1)=real(Yn(li,lj));
end
end
end
- 发电机等效参数并入导纳矩阵
源码如下：
function [ Yn_real_new ] = genMerge( Yn_real , delta )
Yn_real_new=Yn_real;
for li=1:3
[bg,GB]=getGen(li,delta(li));
ii=li*2-1;jj=li*2-1;
Yn_real_new(ii,jj)=Yn_real_new(ii,jj)+GB(1,1);
Yn_real_new(ii,jj+1)=Yn_real_new(ii,jj+1)+GB(1,2);
Yn_real_new(ii+1,jj)=Yn_real_new(ii+1,jj)+GB(2,1);
Yn_real_new(ii+1,jj+1)=Yn_real_new(ii+1,jj+1)+GB(2,2);
end
end
- 改进欧拉法计算程序
源码如下：
function [deltat1 wt1]=doEuler(VxVy,IxIy,Yn_real,deltat0,wt0,dt,E_q)
global Tj P_mech Ra freq
VgeneratorX=zeros(3,1);
VgeneratorX(1,1)=VxVy(1,1);
VgeneratorX(2,1)=VxVy(3,1);
VgeneratorX(3,1)=VxVy(5,1);
VgeneratorY=zeros(3,1);
VgeneratorY(1,1)=VxVy(2,1);
VgeneratorY(2,1)=VxVy(4,1);
VgeneratorY(3,1)=VxVy(6,1);
IgeneratorX=zeros(3,1);
IgeneratorX(1,1)=IxIy(1,1);
IgeneratorX(2,1)=IxIy(3,1);
IgeneratorX(3,1)=IxIy(5,1);
IgeneratorY=zeros(3,1);
IgeneratorY(1,1)=IxIy(2,1);
IgeneratorY(2,1)=IxIy(4,1);
IgeneratorY(3,1)=IxIy(6,1);
Peit0=VgeneratorX.*IgeneratorX+VgeneratorY.*IgeneratorY+(IgeneratorX.^2+IgeneratorY.^2).*Ra;
DerivativeDeltat0=2*pi*freq*(wt0-1);
DerivativeWt0=(P_mech-Peit0)./Tj;
deltat1_0=deltat0+DerivativeDeltat0.*dt;
wt1_0=wt0+DerivativeWt0.*dt;
Yn_generator = genMerge(Yn_real,deltat1_0);
[VxVyt1_0,IxIyt1_0]=getPF(Yn_generator,E_q,deltat1_0);
VgeneratorXt1_0=zeros(3,1);
VgeneratorXt1_0(1,1)=VxVyt1_0(1,1);
VgeneratorXt1_0(2,1)=VxVyt1_0(3,1);
VgeneratorXt1_0(3,1)=VxVyt1_0(5,1);
VgeneratorYt1_0=zeros(3,1);
VgeneratorYt1_0(1,1)=VxVyt1_0(2,1);
VgeneratorYt1_0(2,1)=VxVyt1_0(4,1);
VgeneratorYt1_0(3,1)=VxVyt1_0(6,1);
IgeneratorXt1_0=zeros(3,1);
IgeneratorXt1_0(1,1)=IxIyt1_0(1,1);
IgeneratorXt1_0(2,1)=IxIyt1_0(3,1);
IgeneratorXt1_0(3,1)=IxIyt1_0(5,1);
IgeneratorYt1_0=zeros(3,1);
IgeneratorYt1_0(1,1)=IxIyt1_0(2,1);
IgeneratorYt1_0(2,1)=IxIyt1_0(4,1);
IgeneratorYt1_0(3,1)=IxIyt1_0(6,1);
Peit1_0=VgeneratorXt1_0.*IgeneratorXt1_0+VgeneratorYt1_0.*IgeneratorYt1_0+(IgeneratorXt1_0.^2+IgeneratorYt1_0.^2).*Ra;
DerivativeDeltat1_0=2*pi*freq*(wt1_0-1);
DerivativeWt1_0=(P_mech-Peit1_0)./Tj;
deltat1=deltat0+(DerivativeDeltat0+DerivativeDeltat1_0).*(dt/2);
wt1=wt0+(DerivativeWt0+DerivativeWt1_0).*(dt/2);
end