A=[1 4 0 0.0576 0;
2 7 0 0.0625 0;
3 9 0 0.0586 0;
4 5 0.01 0.085 0.088;
4 6 0.017 0.092 0.079;
5 7 0.032 0.161 0.153;
6 9 0.039 0.17 0.179;
7 8 0.0085 0.072 0.0745;
8 9 0.0119 0.1008 0.1045];
y=zeros(9,9);
for k=1:9
m=A(k,1);
n=A(k,2);
y(m,m)=y(m,m)+1/(A(k,3)+1i*A(k,4))+1i*A(k,5);
y(n,n)=y(n,n)+1/(A(k,3)+1i*A(k,4))+1i*A(k,5);
y(m,n)=y(m,n)-1/(A(k,3)+1i*A(k,4));
y(n,m)=y(n,m)-1/(A(k,3)+1i*A(k,4));
end
U=[1.04,1.025,1.025,1,1,1,1,1,1];
a=[0,0,0,0,0,0,0,0,0];
P=[1,1.63,0.85,0,-1.25,-0.9,0,-1,0];
Q=[0,0,0,0,-0.5,-0.3,0,-0.35,0];
dP=[0,0,0,0,0,0,0,0,0];
dQ=[0,0,0,0,0,0,0,0,0];
G=real(y);B=imag(y);
Pt=zeros(1,9);
Qt=zeros(1,9);
H=zeros(8,8);
N=zeros(8,6);
M=zeros(6,8);
L=zeros(6,6);
Ai=zeros(1,9);
Bi=zeros(1,9);
JJ=zeros(14,14);
I=zeros(1,9);
k=0;precision=1;
while precision>0.0000001
for m=2:9
for n=1:9
pt(n)=U(m)*U(n)*(G(m,n)*cos(a(m)-a(n))+B(m,n)*sin(a(m)-a(n)));
end
dP(m)=P(m)-sum(pt);
end
for m=4:9
for n=1:9
qt(n)=U(m)*U(n)*(G(m,n)*sin(a(m)-a(n))-B(m,n)*cos(a(m)-a(n)));
end
dQ(m)=Q(m)-sum(qt);
end
for m=1:8
for n=1:8
if m==n
else
H(m,n)=-U(m+1)*U(n+1)*(G(m+1,n+1)*sin(a(m+1)-a(n+1))-B(m+1,n+1)*cos(a(m+1)-a(n+1)));
end
end
end
for m=1:8
for n=1:6
if m==n
else
N(m,n)=-U(m+1)*U(n+3)*(G(m+1,n+3)*cos(a(m+1)-a(n+3))+B(m+1,n+3)*sin(a(m+1)-a(n+3)));
end
end
end
for m=1:6
for n=1:8
if m==n
else
M(m,n)=U(m+3)*U(n+1)*(G(m+3,n+1)*cos(a(m+3)-a(n+1))+B(m+3,n+1)*sin(a(m+3)-a(n+1)));
end
end
end
for m=1:6
for n=1:6
if m==n
else
L(m,n)=-U(m+3)*U(n+3)*(G(m+3,n+3)*sin(a(m+3)-a(n+3))-B(m+3,n+3)*cos(a(m+3)-a(n+3)));
end
end
end
%雅克比i不等于j
for m=1:8
for n=1:9
Ai(n)=U(n)*(G(m+1,n)*sin(a(m+1)-a(n))-B(m+1,n)*cos(a(m+1)-a(n)));
end
H(m,m)=U(m+1)*U(m+1)*B(m+1,m+1)+U(m+1)*sum(Ai);
end
for m=1:6
for n=1:9
Bi(n)=U(n)*(G(m+1,n)*cos(a(m+1)-a(n))+B(m+1,n)*sin(a(m+1)-a(n)));
end
N(m,m)=-U(m+1)*U(m+1)*G(m+1,m+1)-U(m+1)*sum(Bi);
end
for m=1:6
for n=1:9
Bi(n)=U(n)*(G(m+3,n)*cos(a(m+3)-a(n))+B(m+3,n)*sin(a(m+3)-a(n)));
end
M(m,m)=U(m+3)*U(m+3)*G(m+3,m+3)-U(m+3)*sum(Bi);
end
for m=1:6
for n=1:9
Ai(n)=U(n)*(G(m+3,n)*sin
(a(m+3)-a(n))-B(m+3,n)*cos(a(m+3)-a(n)));
end
L(m,m)=U(m+3)*U(m+3)*B(m+3,m+3)-U(m+3)*sum(Ai);
end
%雅克比i等于j,Ai(n)=U(n)*(G(m,n)*sin(a(m)-a(n))-B(m,n)*cos(a(m)-a(n)))
%Bi(n)=U(n)*(G(m,n)*cos(a(m)-a(n))+B(m,n)*sin(a(m)-a(n)));
for m=1:8
for n=1:8
JJ(m,n)=H(m,n);
end
end
for m=1:8
for n=1:6
JJ(m,n+8)=N(m,n);
end
end
for m=1:6
for n=1:8
JJ(m+8,n)=M(m,n);
end
end
for m=1:6
for n=1:6
JJ(m+8,n+8)=L(m,n);
end
end
%雅克比矩阵JJ形成
for m=1:8
PQ(m)=dP(m+1);
end
for m=9:14
PQ(m)=dQ(m-5);
end
dUa=-inv(JJ)*PQ';
precision=max(abs(dUa));
for m=1:8
a(m+1)=a(m+1)+dUa(m);
end
for m=9:14
U(m-5)=U(m-5)+dUa(m);
end
k=k+1;
end
u=U.*cos(a)+i*U.*sin(a);
i=y*u.';
I=conj(i);
pq=u.*I.';
P=real(pq);
Q=imag(pq);
format long
XE=[1 0.3 1.137;
2 0.3 1.211;
3 0.3 1.043;
4 0 0;
5 0 0;
6 0 0;
7 0 0;
8 0 0;
9 0 0];
disp('张书槐大作业')
disp('输出修正后的发电机节点和负荷节点的自导纳')
for m=1:9
if abs(P(m))>0.000001 && abs(Q(m))>0.000001
if P(m)>0
y(m,m)=y(m,m)-1i*(1/XE(m,2));
else y(m,m)=y(m,m)-conj(pq(m))/U(m)/U(m);
end
disp(m);
disp(' 点自导纳=');
disp(y(m,m));
end
end
z=inv(y);
line=z(:,4);
If=u(4)/z(4,4);
for m=1:9
uu(m)=u(m)-line(m)*If;
end
for m=1:9
hang=A(m,:);
Y=zeros(2,2);
if hang(1)>3 && hang (2)>3
Y(1,1)=Y(1,1)+1/(A(m,3)+1i*A(m,4))+1i*A(m,5);
Y(1,2)=Y(1,2)-1/(A(m,3)+1i*A(m,4));
Y(2,1)=Y(2,1)-1/(A(m,3)+1i*A(m,4));
Y(2,2)=Y(2,2)+1/(A(m,3)+1i*A(m,4))+1i*A(m,5);
ii(m,1)=Y(1,1)*uu(hang(1))+Y(1,2)*uu(hang(2));
ii(m,2)=Y(2,1)*uu(hang(1))+Y(2,2)*uu(hang(2));
else Y(1,1)=Y(1,1)+1/(A(m,3)+1i*A(m,4));
Y(1,2)=Y(1,2)-1/(A(m,3)+1i*A(m,4));
Y(2,1)=Y(2,1)-1/(A(m,3)+1i*A(m,4));
Y(2,2)=Y(2,2)+1/(A(m,3)+1i*A(m,4));
ii(m,1)=(uu(hang(2))-uu(hang(1)))*Y(1,2);
ii(m,2)=(uu(hang(1))-uu(hang(2)))*Y(2,1);
end
end
a=zeros(9,6);
for m=1:9
a(m,1)=A(m,1);
a(m,2)=A(m,2);
a(m,3)=real(ii(m,1));
a(m,4)=imag(ii(m,1));
a(m,5)=real(ii(m,2));
a(m,6)=imag(ii(m,2));
end
b=zeros(9,6);
for m=1:9
b(m,1)=A(m,1);
b(m,2)=A(m,2);
b(m,3)=abs(ii(m,1));
b(m,4)=angle(ii(m,1))*360/(2*pi);
b(m,5)=abs(ii(m,2));
b(m,6)=angle(ii(m,2))*360/(2*pi);
end
for n=1:3
for m=1:5
s=a(m,:);
a(m,:)=a(m+1,:);
a(m+1,:)=s;
m=m+1;
end
n=n+1;
end
for n=1:3
for m=1:5
s=a(m+3,:);
a(m+3,:)=a(m+4,:);
a(m+4,:)=s;
m=m+1;
end
n=n+1;
end
for n=1:3
for m=1:5
s=b(m,:);
b(m,:)=b(m+1,:);
b(m+1,:)=s;
m=m+1;
end
n=n+1;
end
for n=1:3
for m=1:5
s=b(m+3,:);
b(m+3,:)=b(m+4,:);
b(m+4,:)=s;
m=m+1;
end
n=n+1;
end
disp('输出节点阻抗矩阵的第4列')
line
disp('短路电流If的模值和角度')
disp(' 模值')
abs(If)
disp(' 相角')
angle(If)*360/(2*pi)
disp('短路时各节点电压模值')
sb=abs(uu);
sb.'
disp('短路时各支路电流')
disp(' i j Iij实部 Iij虚部 Iji实部 Iji虚部')
a
disp(' i j Iij模值 Iij相角（角度制） Iji模值 Iji相角（角度制）')
b下载本文
