电力系统潮流计算程序(详细)

潮流计算，电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。

潮流计算程序，相关的原始数据数据数据输入格式如下：

%B1是支路参数矩阵，第一列和第二列是起始节点编号和终点节点编号

%第三列、第四列、第五列、第六列、第七列、第八列分别为：支路电阻、电抗、电导、电纳、变压器变比、是否有变压器（1为有、0为无）。

%B2为节点参数矩阵，其中第一列到第六列为节点编号；为节点类型；注入有功、注入无功、电压幅值、电压相位。

%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

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

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

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

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

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

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

Y=zeros(n);

N=1;

%建立节点导纳矩阵

for i=1:n1

if B1(i,8)==0

p=B1(i,1);

q=B1(i,2);

Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4))-B1(i,5);

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

Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);

Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);

else

p=B1(i,1);

q=B1(i,2);

Y(p,q)=Y(p,q)-1/((B1(i,3)+B1(i,4))*B1(i,7)) -B1(i,5);

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

Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+B1(i,5);

Y(q,q)=Y(q,q)+1/(B1(i,7)^2*(B1(i,3)+B1(i,4)))+B1(i,5);

end

end

Y

G=real(Y);

B=imag(Y);

PriS=zeros(2*n-2,1);

ImbS=zeros(2*n-2,1);

%创建PriS，用于存储初始功率参数

h=0;

j=0;

for i=1:n

if i~=isb&B2(i,2)==1

h=h+1;

for j=1:n

PriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i, j)*B2(j,5));

PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B 2(j,5));

end

end

end

for i=1:n

if i~=isb&B2(i,2)==2

h=h+1;

for j=1:n

PriS(2*h-1,1)=PriS(2*h-1,1)+

B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));

PriS(2*h,1)=PriS(2*h,1)+

B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5))

end

end

end

PriS

U3=zeros(n-h-1,2);%U3存储PV节点的初始电压

t=0;

for i=1:n

if B2(i,2)==2

t=t+1;

U3(t,1)=B2(i,5);

U3(t,2)=B2(i,6);

end

end

U3

%ImbS于存储有功功率、无功功率和电压幅值的不平衡量

h=0;

for i=1:n

if i~=isb&B2(i,2)==1

h=h+1;

ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);

ImbS(2*h,1)=B2(i,4)-PriS(2*h,1);

end

end

t=0;

for i=1:n

if i~=isb&B2(i,2)==2

h=h+1;

t=t+1;

ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);

ImbS(2*h,1)=U3(t,1)^2+U3(t,2)^2-B2(i,5)^2-B2(i,6)^2;

end

end

ImbS

I=zeros(n-1,1);%I，存储节点电流参数

h=0;

for i=1:n

if i~=isb

h=h+1;

I(h,1)=(PriS(2*h-1,1)-PriS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));

end

end

I

Jacbi=zeros(2*n-2);%Jacbi(雅可比矩阵)

h=0;

k=0;

for i=1:n

if B2(i,2)==1

h=h+1;

for j=1:n

if j~=isb

k=k+1;

if i==j

Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));

Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));

else

Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);

Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);

end