%% 输入参数
n=input('请输入节点数:n=');
n1=input('请输入支路数:n1=');
l=input('请输入精度:l=');
A1=input('请输入支路矩阵:A1=');
A2=input('请输入节点矩阵:A2=');
%% 声明变量
Y=zeros(n); %初始节点导纳矩阵
G0=zeros(n);
B0=zeros(n);
JD=zeros(n,5); %重排节点矩阵，行序：先PQ节点，再PV节点，最后平衡节点；列序：节点，有功功率，无功功率，电压幅值，电压相角
B=zeros(n); %重排电纳和电导矩阵
G=zeros(n);
B1=zeros(n-1); %B'矩阵
%% 初始电导和电纳矩阵
Z=A1(:,3)+A1(:,4)*1i; %设置中间变量导纳
Zy=A1(:,5)+A1(:,6)*1i;
for i=1:n1
p=A1(i,1);
q=A1(i,2);
Y(p,q)=Y(p,q)-1./(Z(i,1)*A1(i,7)); %求互导纳,考虑变压器变比
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1./Z(i,1)+Zy(i,1); %求自导纳
Y(q,q)=Y(q,q)+1./(Z(i,1)*A1(i,7)^2)+Zy(i,1);
end
G0=real(Y); %导纳实部为G0
B0=imag(Y); %导纳虚部为B0
%% 重排节点，给P,Q,V,A赋初值，求B'和B''矩阵
j1=0;
for i=1:n
if A2(i,2)==1 %PQ节点，已知有功和无功，电压幅值和相角为1和0
j1=j1+1;
JD(j1,1)=A2(i,1);
JD(j1,2)=A2(i,3);
JD(j1,3)=A2(i,4);
JD(j1,4)=1;
JD(j1,5)=0;
end
end
m=j1; %PQ节点的个数
j2=m;
for i=1:n
if A2(i,2)==2 %PV节点，已知有功和电压幅值，无功和相角为0
j2=j2+1;
JD(j2,1)=A2(i,1);
JD(j2,2)=A2(i,3);
JD(j2,3)=0;
JD(j2,4)=A2(i,5);
JD(j2,5)=0;
else
JD(n,1)=A2(i,1); %平衡节点，已知电压幅值和相角
JD(n,2)=0;
JD(n,3)=0;
JD(n,4)=A2(i,5);
JD(n,5)=A2(i,6);
end
end
for i=1:n %求B矩阵,G矩阵
for j=1:n
B(i,j)=B0(JD(i,1),JD(j,1));
G(i,j)=G0(JD(i,1),JD(j,1));
end
end
B1=B(1:n-1,1:n-1); %B' 矩阵
B2=B(1:m,1:m); %B''矩阵
%disp(B1);
%disp(B2);
%% 迭代计算
delt_P=zeros(n-1,1); %初始化P,Q,V,A为零矩阵
delt_A=zeros(n-1,1);
delt_Q=zeros(m,1);
delt_V=zeros(m,1);
P=zeros(n-1,1); %归一化矩阵
Q=zeros(m,1); %归一化矩阵
X=zeros(n-1,1);
k=0; %迭代次数
Kp=1; %有功迭代收敛标志
Kq=1; %无功迭代收敛标志
%开始迭代
while k<=30 %限制迭代次数最大30
for i=1:n-1
a=0;
for j=1:n
a=a+JD(j,4).*(G(i,j).*cosd(JD(i,5)-JD(j,5))+B(i,j).*sind(JD(i,5)-JD(j,5)));
end
delt_P(i,1)=JD(i,2)-JD(i,4)*a;
P(i,1)=delt_P(i,1)./JD(i,4);
end
if max(abs(delt_P))>=l
X=-(B1\P);
delt_A=(X./JD(1:n-1,4))*360/2/pi; %计算相角差，并转化为角度
JD(1:n-1,5)=JD(1:n-1,5)+delt_A(1:n-1,1);
Kq=1;
else
Kp=0;
end
if Kq~=0
for i=1:m
b=0;
for j=1:n
b=b+JD(j,4).*(G(i,j).*sind(JD(i,5)-JD(j,5))-B(i,j).*cosd(JD(i,5)-JD(j,5)));
end
delt_Q(i,1)=JD(i,3)-JD(i,4)*b;
Q(i,1)=delt_Q(i,1)./JD(i,4);
end
if max(abs(delt_Q))>=l
delt_V=-(B2\Q);
JD(1:m,4)=JD(1:m,4)+delt_V(1:m,1);
Kp=1;
else
Kq=0;
end
if Kp~=0
k=k+1;
else
break
end
else
break
end
end
disp_NVA=zeros(n,3);
for i=1:n
disp_NVA(JD(i,1),1)=JD(i,1);
disp_NVA(JD(i,1),2)=JD(i,4);
disp_NVA(JD(i,1),3)=JD(i,5);
end
disp('共迭代的次数为：');
disp(k);
disp('修正结果：节点 电压幅值 电压相角');
disp(disp_NVA);
%--------------
%计算平衡节点功率和支路功率损耗
S=zeros(n);
V2=zeros(n,1);
delt_P=zeros(n1,1);
delt_Pall=0;
V2=disp_NVA(:,2).*(cosd(disp_NVA(:,3))+sind(disp_NVA(:,3))*1i); %dis_NVA的复数表示形式
for k=1:n1
i=A1(k,1);
j=A1(k,2);
S(i,j)=(disp_NVA(i,2)^2)*(-A1(k,6)*1i)+V2(i,1)*(conj(V2(i,1))-conj(V2(j,1))*(1./A1(k,7)))*conj(1./(A1(k,3)+A1(k,4)*1i));
S(j,i)=(disp_NVA(j,2)^2)*(-A1(k,6)*1i)+(1./A1(k,7))*V2(j,1)*((1./A1(k,7))*conj(V2(j,1))-conj(V2(i,1)))*conj(1./(A1(k,3)+A1(k,4)*1i));
delt_P(k)=abs(real(S(i,j)+S(j,i)));
delt_Pall=delt_Pall+delt_P(k);
end
disp('输电线路的视在功率矩阵为：');
disp(S);
%-------------
%计算平衡节点的视在功率
S_N=0;
c=0;
for j=1:n
c=c+conj(Y(n,j))*conj(V2(j,1));
end
S_N=V2(n,1)*c;
xiaolv=(real(S_N)-delt_Pall)/real(S_N);
disp('平衡节点的视在功率为：');
disp(S_N);
disp('系统的网损为：');
disp(delt_Pall);
disp('系统的传输效率为：');
disp(xiaolv);
