前推回代法计算IEEE33节点潮流

clear
clc
load Bus_33.txt
Bus=Bus_33;
% Bus（:2）=Bus（:2）*0.71;
% Bus（:1）=Bus（:1）*0.71;
Bus（182）=Bus（182）-8000;
% Bus（183）=Bus（183）+800;
load Branch_33.txt
Branch=Branch_33;
[busnumcolumn]=size（Bus）;
[branchnumcolumn]=size（Branch）;
shoubus=Branch（:2）;            %首节点
mobus=Branch（:3）;              %末节点
Vbus=ones（busnum1）;            %节点的初始电压都为1
Vbus（:1）=12.66;                %电压为12.66kV
Vbus1=Vbus;
Ploss=zeros（busnum1）;
Qloss=zeros（busnum1）;
e=1;
k=0;
Branch1=Branch;
n=1;
while ~isempty（Branch1）
m=1;
[srow]=size（Branch1）;
while s>0
t=find（Branch1（:2）== Branch1（s3））;
if isempty（t）
T1（n:）= Branch1（s:）;
n=n+1;
else
T2（m:）= Branch1（s:）;
m=m+1;
end;
s=s-1;
end;
Branch1=T2;
T2=[];
end;
while e>1.0e-05
        P=zeros（busnum1）;
        Q=zeros（busnum1）;
for s=1:branchnum