潮流计算pq分解法概述
clear;clc
n=5;%节点个数
nl=5;%支路条数
max1=100;%最大迭代次数
pr=0.00001;%精度
Z=[1 2 0.04+0.25j 0.25j 0.25j 0;
1 3 0.1+0.35j 0 0 0;
2 3 0.08+0.3j 0.25j 0.25j 0;
2 4 1.05 0.015j 0 1
3 5 1.05 0.03j 0 1];%由支路参数形成的矩阵Z=[节点号1 节点号2 支路阻抗靠近节点1的对地导纳靠近节点2的对地导纳参数]'); 参数的详细资料请查阅本程序附带的文本。Z1=zeros(n);
A=[-1.6-0.8j 1 1;
-2-1j 1 1;
-3.7-1.3j 1 1;
5 1.05 2;
0 1.05 3];%节点的功率，电压矩阵 A=[节点输入功率节点电压节点类型（PQ节点的节点类型是1；PV：2；平衡节点：3）]
%PQ节点的电压初值设为1；平衡节点的功率S初值设为0
YZ=[0.25,0.515,0.28,0,0];
for i=1:n
p=Z(i,1);q=Z(i,2);
Z1(p,q)=Z(i,3);
Z1(q,p)=Z1(p,q);
end
for i=1:nl
if Z(i,6)==1
K=Z(i,3);Zt=Z(i,4);
Z(i,3)=K*Zt;Z(i,4)=(1-K)/(K^2*Zt);Z(i,5)=(K-1)/(K*Zt);
elseif Z(i,6)==2
K=Z(i,3);Zt=Z(i,4);
Z(i,3)=Zt/K;Z(i,4)=K*(K-1)/Zt;Z(i,5)=(1-K)/Zt;
elseif Z(i,6)==3
Zt=Z(i,3);K=Z(i,4);
Z(i,3)=K*Zt;Z(i,4)=(K-1)/(K*Zt);Z(i,5)=(1-K)/(K^2*Zt);
elseif Z(i,6)==4
Zt=Z(i,3);K=Z(i,4);
Z(i,3)=Zt/K;Z(i,4)=(1-K)/Zt;Z(i,5)=K*(K-1)/Zt;
else
end
end%对含有变压器的支路进行处理，详见本程序附带的文本。
Y=zeros(n);
% 计算节点导纳矩阵
for i=1:nl
p=Z(i,1);q=Z(i,2);%p是支路始端编号，q是支路末端编号
Y(p,q)=-1/Z(i,3);
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1/Z(i,3)+Z(i,4);
Y(q,q)=Y(q,q)+1/Z(i,3)+Z(i,5);
end
disp('结点导纳矩阵Y为：'); %输出节点导纳Y矩阵
disp(Y);
Z2=1./(Y+eps);
for i=1:n
for j=1:n
if Z2(i,j)==Inf
Z2(i,j)=0;
end
end
end
G=real(Y);%节点电导G矩阵
B=imag(Y);%节点电纳B矩阵
B11=B(1:4,1:4);
B1=zeros(n);
B2=zeros(n);
R1=real(Z2);
X1=imag(Z2);
for i=1:n
for j=1:n
if i~=j
B1(i,j)=-X1(i,j)./(X1(i,j)^2+R1(i,j)^2+eps);
end
end
end
for i=1:n
for j=1:n
if i~=j
B1(i,i)=B1(i,i)-B1(i,j);
end
end
end
disp(B1);
B1=B1(1:4,1:4);
for i=1:n
for j=1:n
if X1(i,j)==0
B2(i,j)=0;
else
if i~=j
B2(i,j)=-1./(X1(i,j)+eps);
end
end
end
end
for i=1:n
for j=1:n
if i~=j
B2(i,i)=B2(i,i)-B2(i,j);
end
end
B2(i,i)=B2(i,i)+YZ(i);
end
disp(B2);
B2=B2(1:3,1:3);
[n1,m1]=size(B2);
V=[1 1 1 1.05 1.05];
cta=[0 0 0 0 0];
t1=1;%cta标识符
t2=1;%电压幅值标识符
S=zeros(1,n-1);% 除平衡节点以外的节点功率
for i=1:n-1
S(i)=A(i,1);
end
P=real(S);%节点的有功功率
Q=imag(S);%节点的无功功率
for k=1:10
for i=1:n
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+V(i)*V(j1)*G(i,j1)*cos(cta(i)-cta(j1))+V(i)*V(j1)*B(i,j1)*s in(cta(i)-cta(j1));
D(i)=D(i)+V(i)*V(j1)*G(i,j1)*sin(cta(i)-cta(j1))-V(i)*V(j1)*B(i,j1)*c os(cta(i)-cta(j1));
end
if i<=4
DP(i)=P(i)-C(i);
end
if A(i,3)==1
DQ(i)=Q(i)-D(i);
end
end
disp(sprintf('第%d次迭代：',k));
disp('DP:');
disp(DP);
disp('DQ:');
disp(DQ);
for i=1:n-1
dpv(i)=DP(i)/V(i);
end;
invB1=inv(B1);
dpcta=invB1*dpv';
disp('----------------------dpcta----------------------'); disp(dpcta');
for i=1:n-1
cta(i)=cta(i)-dpcta(i);
end
disp('-----------------------cta-----------------------'); disp(cta);
dpvv=DQ/B2;
disp('----------------------dpvv-----------------------'); disp(dpvv);
for i=1:3
V(i)=V(i)-dpvv(i);
end
disp('------------------------v-------------------------'); disp(V);
end