%----------節點導納矩陣程序lfybus.m----------
tic???%?計時開始
i=sqrt（-1）;
j=sqrt（-1）;
nl=linedata（:1）;?
nr=linedata（:2）;?
R=linedata（:3）;
X=linedata（:4）;?
Bc=j*linedata（:5）;?
a=linedata（:?6）;
nbr=length（linedata（:1））;?
nbus=max（max（nl）?max（nr））;
Z=R?+?j*X;?
y=ones（nbr1）./Z;?%支路導納
for?n=1:nbr
????if?a（n）<=0??
???????a（n）=1;?
????else?end
????Ybus=zeros（nbusnbus）;?%?將Ybus初始化為0
%?---非對角元素的形成---
????for?k=1:nbr;
????????Ybus（nl（k）nr（k））=Ybus（nl（k）nr（k））-y（k）/a（k）;
????????Ybus（nr（k）nl（k））=Ybus（nl（k）nr（k））;
????end
end
%---對角元素的形成---
for??n=1:nbus
?????for?k=1:nbr
?????????if?nl（k）==n
????????????Ybus（nn）=Ybus（nn）+y（k）/（a（k）^2）?+?Bc（k）;
?????????elseif?nr（k）==n
????????????????Ybus（nn）=Ybus（nn）+y（k）?+Bc（k）;
?????????else?end
?????end
end


%----------PQ分解法潮流程序decouple.m----------
%1.快速解耦法電力系統潮流計算（PQ分解法）
%2.數據準備;用戶需要在命令窗口編寫busdata和linedata兩個數據文件
basemva=100;??
accuracy=0.0001;?
maxiter=50;
ns=0;?
Vm=0;?
delta=0;?
yload=0;?
deltad=0;
nbus?=?length（busdata（:1））;
kb=[];
Vm=[];?
delta=[];?
Pd=[];?
Qd=[];?
Pg=[];?
Qg=[];?
Qmin=[];?
Qmax=[];
Pk=[];?
P=[];?
Qk=[];?
Q=[];?
S=[];?
V=[];
for?k=1:nbus
????n=busdata（k1）;
????kb（n）=busdata（k2）;?
????Vm（n）=busdata（k3）;?
????delta（n）=busdata（k4）;
????Pd（n）=busdata（k5）;?
????Qd（n）=busdata（k6）;?
????Pg（n）=busdata（k7）;?
????Qg（n）?=?busdata（k8）;
????Qmin（n）=busdata（k?9）;?
????Qmax（n）=busdata（k10）;
????Qsh（n）=busdata（k11）;
????if?Vm（n）?<=?0??
????????Vm（n）?=?1.0;?
????????V（n）?=?1?+?j*0;
????else?delta（n）?=?pi/180*delta（n）;
?????????V（n）?=?Vm（n）*（cos（delta（n））?+?j*sin（delta（n）））;
?????????P（n）=（Pg（n）-Pd（n））/basemva;
?????????Q（n）=（Qg（n）-Qd（n）+?Qsh（n））/basemva;
?????????S（n）?=?P（n）?+?j*Q（n）;
????end
?if?kb（n）?==?1
??????ns?=?ns+1;?
?else?end
nss（n）?=?ns;
end
Ym?=?abs（Ybus）;?
t?=?angle（Ybus）;
ii=0;
for?ib=1:nbus
?????if?kb（ib）?==?0?|?kb（ib）?==?2
?????ii?=?ii+1;
?????jj=0;
?????????for?jb=1:nbus
?????????????if?kb（jb）?==?0?|?kb（jb）?==?2
????????????????jj?=?jj+1;
????????????????B1（iijj）=imag（Ybus（ibjb））;
?????????????elseend
?????????end
?????else?end
end

ii=0;
for?ib=1:nbus
????if?kb（ib）?==?0
???????ii?=?ii+1;
???????jj=0;
???????for?jb=1:nbus
????????????if?kb（jb）?==?0
???????????????jj?=?jj+1;
???????????????B2（iijj）=imag（Ybus（ibjb））;
?????????????elseend
???????end
????else?end
end
B1inv=inv（B1）;?
B2inv?=?inv（B2）;
maxerror?=?1;?
converge?=?1;?
iter?=?0;
mline=ones（nbr1）;
for?k=1:nbr
????for?m=k+1:nbr
????????if（（nl（k）==nl（m））?&?（nr（k）==nr（m）））;
???????????mline（m）=2;
?????????elseif?（（nl（k）==nr（m））?&?（nr（k）==nl（m）））;
?????????????????mline（m）=2;
?????????else?end
????end
end

%?decouple程序開始迭代
while?maxerror?>=?accuracy?&?iter?<=?maxiter???%檢驗不平衡功率
??????iter?=?iter+1;
??????id=0;?
??????iv=0;
for?n=1:nbus
????nn=n-nss（n）;
????J11=0;??
????J33=0;
	for?ii=1:nbr
		if?mline（ii）==1???
?????		if?nl（ii）?==?n?|?nr（ii）?==?n
????????	????if?nl（ii）?==?n?