含风电场的电力系统潮流计算程序

%本程序的功能是用牛拉法进行含风电场的电力系统潮流计算

function s1=pf(a,B1,B2)

n=a(1);%节点数

nl=a(2);%支路数

isb=a(3);%平衡节点号

pr=a(4);%误差精度

for i=1:n

for j=1:n

G(i,j)=0;

B(i,j)=0;

end

end

%求导纳矩阵

%B1为支路参数矩阵，其每一行格式为首节点号，末节点号，支路电导，支路电纳，首节点对地电纳，末节点对地电纳

for i=1:nl

p=B1(i,1);q=B1(i,2)

G(p,q)=G(p,q)-B1(i,3);

B(p,q)=B(p,q)- B1(i,4);

G(q,p)=G(p,q);

B(q,p)=B(p,q);

G(p,p)=G(p,p)+ B1(i,3);

B(p,p)=B(p,p)+B1(i,4);

G(q,q)=G(q,q)+ B1(i,3);

B(q,q)=B(q,q)+B1(i,4);

end

for i=1:n

B(i,i)=B(i,i)+B2(i,5);

end

%求导纳矩阵

%B2为节点参数矩阵，每一行格式为节点注入有功，注入无功，电压实部，电压虚部，对地电纳，节点类型

%节点类型：1为平衡节点，2为PQ节点，3为PV节点，4为风电场节点

%Bf为风电场参数，格式为：有功功率，定子电抗，转子漏抗，转子电阻，励磁电抗

for i=1:n

if B2(i,6)==4

p(i)=Bf(1);

a1=2*p(i)^2*(Bf(2)+Bf(3))^2;

a2=-Bf(4)^3*Bf(5);

a3=Bf(4)^2*Bf(5);

a4=Bf(4)^2;

a5=4*p(i)^2*(Bf(2)+Bf(3))^2*Bf(4)^2;

else

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

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

end

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

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

V(i)=sqrt(e(i)^2+f(i)^2);

end

ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;

while IT2~=0

IT2=0;a=a+1;

for i=1:n

if i~=isb

C(i)=0;

D(i)=0;

for j1=1:n

C(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);

D(i)= D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);

end

P1=C(i)*e(i)+f(i)*D(i);

Q1=f(i)*C(i)-D(i)*e(i);

V2=e(i)^2+f(i)^2;

if B2(i,6)==2

DP=P(i)-P1;

DQ=Q(i)-Q1;

for j1=1:n

if j1~=isb&j1~=i

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

X3=X2;

X4=-X1;

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;

elseif j1==i&j1~=isb

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);

X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);

X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;

end

end

elseif B2(i,6)==3

DP=P(i)-P1;

DV=V(i)^2-V2;

for j1=1:n

if j1~=isb&j1~=i

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

X5=0;

X6=0;

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;

elseif j1==i&j1~=isb

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);

X5=-2*e(i);

X6=-2*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

J(m,q)=X2;

end

end

else

DP=P(i)-P1;

for j1=1:n

if j1~=isb&j1~=i

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

X3=X2;

X4=-X1;

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;

elseif j1==i&j1~=isb

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);

x3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);

x4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);

x5=2*a1*(a2+a3*a4*V2^2/(sqrt(V2^4*a4-a5)))/(a2*V2^2+a3*(sqrt(V2^4*a4-a5)))^2;

X3=x3-e(i)*x5;

X4=x4-f(i)*x5;

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;

end

end

end

end