clear
n=4;
n1=1;
n2=2;
isb=4;
pr=0.00001;
K=[0 0 0 0;0 0 0 0.9625;0 0 0 0;0 0 0 0];
C=[0 0.02+0.06i 0.01+0.03i inf; 0.02+0.06i 0 0.03+0.07i 0.0+0.05i;0.01+0.03i 0.03+0.07i 0 0.02+0.05i;inf 0.0+0.05i 0.02+0.05i 0];
y=[0 0.01i 0.01i 0;0.01i 0 0 0;0.01i 0 0 0;0 0 0 0];
U=[1+0i 1+0i 1.02+0i 1.05+0i];
S=[-0.4-0.3i -0.3-0.2i 0.4 0];
Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros(1,n1+n2); for m=1:n
for R=1:n
C(m,m)=C(m,m)+y(m,R);
if K(m,R)~=0
C(m,m)=C(m,m)+1/(C(m,R) /( K(m,R) * (K(m,R)-1))) ;
C(R,R)=C(R,R)+1/(C(m,R)/(1-K(m,R)));
C(m,R)=C(m,R)/K(m,R);
C(R,m)=C(m,R);
end
end
end
for m=1:n
for R=1:n
if m~=R
Z(m)=Z(m)+1/C(m,R);
end
end
end
for m=1:n
for R=1:n
if m==R
Y(m,m)=C(m,m)+Z(m);
else
Y(m,R)=-1/C(m,R);
end
end
end
disp('结点导纳矩阵:');
disp(Y);
disp('迭代中的雅克比矩阵:');
G=real(Y);
B=imag(Y);
O=angle(U);
U1=abs(U);
k=0;
PR=1;
P=real(S);
Q=imag(S);
while PR>pr
for m=1:n2
UD(m)=U1(m);
end
for m=1:n1+n2
for R=1:n
PT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R )));
end
PT1(m)=sum(PT);
PP(m)=P(m)-PT1(m);
PP1(k+1,m)=PP(m);
end
for m=1:n2
for R=1:n
QT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R )));
end
QT1(m)=sum(QT);
QQ(m)=Q(m)-QT1(m);
QQ1(k+1,m)=QQ(m);
end
PR1=max(abs(PP));
PR2=max(abs(QQ));
PR=max(PR1,PR2);
for m=1:n1+n2
for R=1:n1+n2
if m==R
H(m,m)=U1(m)^2*B(m,m)+QT1(m);
else
H(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O( R)));
end
end
end
for m=1:n1+n2
for R=1:n2
if m==R
N(m,m)=-U1(m)^2*G(m,m)-PT1(m);
else
N(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O (R)));
end
end
end
for m=1:n2
for R=1:n1+n2
if m==R
J(m,m)=U1(m)^2*G(m,m)-PT1(m);
else
J(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O( R)));
end
end
end
for m=1:n2
for R=1:n2
if m==R
L(m,m)=U1(m)^2*B(m,m)-QT1(m);
else
L(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O( R)));
end
end
end
JJ=[H N;J L];
disp(JJ);
PQ=[PP';QQ'];
DA=-inv(JJ)*PQ;
DA1=DA';
for m=1:n1+n2
OO(m)=DA1(m);
end
for m=n:n1+n2+n2
UU1(m-n1-n2)=DA1(m);
end
UD2=diag(UD);
UU=UU1*UD2;
for m=1:n1+n2
O(m)=O(m)+OO(m);
end
for m=1:n2
U1(m)=U1(m)+UU(m);
end
for m=1:n1+n2
o(k+1,m)=180/pi*O(m);
end
for m=1:n2
u(k+1,m)=U1(m);
end
k=k+1;
end
for m=1:n
b(m)=U1(m)*cos(O(m));
c(m)=U1(m)*sin(O(m));
end
U=b+i*c;
for R=1:n
PH1(R)=U(isb)*conj(Y(isb,R))*conj(U(R)); end
PH=sum(PH1);
for m=1:n
for R=1:n
if m~=R
C1(m,R)=1/C(m,R);
else
C1(m,m)=C(m,m);
end
end
end
for m=1:n
for R=1:n
SS(m,R)=U1(m)^2*conj(C1(m,m))+U(m)*(conj(U(m))-conj(U(R)))*conj (C1(m,R));
end
end
SS
disp('迭代中的△P：');disp(PP1);
disp('迭代中的△Q：');disp(QQ1);
disp('迭代中相角:');disp(o);
disp('迭代中电压的模:');disp(u);
disp('平衡结点的功率:');disp(PH);
disp('全部线路功率分布:');disp(SS);。