matlab程序 潮流计算,潮流计算MATLAB程序运行报错，求大神指导！！！

>> clear

G(1,1)=0.1595;B(1,1)=-0.356;

G(1,2)=-0.076;B(1,2)=0.18;G(1,3)=0;B(1,3)=0;

G(1,4)=-0.0345;B(1,4)=0.083;G(1,5)=-0.049;B(1,5)=0.093;

G(2,1)=-0.076;B(2,1)=0.18;

G(2,2)=0.167;B(2,2)=-0.269;G(2,3)=-0.091;B(2,3)=0.089;

G(2,4)=0;B(2,4)=0;G(2,5)=0;B(2,5)=0;

B(3,1)=0;G(3,2)=-0.091;B(3,2)=0.089;G(3,3)=0.091;B(3,3)=-0.089;

G(3,4)=0;B(3,4)=0;G(3,5)=0;B(3,5)=0;G(3,1)=0;

B(4,1)=0.083;G(4,2)=0;B(4,2)=0;G(4,3)=0;B(4,3)=0;

G(4,4)=0.0738;B(4,4)=-0.1572;G(4,5)=-0.0393;B(4,5)=0.0742;G(4,1)=-0.0345;

G(5,1)=-0.049;B(5,1)=0.093;G(5,2)=0;B(5,2)=0;G(5,3)=0;B(5,3)=0;

G(5,4)=-0.0393;B(5,4)=0.0742;G(5,5)=0.0883;B(5,5)=-0.1672;

Y=G+j*B;

delt(2)=0;delt(3)=0;delt(4)=0;delt(5)=0;

u(2)=1.0;u(3)=1.0;u(4)=1.0;u(5)=1.0;

p(2)=36.184;q(2)=20.44;p(3)=20.114;q(3)=11.51;

p(4)=25.138;q(4)=14.378;p(5)=22.13;q(5)=12.82;

k=0;precision=1;

N1=4;

while precision>0.00001

delt(1)=0;u(1)=110;

for m=1:N1

for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

end

pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);

end

for m=1:N1

for n=1:N1+1

h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

L0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

end

H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));

N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));

J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));

L(m,m)=sum(L0)+2*u(m)^2*B(m,m)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));

end

for m=1:N1

for n=1:N1

if m==n

else

H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));

J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));

N(m,n)=-J(m,n);L(m,n)=H(m,n);

JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);

JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);

end

end

end

for m=1:N1

PP(2*m-1)=pp(m);PP(2*m)=qq(m);

end

uu=-inv(JJ)*PP;precision=max(abs(uu));

for n=1:N1

delt(n)=delt(n)+uu(2*n-1);

u(n)=u(n)+uu(2*n);

end

k=k+1;

end

k-1,delt',u'

for n=1:N1+1

U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));

end

for m=1:N1+1

I(m)=Y(1,m)*U(m);

end

S1=U(1)*sum(conj(I))

for m=1:N1+1

for n=1:N1+1

S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n));

end

end

S

QQ截图20150506185836.png