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

%本程序的功能是用牛拉法进行含风电场的电力系统潮流计算
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
end
%求雅克比矩阵
for k=3:N0
k1=k+1;N1=N;
for k2=k1:N1
J(k,k2)=J(k,k2)./J(k,k);
end
J(k,k)=1;
if k~=3;
k4=k-1;
for k3=3:k4
for k2=k1:N1
J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2);
end
J(k3,k)=0;
end
if k==N0,break;end
for k3=k1:N0
for k2=k1:N1
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);
end
J(k3,k)=0;
end
else
for k3=k1:N0
for k2=k1:N1
J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2);
end
J(k3,k)=0;
end
end
end
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
k1=k+1;
f(L)=f(L)-J(k1,N);
end
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
IT2=IT2+1;
end
end
ICT2(a)=IT2;
ICT1=ICT1+1;
end
%用高斯消去法解w=-J*V
fid1=fopen('out1.txt','wt')
fprintf(fid1,‘各节点的实际电压标幺值E和电压大小V为：\n’);
for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
fprintf(fid1,'E(%d)=%8.5f+j%8.5f, V(%d)=%8.5f\n',k,e(k),f(k),k,V(k)); end
for p=1:n
if p==isb
fprintf(fid1,'平衡节点的功率S为:\n');
C(p)=0;D(p)=0;
for q=1:n
C(p)=C(p)+G(p,q)*e(q)-B(p,q)*f(q);
D(p)=D(p)-G(p,q)*f(q)-B(p,q)*e(q);
end
Sp(p)=e(p)*C(p)-f(p)*D(p);
Sq(p)=(e(p)*D(p)+f(p)*C(p));
fprintf(fid1,'S=%8.5f+j%8.5f\n',Sp(p),Sq(p));
end
end
fprintf(fid1,'各条支路的首端功率为:\n');
for i=1:nl
p=B1(i,1);q=B1(i,2);
sip(i)=-e(p)*B1(i,5)+(e(q)-e(p))*G(p,q)-(f(q)-f(p))*B(p,q);
siq(i)=(f(p)*B1(i,5)+(e(p)-e(q))*B(p,q)-(f(q)-f(p))*G(p,q));
Sip(i)=e(p)*sip(i)-f(p)*siq(i);
Siq(i)=e(p)*siq(i)+f(p)*sip(i);
fprintf(fid1,'S(%d,%d)=%8.5f+j%8.5f\n',p,q,Sip(i),Siq(i));
end
fprintf(fid1,'¸各条支路的末端功率为:\n');
for i=1:nl
p=B1(i,1);q=B1(i,2);
sjp(i)=-e(q)*B1(i,6)+(e(p)-e(q))*G(p,q)-(f(p)-f(q))*B(p,q);
sjq(i)=(f(q)*B1(i,6)+(e(q)-e(p))*B(p,q)-(f(p)-f(q))*G(p,q));
Sjp(i)=e(q)*sjp(i)-f(q)*sjq(i);
Sjq(i)=e(q)*sjq(i)+f(q)*sjp(i);
fprintf(fid1,'S(%d,%d)=%8.5f+j%8.5f\n',q,p,Sjp(i),Sjq(i));
end
fprintf(fid1,'各条支路的功率损耗为:\n');
for i=1:nl
DSp(i)=Sip(i)+Sjp(i);
DSq(i)=Siq(i)+Sjq(i);
fprintf(fid1,'DS(%d)=%8.5f+j%8.5f\n',i,DSp(i),DSq(i)); end
fclose(fid1);。