%%%IEEE6算例初始潮流計算
n=6;
nl=7;
zn=[0?0?0?0?0?0];
ph=[3?4?2?3?1?2?4];??
qh=[5?6?5?4?3?6?1];
zz=[0.133i?0.3i?0.723+1.05i?0.097+0.407i?0.08+0.37i?0.282+0.64i?0.123+0.518i];
zx=[0?0?0?0?0?0?0];
gd=[1?1?0?0?0?0?0];
sg=[0?0.501?0?0?0?0];
sl=[0?0?0?0.5+0.05i?0.55+0.13i?0.3+0.18i];
jl=[1?3?2?2?2?2];
pr=0.00001;
ijl=[1?3?2?2?2?2];
uo=[1.05?1.1?1?1?1?1];%1.05?1.1????1.1?1.025????
kk=[1.1??1.025?1?1?1?1?1];?????????????????????????????????
qb=[0?0?0?0?0?0];????????????????????????
y=zeros（n）;
?for?i=1:n
????if?zn（i）~=0;
????????y（ii）=1./zn（i）;
????end
?end
?for?ji=1:nl
????if?gd（ji）==0???
????????p=ph（ji）;q=qh（ji）;%某支路首端號p某支路末端號q???
????else?p=qh（ji）;q=ph（ji）;%反之
????end
????y（pq）=?y（pq）-1./（zz（ji）*kk（ji））;
????y（qp）=?y（pq）;
????y（qq）=?y（qq）+1./（zz（ji）*kk（ji）^2）+zx（ji）;
????y（pp）=?y（pp）+1./zz（ji）+zx（ji）;???
?end
b=imag（y）;
gg=real（y）;bb=imag（y）;
bj=zeros（n-1）;
for?i=1:n??
????if?jl（i）==1??????
?????????b（i:）=[];??
?????????b（:i）=[];?????????????????
????end????
end
bj=b;
jl（jl==1）=[];
ind=find（jl==2）;
?w=length（ind）;
for?i=1:n-1?????
????if?jl（i）==3;?????
????????b（i:）=[];
????????b（:i）=[];
????end
end
bjj=zeros（w）;
bjj=b;
%disp（bjj）
for?i=1:n
???s（i）=sg（i）-sl（i）+qb（i）;
end
p=real（s）;q=imag（s）;
v=zeros（1n）;o=zeros（1n）;
for?i=1:n?
????e（i）=real（uo（i））;
????f（i）=imag（uo（i））;
????v（i）=uo（i）;?
????o（i）=atan（f（i）/e（i））;
end????
????
?jd1=1;jd2=1;dt=0;kp=1;kq=1;
????while??jd1~=0|jd2~=0;
???????????jd1=0;jd2=0;?
???????????dp=zeros（size（ijl））;dpl=zeros（size（ijl））;
???????for?i=1:n
???????????if??ijl（i）~=1
??????????????c（i）=0;
???????????????for?k=1:n
????????????????c（i）=c（i）+v（k）*（gg（ik）*cos（（o（i）-o（k））*pi./180）+bb（ik）*sin（（o（i）-o（k））*pi./180））;?????????????
???????????????end??????????
??????????????dpl（i）=p（i）-v（i）*c（i）;
??????????????dp（i）=dpl（i）./?v（i）;
???????????end
???????end?
????????????dp（ijl==1）=[];
????????????dpl（ijl==1）=[];???????????????
????????????for?i=1:n-1
????????????detp（i）=abs（dpl（i））;?
????????????end
??????????if??detp>=pr
???????????????jd1=jd1+1;???
??????????end??????
????if??jd1~=0??????????
??????????kq=1;???????
?????????dto=-dp*inv（bj）*180./pi;
?????????po=zeros（size（ijl））;
?????????po（ijl~=1）=dto;??????????
??????for?i=1:n
????????deto（i）=po（i）./v（i）;
????????o（i）=o（i）+deto（i）;?
??????end
??????dq=zeros（size（ijl））;dql=zeros（size（ijl））;
??????for?i=1:n
???????if?ijl（i）==2
????????????d（i）=0;
?????????????for?k=1:n
???????????????????d（i）=d（i）+v（k）*（gg（ik）*sin（（o（i