%%?運(yùn)用牛頓-拉夫遜法解IEEE潮流??
%%%%%%把原始文件分割成兩個獨立的純數(shù)據(jù)文件分別保存為Bus.DAT和Line.DAT
?
clear;
clc;
tic;????????????????????????%開始計時
%%?讀取Bus.DAT中的數(shù)據(jù)
load（‘C:\Users\xyphhunj\Desktop\普通潮流\030bus\Bus.DAT‘）;
Node=Bus（:1）?;?????????????%設(shè)定節(jié)點號
Voltbase=Bus（:3）?;?????????%電壓標(biāo)準(zhǔn)值
Type1=Bus（:6）;?????????????%節(jié)點類型
V=Bus（:7）;?????????????????%節(jié)點電壓表幺值
Angle=Bus（:8）;?????????????%電壓相角
Pload=Bus（:9）;?????????????%負(fù)荷有功
Qload=Bus（:10）;????????????%負(fù)荷載無功
Pgen=Bus（:11）;?????????????%發(fā)電有功
Qgen=Bus（:12）;?????????????%發(fā)電無功
Vset=Bus（:14）;?????????????%電壓初始值
Qgmax=Bus（:15）;????????????%Q的最大值
Qgmin=Bus（:16）;????????????%Q的最小值
Qsh=Bus（:18）;??????????????%并聯(lián)電容電納為正，并聯(lián)電抗電納為負(fù)
base=100;???????????????????%功率標(biāo)準(zhǔn)值
%%??讀取Line.DAT中的數(shù)據(jù)
load（‘C:\Users\xyphhunj\Desktop\普通潮流\030bus\Line.DAT‘）;
II=Line（:1）;??????????????????
JJ=Line（:2）;??????????????????%兩端點
Type2=Line（:5）;???????????????%線路類型
R=Line（:6）;???????????????????%兩點間電阻
X=Line（:7）;???????????????????%兩點間電抗
B0=Line（:8）/2;????????????????%線路兩端電納
K=Line（:14）;??????????????????%?變壓器變比
Kmax=Line（:16）;
Kmin=Line（:17）;?
Kset=Line（:18）;???????????????%變比要求
%%?開始求導(dǎo)納?Y
num=size（Bus1）;
y0=zeros（num）?;??????????????????%兩點間的支路導(dǎo)納
yt=zeros（num）?;??????????????????%變壓器的對地導(dǎo)納
for?x=1:size（Line1）?
????switch?Type2（x）?????????????%比較線路類型??
????????case?2??????????????????%含變壓器的
????????????y0（II（x）JJ（x））=1./（K（x）.*（R（x）+1i*X（x）））+y0（II（x）JJ（x））?;???????????????????
????????????y0（JJ（x）II（x））=y0（II（x）JJ（x））????;???????????????
????????????yt（II（x）JJ（x））=（1-K（x））./（（R（x）+1i*X（x））.*K（x）^2）+yt（II（x）JJ（x））;?
????????????yt（JJ（x）II（x））=（K（x）-1）./（（R（x）+1i*X（x））.*K（x））+yt（JJ（x）II（x））;????
????????case?1?????????????????%不含變壓器的
????????????y0（II（x）JJ（x））=1/（R（x）+1i*X（x））+y0（II（x）JJ（x））;?????????????????
????????????y0（JJ（x）II（x））=y0（II（x）JJ（x））;??????????????????
????end
end
Y=zeros（num）;??????????????????%導(dǎo)納矩陣
clear?x
for?x=1:num
????Y（xx）=sum（y0（x:）+yt（x:））+1i*Qsh（x）?;?????%求出對角導(dǎo)納（不含B/2）
end
Y=Y-y0;
clear?x
for?x=1:size（Line1）
????Y（II（x）II（x））=Y（II（x）II（x））+1i*B0（x）?;
????Y（JJ（x）JJ（x））=Y（JJ（x）JJ（x））+1i*B0（x）?;??????????%加B
end
clear?x???????????????????????????%?Y導(dǎo)納矩陣結(jié)束
G=real（Y）;
B=imag（Y）;
%%?
D=ones（2*num1）;
%給定節(jié)點電壓值，a為電壓相角，u為電壓幅值
for?m=1:num
????a（m）=0;??????????%初始設(shè)定角度
????u（m）=Bus（m14）;??%初始迭代電壓值
????e（m）=u（m）*cos（a（m））;
????f（m）=u（m）*sin（a（m））;??%電壓實部和虛部
end
t=0;????????????????????%設(shè)定迭代次數(shù)變量t
while?（max（abs（D））>1e-8）

???for?m=1:num
??????Ps（m1）=（Pgen（m1）-Pload（m1））/base;
??????Qs（m1）=（Qgen（m1）-Qload（m1））/base;??????%節(jié)點注入有功和注入無功????????????????????????????????????????????????????????????????????????????
???end
%計算功率不平衡量有功不平衡量為dP，無功不平衡量為dQ電壓平方的不平衡量為dUU
???sum1=zeros（1num）;
???sum2=zeros（1num）;
???sum3=zeros（1num）;
???for?m=1:num
???????for?n=1:num
??????????sum1（m）=sum1（m）+G（mn）*e（n）-B（mn）*f（n）;
??????????sum2（m）=sum2（m）+G（mn）*f（n）+B（mn）*e（n）;?????
???????end
???????sum3（m）=e（m）^2+f（m）^2;
???end
???for?m=1:nu