三机九节点潮流计算

《三机九节点潮流计算》由会员分享，可在线阅读，更多相关《三机九节点潮流计算（4页珍藏版）》请在装配图网上搜索。

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);J=zeros(2*n1); Y(p,q)=-1./(B2(j,4)+1i*B2

2、(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 (abs(max(dp)=1.0e-5)|abs(max(dq)=1.0e-5

3、)|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*B(i,i)-qs(i); else H1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(

4、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:n1);J=H N;J L;%雅克比矩阵J for i=2:n1 a(i)=0;b(i)=0; for j=1:

5、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);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

6、(各节点电压);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+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(

7、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(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:

8、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.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.

9、*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*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

10、);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); E_q(i)=Vq(i)+Ra(i)*Iq(i)+X_d(i)*Id(i); end disp(E_q=);disp(E_q); for i=1:3 Pm(i)=real(s(i)+(Ix(i)2+Iy(i)2)*Ra(i); end disp(Pm=);disp(Pm);专心-专注-专业