求MATLAB简单潮流计算程序。。。。

function lianxuchaoliu
clear;
clc;
n=9;%节点数；
nl=9;%支路数；
isb=1;%平衡节点号；
pr=0.00001;%误差精度；
b1=[1 4 0.0576i 0 1.05 1; 4 5 0.017+0.092i 0.158i 1 0; 5 6 0.039+0.17i 0.358i 1 0; 3 6 0.0586i 0 1.05 1; 6 7 0.0119+0.1008i 0.209i 1 0; 7 8 0.0085+0.072i 0.149i 1 0; 8 2 0.0625i 0 1.05 1; 8 9 0.032+0.161i 0.306i 1 0; 9 4 0.01+0.085i 0.176i 1 0];
%依次是支路首端；末端，支路阻抗；对地电纳；支路变比；折算到哪一侧标志（高压侧为1；低压侧为0）；
b2=[0 0 1.05 1.05 0 1; 1.63 0 1.05 1.05 0 3; 0.85 0 1.05 1.05 0 3; 0 0 1 0 0 2; 0 0.9+0.3i 1 0 0 2; 0 0 1 0 0 2; 0 1+0.35i 1 0 0 2; 0 0 1 0 0 2; 0 1.25+0.5i 1 0 0 2];
%依次是节点的发电机功率；负荷功率；节点电压初值；PV节点电压V给定值；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；
Y=zeros(n);%求导纳阵；
for i=1:nl
if b1(i,6)==0
p=b1(i,1);q=b1(i,2);
else
p=b1(i,2);q=b1(i,1);
end
Y(p,q)=Y(p,q)-1./(b1(i,3)*b1(i,5));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./(b1(i,3)*b1(i,5)^2)+b1(i,4)./2;
Y(p,p)=Y(p,p)+1./b1(i,3)+b1(i,4)./2;
end
%disp('系统的导纳阵为：');
%disp(Y);
G=real(Y);B=imag(Y);
for i=1:n
e(i)=real(b2(i,3));
f(i)=imag(b2(i,3));
v(i)=b2(i,4);
end
for i=1:n
S(i)=b2(i,1)-b2(i,2);
B(i,i)=B(i,i)+b2(i,5);
end
P=real(S);Q=imag(S);
w=zeros(2*n,1);Jac=zeros(2*n);
t=0;
while t==0
for i=1:n
if b2(i,6)~=isb
C=0;D=0;
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j);
D=D+G(i,j)*f(j)+B(i,j)*e(j);
end
if b2(i,6)==2%P，Q节点；
w(2*i)=P(i)-e(i)*C-f(i)*D;
w(2*i-1)=Q(i)-f(i)*C+e(i)*D;
else if b2(i,6)==3%P，V节点；
w(2*i)=P(i)-e(i)*C-f(i)*D;
w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2);
end
end
else
w(2*i-1)=0;
w(2*i)=0;
end
end
%disp(w);
w1=w(3:2*n);
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P，Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
end
end
else if b2(i,6)==3%P，V节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i-1,2*j)=0;
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
Jac(2*i-1,2*j-1)=0;
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
Jac(2*i-1,2*j)=-2*f(i);
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
Jac(2*i-1,2*j-1)=-2*e(i);
end
end
end
end
else
Jac(2*i-1,2*j-1)=0;
Jac(2*i,2*j)=0;
Jac(2*i-1,2*j)=0;
Jac(2*i,2*j-1)=0;
end
end
end
%disp(Jac);
Jac1=Jac(3:2*n,3:2*n);
for k=1:2*n-2
m=0;
for i=k+1:2*n-2
m=Jac1(i,k)./Jac1(k,k);
w1(i)=w1(i)-m*w1(k);
for j=k+1:2*n-2
Jac1(i,j)=Jac1(i,j)-m*Jac1(k,j);
end
end
end
E(2*n-2)=w1(2*n-2)./Jac1(2*n-2,2*n-2);
for i=2*n-3:-1:1
c=0;
for k=i+1:2*n-2
c=c+Jac1(i,k)*E(k);
E(i)=(w1(i)-c)./Jac1(i,i);
end
end
%disp(E);
for i=1:n-1
e(i+1)=e(i+1)-E(2*i-1);
f(i+1)=f(i+1)-E(2*i);
end
%disp(e);
%disp(f);
for i=1:n-1
b(2*i-1)=abs(E(2*i-1));
b(2*i)=abs(E(2*i));
end
KB=max(b);
%disp(KB);
if KBt=1;
else
t=0;
end
end
%disp(e);
%disp(f);
for i=1:n
fz(i)=sqrt(e(i)^2+f(i)^2);
end
disp(fz);
JJ1=zeros(2*n-1,2*n-1);Jac1=zeros(2*n-2,2*n-2);ttt=0;ccc=1;%迭代次数；
L=0.05;T=0;%L是步长，T是参数；
while ttt==0
%确定扩展矩阵
for i=1:n-1
if b2(i+1,6)==3
JJ1(2*i,2*n-1)=1.5*1.05;
JJ1(2*i-1,2*n-1)=0;
else if b2(i+1,6)==2
if b2(i+1,2)==0
JJ1(2*i-1,2*n-1)=0;
JJ1(2*i,2*n-1)=0;
else
JJ1(2*i,2*n-1)=-real(b2(i+1,2));
JJ1(2*i-1,2*n-1)=-imag(b2(i+1,2));
end
end
end
end
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2% P,Q节点；
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
end
end
else if b2(i,6)==3%P,V节点，
if j~=i
Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
Jac(2*i-1,2*j)=0;
Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
Jac(2*i-1,2*j-1)=0;
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
Jac(2*i-1,2*j)=-2*f(i);
Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
Jac(2*i-1,2*j-1)=-2*e(i);
end
end
end
end
else
Jac(2*i-1,2*j-1)=0;
Jac(2*i,2*j)=0;
Jac(2*i-1,2*j)=0;
Jac(2*i,2*j-1)=0;
end
end
end
Jac1=Jac(3:2*n,3:2*n);
for i=1:2*n-2
for j=1:2*n-2
JJ1(i,j)=Jac1(i,j);
end
end
if ccc==1%以负荷为连续参数；
for j=1:2*n-2
JJ1(2*n-1,j)=0;
end
JJ1(2*n-1,2*n-1)=1;
w2=zeros(2*n-1,1);
for i=1:2*n-2
w2(i,1)=0;
end
w2(2*n-1,1)=1;
end
%disp(JJ1);
if ccc>1%以切向量中分量最大（绝对值最大）的状态变量选定为连续参数；
for i=1:2*n-2
Jd(i)=abs(d(i));
end
for i=1:2*n-3
if Jd(i)>=Jd(i+1)
zd=Jd(i);ek=i;
else if Jd(i)<=Jd(i+1)
zd=Jd(i+1);ek=i+1;
end
end
end
for j=1:2*n-1
JJ1(2*n-1,j)=0;
end
JJ1(2*n-1,ek)=1;
w2=zeros(2*n-1,1);
for i=1:2*n-1
w2(i,1)=0;
end
if d(ek)>0
w2(ek,1)=1;
else if d(ek)<0
w2(ek,1)=-1;
end
end
end
for k=1:2*n-1
m=0;
for i=k+1:2*n-1
m=JJ1(i,k)./JJ1(k,k);
w2(i)=w2(i)-m*w2(k);
for j=k+1:2*n-1
JJ1(i,j)=JJ1(i,j)-m*JJ1(k,j);
end
end
end
d(2*n-1)=w2(2*n-1)./JJ1(2*n-1,2*n-1);
for i=2*n-2:-1:1
c=0;
for k=i+1:2*n-1
c=c+JJ1(i,k)*d(k);
d(i)=(w2(i)-c)./JJ1(i,i);
end
end
%disp(d);
for i=1:n-1
e(i+1)=e(i+1)+L*d(2*i-1);
f(i+1)=f(i+1)+L*d(2*i);
end
T=T+L*d(2*n-1);
disp(d(2*n-1));
%disp(e);
%disp(f);
%disp(T);
%对预估的近似解进行校正；
tt=0;
while tt==0
for i=1:n
if b2(i,6)~=isb
C=0;D=0;
for j=1:n
C=C+G(i,j)*e(j)-B(i,j)*f(j);
D=D+G(i,j)*f(j)+B(i,j)*e(j);
end
if b2(i,6)==2% P,Q节点；
if b2(i,2)~=0
wj(2*i)=P(i)-T*real(b2(i,2))-e(i)*C-f(i)*D;
wj(2*i-1)=Q(i)-T*imag(b2(i,2))-f(i)*C+e(i)*D;
else
wj(2*i)=P(i)-e(i)*C-f(i)*D;
wj(2*i-1)=Q(i)-f(i)*C+e(i)*D;
end
else if b2(i,6)==3%P,V节点；
wj(2*i)=P(i)+T*1.5*1.05-e(i)*C-f(i)*D;
wj(2*i-1)=0;
end
end
else
wj(2*i)=0;
wj(2*i-1)=0;
end
end
wj1=wj(3:2*n);
for i=1:n
for j=1:n
if b2(i,6)~=isb
if b2(i,6)==2%P,Q节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
JJac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
JJac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
JJac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
JJac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
end
end
else if b2(i,6)==3%P,V节点；
if j~=i
JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
JJac(2*i-1,2*j)=0;
JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
JJac(2*i-1,2*j-1)=0;
else if j==i
m=0;h=0;
for r=1:n
m=m+G(i,r)*e(r)-B(i,r)*f(r);
h=h+G(i,r)*f(r)+B(i,r)*e(r);
end
JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
JJac(2*i-1,2*j)=-2*f(i);
JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
JJac(2*i-1,2*j-1)=-2*e(i);
end
end
end
end
else
JJac(2*i-1,2*j-1)=0;
JJac(2*i,2*j)=0;
JJac(2*i-1,2*j)=0;
JJac(2*i,2*j-1)=0;
end
end
end
JJac1=JJac(3:2*n,3:2*n);
for k=1:2*n-2
m=0;
for i=k+1:2*n-2
m=JJac1(i,k)./JJac1(k,k);
wj1(i)=wj1(i)-m*wj1(k);
for j=k+1:2*n-2
JJac1(i,j)=JJac1(i,j)-m*JJac1(k,j);
end
end
end
E1(2*n-2)=wj1(2*n-2)./JJac1(2*n-2,2*n-2);
for i=2*n-3:-1:1
c=0;
for k=i+1:2*n-2
c=c+JJac1(i,k)*E1(k);
E1(i)=(wj1(i)-c)./JJac1(i,i);
end
end
%disp(E1);
for i=1:n-1
e(i+1)=e(i+1)-E1(2*i-1);
f(i+1)=f(i+1)-E1(2*i);
end
for i=1:n-1
bx(2*i-1)=abs(E1(2*i-1));
bx(2*i)=abs(E1(2*i));
end
KB1=max(bx);
if KB1tt=1;
else
tt=0;
end
end
%disp(e);
%disp(f);
if d(2*n-1)>0
ttt=0;ccc=ccc+1;
else
ttt=1;
end
end
disp(T);
disp(ccc);
%disp(e);
%disp(f);
for i=1:n
fz1(i)=sqrt(e(i)^2+f(i)^2);
end
disp(fz1);