33节点潮流计算matlab

%形成邻接矩阵
data=textread（‘IEEE33.dat‘）; 
a=zeros（3333）;
for j=1:32;
  a（data（j2）+1data（j3）+1）=1;
end
%%电压赋初值
v=ones（331）;
v（1）=1.05;
delt=zeros（331）;
V=v.*cos（delt）+i*v.*sin（delt）;
%%节点负荷功率
SL=zeros（331）;
for k=2:33
   SL（k）=（data（k-16）+i*data（k-17））/10;
end
SL（1）=0;
%%线路阻抗
Z=zeros（331）;
for m=2:33
    Z（m）=（data（m-14）+i*data（m-15））/100;
end
Z（1）=0;
S=zeros（331）;
%迭代
for dd=1:100
     V=v.*cos（delt）+i*v.*sin（delt）;
     U=V;
     u=abs（U）
    for j=1:33
        if a（j:）==0
           S（j）=SL（j）+conj（SL（j））*SL（j）/（conj（V（j））*V（j））*Z（j）;
        end
    end