IEEE39节点系统下分析各支路暂态稳定性指标（IEEE39 Transient stability analysis）

IEEE39节点系统下设置单点三相短路故障，运用能量函数法分析各支路暂态稳定性指标

function BPEM（）
%%%%%%%%%%%%%%%数据输入%%%%%%%%%%%%%%%%%%%%%
global SS S;
a=‘1‘
SS=xlsread（a‘4‘）;%系统可变参数 随故障类型或地点发生变化的数据
S=xlsread（a‘S‘）;%系统的不变参数
% [SBIvpb]=BPAM（SSS）;%支路势能稳定指标
% [vk_ratio]=index（SSS）;
% %%%%%%%%%%%%%%%%%数据输入%%%%%%%%%%%%%%%%%%%%%
% XX=xlsread（‘111111.xls‘‘XX‘）;%随故障发生变化的参数
% X=xlsread（‘111111.xls‘‘X‘）;%随故障发生变化的参数
% Pt1=cut_machine（XXX）;%确定切机量
%  function[SBIvpb]=BPEM（SSS）
%BPEM（Branch Potential Energy Method）
%function[SBIvpb]=PROJECT（nbGMSnSbnum_Tnum_LIJnum_GLstept）;
%program
% n=input（‘请输入节点的个数：n=‘）;
% b=input（‘请输入支路数：b=‘）;
% G=input（‘请输入发电机与所连母线的编号：G=‘）;
% M=input（‘请输入每台发电机的惯性时间常数：M=‘）;
% Sn=input（‘请输入每台发电机的额定容量：Sn=‘）;
% Sb=input（‘请输入系统的基准容量：Sb=‘）;
% num_T=input（‘请输入变压器的台数：num_T=‘）;
% num_L=input（‘请输入系统负荷的个数：num_L=‘）;
% I=input（‘请输入各支路的起始节点：I=‘）;
% J=input（‘请输入各支路的起始节点：J=‘）;
% num_G=input（‘请输入发电机的台数：num_G=‘）;
% L=input（‘请输入各负荷所在母线的编号：L=‘）;
% step=input（‘请输入仿真步长：step=‘）;
% t=input（‘请输入故障结束时刻：t=‘）;%
n=S（11）;
b=S（12）;
num_G=S（18）;
num_T=S（17）;
step=S（19）;
Sb=S（16）;
G=S（1:num_G3）;
M=S（1:num_G4）;
Sn=S（1:num_G5）;
I=S（1:b11）;
J=S（1:b12）;
t=SS（11+b+n+2*num_G+1）;
t0=SS（11+b+n+2*num_G+2）/step+1;
for i=1:b
    p（:i）=SS（:i+1）;
end
% w-各节点的角频率
for i=1:n
    w（:i）=SS（:i+b+1）;
end
% wg-各发电机的角频率;
for i=1:num_G
    wg（:i）=SS（:i+b+n+num_G+1）;
end
% pg-各发电机的有功出力
for i=1:num_G
    pg（:i）=SS（:i+1+b+n）;
end
[sy]=size（p）;
tc=t/0.01+1;
%***************************************************惯量中心频率
new_M=diag（M）*Sn/Sb;
for i=1:size（wg1）
 wcoi（i1）=wg（i:）* new_M/（sum（new_M））;
end
for i=1:n
new_w（:i）=w（:i）-wcoi（:1）;
end
for i=1:num_G
new_wg（:i）=wg（:i）-wcoi（:1）;
end
wn=314;
for i=1:b
    ww（:i）=w（:I（i））-w（:J（i））;
end
for i=1:num_G
    wgw（:i）=wg（:i）-w（:G（i））;
end
 %发电机内节点的暂态势能
vpn=zeros（snum_G）;           
for i=tc:s
    for j=1:num_G
        if i==tc
            vpn（ij）=wgw（ij）*（pg（ij）-pg（1j））*wn*step;
        else
            vpn（ij）=vpn（i-1j）+（wgw（ij）*（pg（ij）-pg（1j））+wgw（i-1j）*（pg（i-1j）-pg（1j）））/2*wn*step;
        end
    end
end 
%同步机动能
vk1=zeros（snum_G）;             
for i=tc:s
    for j=1:num_G
    vk1（ij）=1/2*（new_wg（ij））^2*new_M（j）*wn;
    end
end
 %支路的暂态势能
vpb=zeros（sb）;           
for i=tc:s
    for j=1:b
        if i==tc
            vpb（ij）=ww（ij）*（p（ij）-p（1j））*wn*step;
        else
            vpb（ij）=vpb（i-1j）+（ww（ij）*（p（ij）-p（1j））+ww（i-1j）*（p（i-1j）-p（1j）））/2*wn*step;
        end
    end
end 
 %****************支路的势能和
vps=zeros（s1）;             
for i=1:b
    vps（:1）=vps（:1）+vpb（:i）;
end
%***************发电机内节点势能和
vpns=zeros（s1）;
for i=1:num_G
    vpns（:1）=vpns（:1）+vpn（:i）;
end
%*****************系统总势能和
vp=zeros（s1）;
vp（:1）=vps（:1）+vpns（:1）;

vk=zeros（s1）;              
for i=1:num_G
    vk（:1）=vk（:1）+vk1（:i）;
end
v=vp+vk;
%*************************************************************
  i=0.01:0.01:5.01;
  figure（1）;
  plot

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件    3947520  2013-12-17 10:51  1.xls

     文件       5476  2014-02-28 20:17  BPEM.m

----------- ---------  ---------- -----  ----

              3952996                    2