暂态稳定性分析大程

暂态稳定性分析大程

作者 : 李志
学号 : 21410076
Tel : 13732254786

目录

[TOC]

1. 引言

对于地区性的电力系统来说,一般失去暂态稳定过程发展快速,通常分析系统遭到干扰后先是摆摆周期(1~1.5s)的机电暂态过程就得看清系是否会保障平稳运转。这种状况下的暂态稳定分析着,由于调速系统的惯性,使得在短缺日内原动机的功率不会见有很老变迁,因此好忽略调速系统的用意而而原动机的功率保持不更换;此外由于发电机励磁绕组的时间常数较生,这样于缺乏日内该磁链也无会见产生强烈扭转,而针对励磁调节系统的图,可以用发电机暂态电势$E_q{’}$或$E{’}$保持稳来类模拟,即当在第一摇摆周期内,励磁绕组中自由电流分量的转移是因为励磁调节体系的来意所加,从而使鼓励磁绕组的磁链$Ψ_f$在即时段时内保障不转换。相应地,阻尼绕组的影响将略微去不计算。

简模型下之暂态稳定分析程序在电力系统运行及计划被获了科普的应用,用其好印证电力系统接线方式及运作方式的合理性,计算输电线路的极度老输送功率,确定系故障切除的逼时间,以及研究一些提高电力系统稳定措施之效能,等等。

2. 虚原理

针对发电机、负荷和电力网络利用不同之数学模型可以构成各种简化的暂态稳定分析程序,采用何种组合而基于所研究问题的性而肯定。本文采用的数学模型和计算方法如下所示。

内容 方法
发电机 发电机暂态电势Eq’保持恒定
负荷 较小的负荷用恒定阻抗
电力网络 导纳矩阵
微分方程 改进欧拉法
网络方程 直接法

3. 流程

亚洲城误乐城ca88网站 1

流程图\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. 程序结果

亚洲城误乐城ca88网站 2

对立摇摆角和岁月之关联曲线\label{flow}

5. 源码分析

PSTS选择界面主要进行了序的可视化输出,让用户选择得展开的计算模式。具体模式分析如下:

  1. 设想同切除时间,是否计及突极效应,在模式1下蛋计算了在$ct=0.08333s$时候考虑同莫考虑突极效应的情景。

  2. 考虑同情形,不同切除时间,在模式2下计算了非划算突极效应$ct=0.163s$和$ct=0.162s$的星星点点种情形。

  3. 设想同情况,不同切除时间,在模式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
  • 导纳矩阵转化程序

源码如下亚洲城误乐城ca88网站:

%% 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