%%程序功能：IEEE57节点潮流计算 PQ分解法
%2017.7.6
%%
clear all,clc %清除工作空间
%%————————————导入阻抗数据并得到阻抗矩阵—————————————%%
A=xlsread('admittance','case57'); %从 admittance 中 case57 工作表中读入数据，格式为： 线路起点 线路终点 r x b/2 k
[m,n]=size(A); %m,n存放A矩阵的行数、列数
w=1i; %产生虚部 i
u=1; %循环系数 u
y=zeros(57,57); %产生57阶导纳0矩阵
for u=1:m; %循环求解导纳
hnode=A(u,1); %hnode为支路起点
enode=A(u,2); %enode为支路终点
z=A(u,3)+A(u,4)*w; %该线路阻抗 z
b=A(u,5)*w; %该线路对地导纳 b
k=A(u,6); %变比k，线路时为0，变压器时为k
y(hnode,enode)=y(hnode,enode)-1/(k*z); %互导纳迭代计算
y(enode,hnode)=y(hnode,enode); %互导纳y(i,j)=y(j,i)
y(hnode,hnode)=y(hnode,hnode)+1/(k*z)+(k-1)/(k*z); %自导纳计算，根据线路及变压器等效电路计算
y(enode,enode)=y(enode,enode)+1/(k*z)+(1-k)/(k*k*z); %自导纳计算，由于变压器位置而产生的修正量
y(hnode,hnode)=y(hnode,hnode)+b; %自导纳与对地导纳相加
y(enode,enode)=y(enode,enode)+b; %自导纳与对地导纳相加
end %y求解结束
N=57; %节点数 N
NPQ=50; %PQ节点数 NPQ
NPV=6; %PV节点数 NPV
D=10^-4; %精度 D
G=real(y); %电导矩阵 G
B=imag(y); %电纳矩阵 B
PQPV=xlsread('PQPV','PQPV'); %导入节点数据：节点编号 有功 无功 电压幅值，且按PQ,PV,平衡节点顺序排列
%%
%%———————迭代准备，获得B' B''及其逆矩阵，以及迭代初值dP dQ P Q————%%
for i=1:(N-1); %获得B'矩阵 B1及其逆矩阵 B1_inv
for j=1:(N-1);
B1(i,j)=imag(y(i,j));
end
end
B1_inv=inv(B1);
for i=1:NPQ %获得B''矩阵 B2及其逆矩阵B2_inv
for j=1:NPQ
B2(i,j)=imag(y(i,j));
end
end
B2_inv=inv(B2);
for i=1:(N-1); %形成P初值矩阵
Pnode(i)=PQPV(i,2);
end
for i=1:(N-1); %形成Q初值矩阵
Qnode(i)=PQPV(i,3);
end
for i=(NPQ+1):(N-1); %形成V初值矩阵
Vnode(i)=PQPV(i,4);
end
for i=1:NPQ; %形成迭代电压初值
e(i)=1;
end
for i=(NPQ+1):(N-1);
e(i)=Vnode(i);
end
e(N)=1.04;
f=zeros(1,57);
V=complex(e,f); %构成迭代电压初值矩阵复数形式V
theta=angle(V); %获得电压相角矩阵 theta
V=abs(V); %获得电压幅值矩阵 V
dP=zeros(N-1,1); %构造有功P 无功Q不平衡量矩阵 dP dQ
dQ=zeros(NPQ,1);
k=0; %迭代系数k
kp=0; %dP收敛判据 =1时，有功P收敛
kq=0; %dQ收敛判据 =1时，无功Q收敛
kmax=50; %最大迭代次数，防止死循环
%%
%%————————————迭代求解V theta——————————————————%%
for k=1:kmax;
for i=1:(N-1); %不断地累减,得dP
sum1=0;
for j=1:N;
sum1=sum1+V(i)*V(j)*(G(i,j)*cos(theta(i)-theta(j))+B(i,j)*sin(theta(i)-theta(j)));
end
dP(i)=Pnode(i)-sum1;
end
Vp=V; %构造修正相角所用的电压矩阵Vp，不含平衡节点电压
Vp(N)=[];
dtheta=(-B1_inv*(dP./Vp')); %获得相角修正量 dtheta,忽略电压幅值对相角的影响
for i=1:(N-1);
theta(i)=theta(i)+dtheta(i);
end
max1=abs(dP(1)/V(1)); %获得不平衡有功电压比最大值 dP/V
for m=1:(N-2);
if abs(dP(m)/V(m))<abs(dP(m+1)/V(m+1))
max1=abs(dP(m+1)/V(m+1));
end
end
if max1<D %如果不平衡量小于给定误差则置kp为1，表示dP已经收敛
kp=1;
end
for i=1:NPQ; %不断地累减,得dQ
sum2=0;
for j=1:N;
sum2=sum2+V(i)*V(j)*(G(i,j)*sin(theta(i)-theta(j))-B(i,j)*cos(theta(i)-theta(j)));
end
dQ(i)=Qnode(i)-sum2;
end
Vq=V; %构建修正相角所用的电压矩阵Vq，不含PV节点和平衡节点电压，
Vq(:,[51 52 53 54 55 56 57])=[]; %由于本例57节点中把50个PQ节点数据放在前面，故将51列及以后置空
dV=(-B2_inv*(dQ./Vq')); %获得电压幅值修正量 dV
max2=max(abs(dQ./Vq')); %获得不平衡无功电压比最大值
if max2<D; %如果不平衡量小于给定误差则置kq为1，表示dQ已经收敛
kq=1;
end
for i=1:NPQ; %得到修正后的电压幅值矩阵
V(i)=V(i)+dV(i);
end
format long g; %显示迭代结果并判断是否收敛，如果收敛则跳出循环，结束迭代
fprintf('\n 第%d次迭代',k);
V
theta
if (kp==1&&kq==1)
fprintf('迭代是收敛的。第%d次迭代后终止迭代。\n',k);
break;
end
end
%%
%%——————————处理节点电压，将其分解为直角坐标数据——————————%%
for i=1:57;
e(i)=V(i)*cos(theta(i));
f(i)=V(i)*sin(theta(i));
end
V=complex(e,f);
%%—————————由节点电压得节点功率，电压，线路功率———————————%%
for i=(NPQ+1):(N-1); %计算PV节点无功功率Q
for j=1:N;
Qnode(i)=Qnode(i)+f(i)*G(i,j)*e(j)-f(i)*B(i,j)*f(j)-e(i)*G(i,j)*f(j)-e(i)*B(i,j)*e(j);
end
end
Pnode(57)=0;
Qnode(57)=0;
Snode=zeros(57,57);
for j=1:N; %计算平衡节点出力
Snode(57)=Snode(57)+V(57)*(conj(y(57,j)*conj(V(j))));
end
Snode=complex(Pnode,Qnode); %将P,Q复合成S
%%计算线路功率
yi0=zeros(57,57); %生成线路对地导纳矩阵 yi0
for i=1:57;
yi0(A(i,1),A(i,2))=A(i,5); %线路对地导纳矩阵 yi0（起点，终点，b/2)
end
yi0=yi0+yi0'; %使yi0成为对称矩阵
for i=1:57; %带入计算线路功率S(i,j)
for j=(1+i):57;
S(i,j)=V(i)^2*conj(yi0(i,j))+V(i)*(conj(V(i))-conj(V(j)))*conj(-y(i,j));
S(j,i)=V(j)^2*conj(yi0(j,i))+V(j)*(conj(V(j))-conj(V(i)))*conj(-y(j,i));
end
end
V=[abs(V'),theta']; %将矩阵拼接起来V[幅值 相位]
%%————————显示最终结果 节点功率 节点电压 线路功率—————————%%
disp('节点功率：');
Snode
disp('节点电压：');
V
disp('线路功率：');
S
