MATLAB,MATPOWER 生成导纳矩阵?

clc
G(1,1)=3.750;B(1,1)=-11.250;G(1,2)=-2.500;B(1,2)=7.500;
G(1,3)=0;B(1,3)=0;G(1,4)=-1.250;B(1,4)=3.750;
G(1,5)=0;B(1,5)=0;
G(2,1)=-2.500;B(2,1)=7.500;G(2,2)=10.834;B(2,2)=-32.500;G(2,3)=-1.667;B(2,3)=5.000;
G(2,4)=-1.667;B(2,4)=5.000;G(2,5)=-5.000;B(2,5)=15.000;
G(3,1)=0;B(3,1)=0;G(3,2)=-1.667;B(3,2)=5.000;G(3,3)=12.917;B(3,3)=-38.75;
G(3,4)=-10.000;B(3,4)=30.000;G(3,5)=-1.250;B(3,5)=3.750;
G(4,1)=-1.250;B(4,1)=3.750;G(4,2)=-1.667;B(4,2)=5.000;G(4,3)=-10.000;B(4,3)=30.000;
G(4,4)=12.917;B(4,4)=-38.750;G(4,5)=0;B(4,5)=0;
G(5,1)=0;B(5,1)=0;G(5,2)=-5.000;B(5,2)=15.000;G(5,3)=-1.250;B(5,3)=3.750;
G(5,4)=0;B(5,4)=0;G(5,5)=6.25;B(5,5)=-18.750;
Y=G+j*B
delt(1)=0; delt(2)=0; delt(3)=0; delt(4)=0;
u(1)=1.0; u(2)=1.0; u(3)=1.0; u(4)=1.0;
ps(1)=-0.60;qs(1)=-0.10;ps(2)=0.20;qs(2)=0.20;
ps(3)=-0.45;qs(3)=-0.15; ps(4)=-0.40;qs(4)=-0.05;
k=1;precision=1;
N1=4;
while precision>0.00001
delt(5)=0;u(5)=1.06;
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


pi(m)=sum(pt);qi(m)=sum(qt);


dp(m)=ps(m)-pi(m);


dq(m)=qs(m)-qi(m);
end
for m=1:N1


for n=1:N1




if m==n




H(m,m)=-qi(m)-u(m)^2*B(m,m); N(m,m)=pi(m)+u(m)^2*G(m,m);




J(m,m)=pi(m)-u(m)^2*G(m,m); L(m,m)=qi(m)-u(m)^2*B(m,m);




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





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




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
JJ


for m=1:N1




PP(2*m-1)=dp(m);




PP(2*m)=dq(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)*u(n);
end
k=k+1;
end
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(5,m)*U(m);
end
S5=U(5)*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
end
S



%
用过的例子,若有标点问题,请调试后再用.