matlab 6节点潮流计算程序

%%%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