A=[?1	2	0.01938	0.05917	0.0528	0	0	0	0	0	1	-360	360;
	1	5	0.05403	0.22304	0.0492	0	0	0	0	0	1	-360	360;
	2	3	0.04699	0.19797	0.0438	0	0	0	0	0	1	-360	360;
	2	4	0.05811	0.17632	0.034	0	0	0	0	0	1	-360	360;
	2	5	0.05695	0.17388	0.0346	0	0	0	0	0	1	-360	360;
	3	4	0.06701	0.17103	0.0128	0	0	0	0	0	1	-360	360;
	4	5	0.01335	0.04211	0	0	0	0	0	0	1	-360	360;
	4	7	0	0.20912	0	0	0	0	0.978	0	1	-360	360;
	4	9	0	0.55618	0	0	0	0	0.969	0	1	-360	360;
	5	6	0	0.25202	0	0	0	0	0.932	0	1	-360	360;
	6	11	0.09498	0.1989	0	0	0	0	0	0	1	-360	360;
	6	12	0.12291	0.25581	0	0	0	0	0	0	1	-360	360;
	6	13	0.06615	0.13027	0	0	0	0	0	0	1	-360	360;
	7	8	0	0.17615	0	0	0	0	0	0	1	-360	360;
	7	9	0	0.11001	0	0	0	0	0	0	1	-360	360;
	9	10	0.03181	0.0845	0	0	0	0	0	0	1	-360	360;
	9	14	0.12711	0.27038	0	0	0	0	0	0	1	-360	360;
	10	11	0.08205	0.19207	0	0	0	0	0	0	1	-360	360;
	12	13	0.22092	0.19988	0	0	0	0	0	0	1	-360	360;
	13	14	0.17093	0.34802	0	0	0	0	0	0	1	-360	360;
];
Bus=[
	1	3	0	0	0	0	1	1.06	0	0	1	1.06	0.94;
	2	2	21.7	12.7	0	0	1	1.045	-4.98	0	1	1.06	0.94;
	3	2	94.2	19	0	0	1	1.01	-12.72	0	1	1.06	0.94;
	4	1	47.8	-3.9	0	0	1	1.019	-10.33	0	1	1.06	0.94;
	5	1	7.6	1.6	0	0	1	1.02	-8.78	0	1	1.06	0.94;
	6	2	11.2	7.5	0	0	1	1.07	-14.22	0	1	1.06	0.94;
	7	1	0	0	0	0	1	1.062	-13.37	0	1	1.06	0.94;
	8	2	0	0	0	0	1	1.09	-13.36	0	1	1.06	0.94;
	9	1	29.5	16.6	0	19	1	1.056	-14.94	0	1	1.06	0.94;
	10	1	9	5.8	0	0	1	1.051	-15.1	0	1	1.06	0.94;
	11	1	3.5	1.8	0	0	1	1.057	-14.79	0	1	1.06	0.94;
	12	1	6.1	1.6	0	0	1	1.055	-15.07	0	1	1.06	0.94;
	13	1	13.5	5.8	0	0	1	1.05	-15.16	0	1	1.06	0.94;
	14	1	14.9	5	0	0	1	1.036	-16.04	0	1	1.06	0.94;
];
n=size（Bus1）;%求系統節點數
m=size（A1）;%確定支路數
%%
%做一些基礎工作
%按PQ?PV排列的節點信息矩陣
Bus2=sortrows（Bus2）;
list=zeros（n1）;
?for?i=1:n;
?????a=Bus2（i1）;
?????list（i）=a;
?end
?PQn=length（find（Bus2（:2）==1））;%PQ節點的個數
?PVn=length（find（Bus2（:2）==2））;%PV節點的個數
?%求標幺值
?sb=100;
?for?i=1:n
?????Bus2（i3）=Bus2（i3）/sb;
?????Bus2（i4）=Bus2（i4）/sb;
?end
%%
%節點導納矩陣的求解
Y=zeros（n）;%形成一個n階的零矩陣
M=zeros（n）;%形成一個n階的零矩陣用于儲存自導的臨時結果
K=zeros（n）;%形成一個n階的零矩陣用于儲存自導的臨時結果
for?i=1:m;
????Y（A（i1）A（i2））=（A（i3）+j*A（i4））^（-1）;
????a=A（i1）;
????b=A（i2）;
????r=A（i3）;
????x=A（i4）;
????Y（ab）=-1/（r+j*x）;
????Y（ba）=Y（ab）;
end%求互導
for?i=1:m;
??a=A（i1）;
??M（aa）=j*A（i5）/2+M（aa）;
end%掃描第一列求對地導納
for?i=1:m;
??b=A（i2）;
??K（bb）=A（b5）*j/2+K（bb）;
end%掃描第二列求對地導納
for?i=1:n;
Y（ii）=-sum（Y（i:））+M（ii）+K（ii）;
end%求每個節點的自導
G=real（Y）;%電導矩陣
B=imag（Y）;%電納矩陣
G=G（list:）;%按照PQ、PV順序的電導矩陣
B=B（list:）;%按照PQ、PV順序的電納矩陣
%%
%迭代開始
COUNT=0;
Tmax=100;%迭代次數和最大迭代次數。
V=10^（-6）;%設置精度。
%設電壓初始值
e=ones（n1）;
f=zeros（n1）;
%求不平衡量dp、dq、dv
Pi=zeros（n1）;
Qi=zeros（n1）;
Vi=zeros（n1）;
for?i=1:Tmax;
for?i=1:n;
????Pi（i1）=Bus2（i3