三机九节点潮流计算

clc;clear;closeall;n1=9;n2=9;B1=load('b1.txt','ascii');B2=load('b2.txt','ascii');m=0;r=0;fori=1:n1ifB1(i,8)==1m=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);forj=1:n2p=B2(j,2);q=B2(j,3);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(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))fori=1:n1forj=1:n1si(i,j)=sin(f(i)-f(j));co(i,j)=cos(f(i)-f(j));ifi==jH1(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);elseH1(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));endendendH=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=[HN;JL];%雅克比矩阵Jfori=2:n1a(i)=0;b(i)=0;forj=1:n1a(i)=a(i)+e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));ifi=(r+2)b(i)=b(i)+e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endenddp(i-1)=ps(i)-e(i)*a(i);ifi=(r+2)dq(i-(r+1))=qs(i)-e(i)*b(i);endendF=[dp';dq'];dx=-J\F;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('各节点电压');disp('e=');disp(e');disp('单位：度f=');disp(f'*180/3.1415926)%求出各节点电压幅值及相角a=0;fori=1:n1a=a+conj(Y(1,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i)));enda=a*(e(1)*cos(f(1))+1i*e(1)*sin(f(1)));disp('平衡节点P+jQ=');disp(a);disp('输电线路有功无功:');fori=1:n2p=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(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%初值计算forj=1:n1s(j)=0;fori=1:n1s(j)=s(j)+conj(Y(j,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i)));ends(j)=s(j)*(e(j)*cos(f(j))+1i*e(j)*sin(f(j)));endY1=conj(-s)./(e'.^2);%负荷等值导纳disp('各节点等值导纳');disp('Y1=');disp(Y1);disp('各节点功率');disp('S=');disp(s);Ra=[000];X_d=[0.06080.11980.1813];Xq=[0.09690.86451.2578];%Xq=X_d;%不计凸极效应v=e.*cos(f)+1i*e.*sin(f);fori=1:3I(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));enddisp('各发电机的暂态电动势，功角和输入机械功率初值');disp('delta=');disp(delta*180/3.1415926);fori=1:3Vx(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);E_q(i)=Vq(i)+Ra(i)*Iq(i)+X_d(i)*Id(i);enddisp('E_q=');disp(E_q);fori=1:3Pm(i)=real(s(i))+(Ix(i)^2+Iy(i)^2)*Ra(i);enddisp('Pm=');disp(Pm);