电力系统数字仿真作业1

电力系统数字仿真作业

电力系统及其自动化研硕10-13 韩暘

1 . Matlab潮流计算作业

本程序的功能是用牛顿——拉夫逊法进行潮流计算。其中：

B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳；5、支路的变比；6、支路首端处于K 侧为1，1侧为0

B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值；4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量；6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；

3为PV节点；

程序如下：

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

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

isb=input('请输入平衡母线节电号:isb=');

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

B1=input('请输入由支路参数形成的矩阵:B1=');%变压器侧为1，否则为0

B2=input('请输入各节点参数形成的矩阵:B2=');

Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);

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

for i=1:nl %支路数

if B1(i,6)==0 %左节点处于1侧

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

else %左节点处于K侧

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

end

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

Y(q,p)=Y(p,q); %非对角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧

end

%求导纳矩阵

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

disp(Y)

e(1)=1.00;e(2)=1.00;e(3)=1.10;e(4)=1.05;

f(1)=0;f(2)=0;f(3)=0;f(4)=0;

G=real(Y);

B=imag(Y);

%设S=P+Qj；

S(1)=-0.30-0.18j;

S(2)=-0.55-0.13j;

S(3)=0.50;

P=real(S);

Q=imag(S);

V(3)=1.10;

V(4)=1.05+0j;

k=0;

precision=1;

N1=3;

while precision > 0.00001

for m=1:N1

for n=1:4

Pt(n)=(e(m)*(G(m,n)*e(n)-B(m,n)*f(n))+f(m)*(G(m,n)*f(n)+B(m,n)*e(n))); end

dP(m)=P(m)-sum(Pt);

end

for m=1:3

for n=1:4

Qt(n)=(f(m)*(G(m,n)*e(n)-B(m,n)*f(n))-e(m)*(G(m,n)*f(n)+B(m,n)*e(n))); end

dQ(m)=Q(m)-sum(Qt);

end

Vt=e(3)*e(3)+f(3)*f(3);

dV=(V(3)^2)-Vt;

for m=1:3

for n=1:4

Bi(n)=G(m,n)*f(n)+B(m,n)*e(n);

Ai(n)=G(m,n)*e(n)-B(m,n)*f(n);

end

end

for m=1:3

for n=1:3

if m==n

H(m,n)=sum(Ai)-G(m,n)*e(m)-B(m,n)*f(m);

N(m,n)=-sum(Bi)+B(m,n)*e(m)-G(m,n)*f(m);

J(m,n)=sum(Bi)+B(m,n)*e(m)-G(m,n)*f(m);

L(m,n)=-sum(Ai)+G(m,n)*e(m)+B(m,n)*f(m);

J(3,3)=-2*e(3);

L(3,3)=-2*f(3);

else

H(m,n)=-G(m,n)*e(m)-B(m,n)*f(m);

N(m,n)=B(m,n)*e(m)-G(m,n)*f(m);

J(m,n)=B(m,n)*e(m)-G(m,n)*f(m);

L(m,n)=G(m,n)*e(m)+B(m,n)*f(m);

J(3,1)=0;

L(3,1)=0;

J(3,2)=0;

L(3,2)=0;

end

end

end

%生成雅克比矩阵。

%对角线上的元素

for m=1:3

Ja(2*m-1,2*m-1)=H(m,m);

Ja(2*m-1,2*m)=N(m,m);

Ja(2*m,2*m-1)=J(m,m);

Ja(2*m,2*m)=L(m,m);

end

%非对角线上的元素

for m=1:3

for n=1:3

if m==n

else

Ja(2*m-1,2*n-1)=H(m,n);

Ja(2*m-1,2*n)=N(m,n);

Ja(2*m,2*n-1)=J(m,n);

Ja(2*m,2*n)=L(m,n);

end

end

end

PQ=[dP(1);dQ(1);dP(2);dQ(2);dP(3);dV];

dU=-inv(Ja)*PQ;

precision=max(abs(dU));

for n=1:3

e(n)=e(n)+dU(2*n-1);

f(n)=f(n)+dU(2*n);

end

k=k+1;

for n=1:4

Un(n)=e(n)+f(n)*j;

end

disp('第')

k-1

disp('次迭代')

disp('各节点电压')

Un

disp('迭代过程中节点不平衡量的变化情况') PQ

end

disp('计算平衡节点的功率')

for m=1:4

In(m)=Y(m,4)*Un(m);

end

S=Un(4)*sum(conj(In))