IEEE33节点配电网，模型和程序都有

clear all; Bus=[1,0,0 ; 2,100, 60; 3,90,40; 4,120,80; 5,60,30; 6,60,20; 7, 200, 100 ; 8,200, 100 ; 9, 60, 20; 10,60, 20; 11,45,30; 12,60, 35; 13,60, 35; 14,120,80; 15,60, 10; 16,60, 20; 17,60, 20; 18,90, 40; 19,90, 40; 20,90, 40; 21,90, 40; 22,90, 40; 23,90, 40; 24,420, 200; 25,420, 200 ; 26,60, 25; 27,60 ,25; 28,60, 20; 29,120, 70; 30,200,600 ; 31,150, 70; 32,210, 100 ; 33,60,40; ]; Branch=[1 ,1 ,2 ,0.0922,0.0407; 2 ,2 ,3 ,0.4930,0.2511; 3 ,3 ,4 ,0.3660,0.1864; 4 ,4 ,5 ,0.3811,0.1941; 5 ,5 ,6 ,0.8190,0.7070; 6 ,6 ,7 ,0.1872,0.6188; 7 ,7 ,8 ,0.7144,0.2351; 8 ,8 ,9 ,1.0300,0.7400; 9 ,9 ,10,1.0440,0.7400; 10,10,11,0.1966,0.065; 11,11,12,0.3744,0.1238; 12,12,13,1.4680,1.1550; 13,13,14,0.5416,0.7129; 14,14,15,0.5910,0.5260; 15,15,16,0.7463,0.5450; 16,16,17,1.2890,1.7210; 17,17,18,0.7320,0.5740; 18,2, 19,0.1640,0.1565; 19,19,20,1.5042,1.3554; 20,20,21,0.4095,0.4784; 21,21,22,0.7089,0.9373; 22,3, 23,0.4512,0.3083; 23,23,24,0.8980,0.7091; 24,24,25,0.8960,0.7011; 25,6, 26,0.2030,0.1034; 26,26,27,0.2842,0.1447; 27,27,28,1.0590,0.9337; 28,28,29,0.8042,0.7006; 29,29,30,0.5075,0.2585; 30,30,31,0.9744,0.9630; 31,31,32,0.3105,0.3619; 32,32,33,0.3410,0.5302; ]; [busnum,row]=size(Bus); [branchnum,row]=size(Branch); soubus=Branch(:,2); mobus=Branch(:,3); Vbus=ones(busnum,1); Vbus(:,1)=12.66; Vbus1=Vbus; Ploss=zeros(busnum,1); Qloss=zeros(busnum,1); e=1; k=0; Branch1=Branch; n=1; while ~isempty(Branch1) m=1; [s,row]=size(Branch1); while s>0 t=find(Branch1(:,2)== Branch1(s,3)); 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(busnum,1); Q=zeros(busnum,1); for s=1:branchnum i=T1(s,2); j=T1(s,3); R=T1(s,4); X=T1(s,5); Pload=Bus(j,2); Qload=Bus(j,3); II=((Pload+P(j))^2+(Qload+Q(j))^2)/(Vbus(j)^2*1000); Ploss(i,j)=II*R; Qloss(i,j)=II*X; P(i,j)=Pload+Ploss(i,j)+P(j); Q(i,j)=Qload+Ploss(i,j)+Q(j); P(i)=P(i)+P(i,j); Q(i)=Q(i)+Q(i,j); end; for s=branchnum:-1:1 i=T1(s,2); j=T1(s,3); R=T1(s,4); X=T1(s,5); Vbus(j)=(Vbus(i)-(P(i,j)*R+Q(i,j)*X)/(Vbus(i)*1000))^2+((P(i,j)*X-Q(i,j)*R)/(Vbus(i)*1000))^2; Vbus(j)=sqrt(Vbus(j)); end; e=max(abs(Vbus1-Vbus)); Vbus1=Vbus; k=k+1; end