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