潮流计算的计算机算法

高等电力系统分析

（潮流计算的计算机算法）PQ分解法潮流计算（IEEE14)

目录

一、MATLAB源程序

二、对支路参数（B1）、节点参数（B2）的说明

三、带入数据，运行结果

一、MATLAB源程序

clear

close all

n=input('请输入节点数:n=');

n1=input('请输入支路数:n1=');

isb=input('请输入平衡节点号:isb=');

pr=input('请输入误差精度:pr=');

B1=input('请输入支路参数:B1=');

B2=input('请输入节点参数:B2=');

n2=input('请输入PQ节点个数:n2=');

Y=zeros(n);

for i=1:n1

p=B1(i,1);

q=B1(i,2);

Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4)*1j); %非对角元

Y(q,p)=Y(p,q);

Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; %对角元 Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j;

end

disp('导纳矩阵Y=');

disp(Y)

%---------------------------------------------

%---------------下面是求P,Q,V,O矩阵---------------

V=zeros(1,n);O=zeros(1,n);P=zeros(1,n);Q=zeros(1,n);

G=real(Y);B=imag(Y);

for i=1:n

P(i)=B2(i,3);

Q(i)=B2(i,4);

V(i)=B2(i,5);

O(i)=B2(i,6);

end

B3=B(1:n-1,1:n-1); %不含平衡节点，由节点导纳虚部构成

B4=B(1:n2,1:n2); %所有PQ节点

%----------------------------------------------

%---------------下面是求ΔP,ΔQ矩阵---------------

DX=0;ICT=1;Mp=1;Mq=1;

while ICT~=0

m1=1;m2=1;

for i=1:n

if i~=isb

C(i)=0;

D(i)=0;

for j1=1:n

C(i)=C(i)+V(i)*V(j1)*(G(i,j1)*cos(O(i)-O(j1))+B(i,j1)*sin(O(i)-O( j1)));

D(i)=D(i)+V(i)*V(j1)*(G(i,j1)*sin(O(i)-O(j1))-B(i,j1)*cos(O(i)-O( j1)));

end

DP(m1)=P(i)-C(i);

m1=m1+1;

if B2(i,2)==1

DQ(m2)=Q(i)-D(i);

m2=m2+1;

end

end

end

m1=m1-1; %所有节点数

m2=m2-1; %PQ节点数

DPQ=[DP';DQ']; %求DP,DQ

V1=V(:,1:m1);

V2=diag(V1);

V3=inv(V2); %对V矩阵求逆

H=V3*DP'; %ΔP/V

K=-inv(B3)*H; %-ΔP/V/B3

deltO=V3*K; %Δ角=-ΔP/V/V/B3

max1=max(abs(DP));

for i=1:m1

if max1

Mp=0;

else

O(i)=O(i)+deltO(i)';

Mq=1;

end

end

V4=V(:,1:m2);

V5=diag(V4);

V6=inv(V5);

L=V6*DQ';

N=-inv(B4)*L;

deltV=N; %ΔV=-ΔQ/V/B

max2=max(abs(DQ));

for i=1:m2

if max2

Mq=0;

else

if B2(i,2)==1;

V(i)=V(i)+deltV(i)';

Mp=1;

end

end

end

if Mp==0&&Mq==0

ICT=0;

else

ICT=1;

end

DX=DX+1;

end

%-------------------------------------------

%----------------迭代结束，开始输出结果----------------

disp('----------------------------------------------'); disp('迭代次数为：');

disp(DX);

for i=1:n

E(i)=V(i)*cos(O(i))+1j*V(i)*sin(O(i));

o(i)= 180*angle(E(i))/pi;

end

disp('---------------------------------------------'); disp('修正后各节点电压标么值为（节点号从小到大排列）:');

disp(V);

disp('---------------------------------------------'); disp('修正后各节点电压相角为（节点号从小到大排列）:');

disp(o);

%------------计算各个节点的功率----------------------

disp('----------------------------------------------'); disp('各节点的功率为:');

for p=1:n

C(p)=0;

for q=1:n

C(p)=C(p)+conj(Y(p,q)*conj(E(q)));

end

S(p)=E(p)*C(p);

end