三机九节点潮流计算

1、精选优质文档-倾情为你奉上clc;clear;close all;n1=9;n2=9;B1=load('b1.txt','ascii');B2=load('b2.txt','ascii');m=0;r=0;for i=1:n1 if B1(i,8)=1 m=m+1; endendr=n1-m-1;e=B1(1:n1,2);f=B1(1:n1,3);ps=B1(1:n1,4)-B1(1:n1,6);qs=B1(1:n1,5)-B1(1:n1,7);Y=zeros(n1);for j=1:n2 p=B2(j,2);q=B2(j,3);

2、J=zeros(2*n1); Y(p,q)=-1./(B2(j,4)+1i*B2(j,5)*B2(j,7); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/(B2(j,4)+1i*B2(j,5)+1i*B2(j,6); Y(q,q)=Y(q,q)+1/(B2(j,4)+1i*B2(j,5)*B2(j,7)2)+1i*B2(j,6);end%导纳矩阵YH1=zeros(n1);N1=zeros(n1);J1=zeros(n1);L1=zeros(n1);G=real(Y);B=imag(Y);dp(1:(n1-1)=1;dq(1:m)=1;dx(1:(2*m+r)=1;while

3、 (abs(max(dp)>=1.0e-5)|abs(max(dq)>=1.0e-5)|abs(max(dx(1:(m+r)>=1.0e-5)|abs(max(dx(m+r+1):(2*m+r).*e(1:m)>=1.0e-5)for i=1:n1 for j=1:n1 si(i,j)=sin(f(i)-f(j); co(i,j)=cos(f(i)-f(j); if i=j H1(i,j)=e(i)2*B(i,i)+qs(i); N1(i,j)=-e(i)2*G(i,i)-ps(i); J1(i,j)=e(i)2*G(i,i)-ps(i); L1(i,j)=e(i)2*

4、B(i,i)-qs(i); else H1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j); N1(i,j)=-e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j); J1(i,j)=e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j); L1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j); end endendH=H1(2:n1,2:n1);N=N1(2:n1,r+2:n1);J=J1(r+2:n1,2:n1);L=L1(r+2:n1,r+2:

5、n1);J=H N;J L;%雅克比矩阵J for i=2:n1 a(i)=0;b(i)=0; for j=1:n1 a(i)=a(i)+e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j); if i>=(r+2) b(i)=b(i)+e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j); end end dp(i-1)=ps(i)-e(i)*a(i); if i>=(r+2) dq(i-(r+1)=qs(i)-e(i)*b(i); endendF=dp'dq'dx=-JF;f(2:n1)=f(2:n1)+dx(1:(m+r);

6、e1=e(r+2):n1)+dx(m+r+1):(2*m+r).*e(r+2):n1);e2=e(1:(r+1);e=e2;e1;enddisp('各节点电压');disp('e=');disp(e');disp('单位：度 f=');disp(f'*180/3.)%求出各节点电压幅值及相角 a=0; for i=1:n1 a=a+conj(Y(1,i)*(e(i)*cos(f(i)-1i*e(i)*sin(f(i); end a=a*(e(1)*cos(f(1)+1i*e(1)*sin(f(1);disp('平衡节点P+

7、jQ=');disp(a);disp('输电线路有功无功:');for i=1:n2 p=B2(i,2);q=B2(i,3); s1(p,q)=e(p)2*(-1i*B2(i,6)+(e(p)*cos(f(p)+1i*e(p)*sin(f(p)*(e(p)*cos(f(p)-1i*e(p)*sin(f(p)-e(q)*cos(f(q)+1i*e(q)*sin(f(q)*conj(-Y(q,p); s1(q,p)=e(q)2*(-1i*B2(i,6)+(e(q)*cos(f(q)+1i*e(q)*sin(f(q)*(e(q)*cos(f(q)-1i*e(q)*sin(f(

8、q)-e(p)*cos(f(p)+1i*e(p)*sin(f(p)*conj(-Y(q,p); fprintf('S%d%d=',p,q);disp(s1(p,q); fprintf('S%d%d=',q,p);disp(s1(q,p);end%初值计算 for j=1:n1 s(j)=0; for i=1:n1 s(j)=s(j)+conj(Y(j,i)*(e(i)*cos(f(i)-1i*e(i)*sin(f(i); end s(j)=s(j)*(e(j)*cos(f(j)+1i*e(j)*sin(f(j); end Y1=conj(-s)./(e'

9、.2);%负荷等值导纳 disp('各节点等值导纳'); disp('Y1=');disp(Y1); disp('各节点功率'); disp('S=');disp(s); Ra=0 0 0;X_d=0.0608 0.1198 0.1813;Xq=0.0969 0.8645 1.2578; %Xq=X_d;%不计凸极效应 v=e.*cos(f)+1i*e.*sin(f); for i=1:3 I(i)=conj(s(i)/(e(i)*cos(f(i)+1i*e(i)*sin(f(i); EQ(i)=e(i)*cos(f(i)+1i*

10、e(i)*sin(f(i)+(Ra(i)+1i*Xq(i)*I(i); EQx(i)=real(EQ(i);EQy(i)=imag(EQ(i); delta(i)=atan(EQy(i)/EQx(i); end disp('各发电机的暂态电动势，功角和输入机械功率初值'); disp('delta=');disp(delta*180/3.); for i=1:3 Vx(i)=e(i)*cos(f(i);Vy(i)=e(i)*sin(f(i); Vdq=sin(delta(i) -cos(delta(i);cos(delta(i) sin(delta(i)*Vx(i);Vy(i); Vd(i)=Vdq(1);Vq(i)=Vdq(2); Ix(i)=real(I(i);Iy(i)=imag(I(i); Idq=sin(delta(i) -cos(delta(i);cos(delta(i) sin(delta(i)*Ix(i);Iy(i); Id(i)=Idq(1);Iq(i)=Idq(2);