clc;
clear;
global B1;
global B2;
global B4;
global B5;
pr=0.001;
case2_5;
n=max(B1(:,1));
nl=max(B2(:,1));
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);
P=zeros(1,n);Q=zeros(1,n);
%=====================求导纳矩阵Y==========================================
for j=1:nl
if B2(j,8)==0
p=B2(j,2);q=B2(j,3);
else q=B2(j,2);p=B2(j,3);
end
Y(p,q)=Y(p,q)-1./((B2(j,4)+B2(j,5)*i)*B2(j,7));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./((B2(j,4)+B2(j,5)*i)*B2(j,7)^2)+B2(j,6)*i;
Y(p,p)=Y(p,p)+1./(B2(j,4)+B2(j,5)*i)+B2(j,6)*i;
end
%disp('导纳矩阵Y=')
%disp(Y);
%==========================================================================
%=============生成G、B、e、f、v、P、Q=======================================
G=real(Y);B=imag(Y);%导纳矩阵的实部和虚部G B
for j=1:n
e(j)=B1(j,3)*cos(B1(j,4));%节点电压的实部e
f(j)=B1(j,3)*sin(B1(j,4));%节点电压的虚部f
V(j)=B1(j,5);%PV节点电压的幅值V
end
for j=1:n
P(j)=B1(j,6)-B1(j,8);%节点注入的有功功率P
Q(j)=B1(j,7)-B1(j,9);%节点注入的无功功率Q
B(j,j)=B(j,j)+B4(1,j);%节点无功补偿后的虚部
end
%==========================================================================
CT1=0;T2=1;N0=2*n;N=N0+1;a=0;
while T2~=0
T2=0;a=a+1;
%============================生成雅可比矩阵J==================================
for j=1:n
if j~=B(1,1)
C(j)=0;
D(j)=0;
for j1=1:n
C(j)=C(j)+G(j,j1)*e(j1)-B(j,j1)*f(j1);
D(j)=D(j)+G(j,j1)*f(j1)+B(j,j1)*e(j1);
end
P1=e(j)*C(j)+f(j)*D(j);
Q1=f(j)*C(j)-e(j)*D(j);
V2=e(j)^2+f(j)^2;
switch B1(j,2)
case 1||3
DP=P(j)-P1;
DQ=Q(j)-Q1;
for j1=1:n
if j1~=B1(1,1)&&j1~=j
H1=-B(j,j1)*e(j)+G(j,j1)*f(j);
N1=B(j,j1)*f(j)+G(j,j1)*e(j);
J1=-N1;
L1=H1;
p=2*j-1;q=2*j1-1;p2=p+1;q2=q+1;
J(p,q)=-H1;J(p,q2)=-N1;J(p,N)=DP;
J(p2,q)=-J1;J(p2,q2)=-L1;J(p2,N)=DQ;
elseif j1==j&&j1~=B(1,1)
H1=G(j,j1)*f(j)-B(j,j1)*e(j)+D(j);
N1=G(j,j1)*e(j)+B(j,j1)*f(j)+C(j);
J1=-G(j,j1)*e(j)-B(j,j1)*f(j)+C(j);
L1=G(j,j1)*f(j)-B(j,j1)*e(j)-D(j);
p=2*j-1;q=2*j1-1;p2=p+1;q2=q+1;
J(p,q)=-H1;J(p,q2)=-N1;J(p,N)=DP;
J(p2,q)=-J1;J(p2,q2)=-L1;J(p2,N)=DQ;
end
end
case 2
DP=P(j)-P1;
DV=V(j)^2-V2;
for j1=1:n
if j1~=B1(1,1)&&j1~=j
H1=-B(j,j1)*e(j)+G(j,j1)*f(j);
N1=B(j,j1)*f(j)+G(j,j1)*e(j);
R1=0;
S1=0;
p=2*j-1;q=2*j1-1;p2=p+1;q2=q+1;
J(p,q)=-H1;J(p,q2)=-N1;J(p,N)=DP;
J(p2,q)=-R1;J(p2,q2)=-S1;J(p2,N)=DV;
elseif j1==j&&j1~=B(1,1)
H1=G(j,j1)*f(j)-B(j,j1)*e(j)+D(j);
N1=G(j,j1)*e(j)+B(j,j1)*f(j)+C(j);
R1=2*f(j);
S1=2*e(j);
p=2*j-1;q=2*j1-1;p2=p+1;q2=q+1;
J(p,q)=-H1;J(p,q2)=-N1;J(p,N)=DP;
J(p2,q)=-R1;J(p2,q2)=-S1;J(p2,N)=DV;
end
end
end
end
end
%==========================================================================
%====================用高斯消去法解w=-J*V===================================
for k=3:N0
k1=k+1;N1=N;k5=k-1;
for k2=k1:N1
J(k,k2)=J(k,k2)./J(k,k);
end
J(k,k)=1;
for k3=k1:N0
for k4=k1:N1
J(k3,k4)=J(k3,k4)-J(k3,k)*J(k,k4);
end
J(k3,k)=0;
end
if k~=3
for k3=3:k5
for k4=k1:N1
J(k3,k4)=J(k3,k4)-J(k3,k)*J(k,k4);
end
J(k3,k)=0;
end
end
end
% disp(J);
%==========================================================================
for k=3:2:N0-1
L=(k+1)./2;
f(L)=f(L)-J(k,N);
k1=k+1;
e(L)=e(L)-J(k1,N);
end
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
T2=T2+1;
end
end
CT2(a)=T2+1;
CT1=CT1+1;
% for k=1;n
% DET1=abs(atan(f(k)./e(k))*180./pi);
% if DET1>1.05||DET1<0.9
% T2=0;
% end
% end
end
for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
O(k)=atan(f(k)./e(k))*180./pi;
E(k)=e(k)+f(k)*j;
end
for k=1:n
I(k)=0;
for k1=1:n
I(k)=I(k)+conj(Y(k,k1))*conj(E(k1));
end
S(k)=E(k)*I(k);
SE=real(S);SF=imag(S);
end
%=======================求首末端功率和线路损耗===============================
for j=1:nl
if B2(j,8)==0
p=B2(j,2);q=B2(j,3);
else q=B2(j,2);p=B2(j,3);
end
Si(p,q)=E(p)*(conj(E(p)*conj(B2(j,6)))+(conj(E(p)*B2(j,7))-conj(E(q)))*conj(1./((B2(j,4)+B2(j,5)*i)*B2(j,7))));
SiES(j)=real(Si(p,q));SiFS(j)=imag(Si(p,q));
end
for j=1:nl
if B2(j,8)==0
p=B2(j,2);q=B2(j,3);
else q=B2(j,2);p=B2(j,3);
end
Sj(q,p)=E(q)*(conj(E(q)*conj(B2(j,6)))+(conj(E(q)*B2(j,7))-conj(E(p)))*conj(1./((B2(j,4)+B2(j,5)*i)*B2(j,7))));
SjEM(j)=real(Sj(q,p));SjFM(j)=imag(Sj(q,p));
end
for j=1:nl
if B2(j,8)==0
p=B2(j,2);q=B2(j,3);
else q=B2(j,2);p=B2(j,3);
end
DS(j)=Si(p,q)+Sj(q,p);
DSE(j)=real(DS(j));DSF(j)=imag(DS(j));
end
%========================输出Busdate数据====================================
Busdate=zeros(n,8);
for k=1:n
Busdate(k,[1 2])=B1(k,[1 2]);
Busdate(k,[5 6 7 8])=100*B1(k,[6 7 8 9]);
Busdate(k,[3 4])=[V(k);O(k)];
end
%========================输出Branchdate数据=================================
Branchdate=zeros(nl,9);
for k=1:nl
Branchdate(k,[1,2,3])=B2(k,[1,2 3]);
Branchdate(k,4)=SiES(k);
Branchdate(k,5)=SiFS(k);
Branchdate(k,6)=SjEM(k);
Branchdate(k,7)=SjFM(k);
Branchdate(k,8)=DSE(k);
Branchdate(k,9)=DSF(k);
end
disp('精度未达到要求的个数为：');
disp(CT2);
disp('================================================================================')
disp('| Bus Data |')
disp('================================================================================')
disp(' Bus Voltage Generation Load ')
disp(' # Mag(pu) Ang(deg) P (MW) Q (MVAr) P (MW) Q (MVAr)')
disp('----- ------- -------- -------- -------- -------- --------')
disp(Busdate);
disp('')
disp('================================================================================')
disp('| Branch Data |')
disp('================================================================================')
disp('Brnch From To From Bus Injection To Bus Injection Loss (I^2 * Z) ')
disp(' # Bus Bus P (MW) Q (MVAr) P (MW) Q (MVAr) P (MW) Q (MVAr)')
disp('----- ----- ----- -------- -------- -------- -------- -------- --------')
disp(Branchdate);
