《matlab潮流计算》word版.doc

%线路数据 Line=1-节点编号 2-线路首端节点号 3-线路末端节点号 4-支路电阻 5-支路电抗 6-支路电纳(注意：此处取的是B/2)Line=1 1 2 0.04 0.25 0.25 2 1 3 0.1 0.35 0 3 2 3 0.08 0.30 0.25 ; %变压器数据 transform=1-支路编号 2-支路首节点编号 3-支路末节点编号 4-支路电阻(p.u.) 5-支路电抗(p.u.) 6-变压器变比(p.u.)transform= 1 2 4 0 0.015 1.052 3 5 0 0.03 1.05;%数据预处理 Nbus=5; % 节点数 nline=size(Line,1); % 线路个数 ntrans=size(transform,1); % 变压器个数 slack=5; % 平衡节点号 Npq=4; % PQ节点的个数 Ss=0;%计算节点导纳矩阵Y=zeros(Nbus);if nline=1%判断是否存在线路 for k=1:nline %以下处理线路 t1=Line(k,2); t2=Line(k,3); b2=Line(k,6);%分别取出线路的首端节点编号t1、末端节点编号t2和对地电纳b2 Yl=1/(Line(k,4)+j*Line(k,5);%计算线路的支路电导Yl Y(t1,t1)=Y(t1,t1)+Yl+j*b2;%修正第k条线路首端节点的自导纳 Y(t1,t2)=Y(t1,t2)-Yl;%修正第k条线路首端节点与末端节点之间的互导纳 Y(t2,t1)=Y(t2,t1)-Yl;%修正第k条线路末端节点与首端节点之间的互导纳 Y(t2,t2)=Y(t2,t2)+Yl+j*b2;%修正第k条线路末端节点的自导纳 endendif ntrans=1%判断是否存在变压器 for k=1:ntrans %以下处理变压器 t1=transform(k,2);t2=transform(k,3);t3=transform(k,6);%分别取出变压器的首端节点编号t1、末端节点编号t2和变比t3 Yt=1/(transform(k,4)+j*transform(k,5); Yt1=Yt/t3; Yt2=Yt*(1-t3)/(t3*t3); Yt3=Yt*(t3-1)/t3; Y(t1,t1)=Y(t1,t1)+Yt1+Yt2; Y(t1,t2)=Y(t1,t2)-Yt1; Y(t2,t1)=Y(t2,t1)-Yt1; Y(t2,t2)=Y(t2,t2)+Yt1+Yt3; endendG=real(Y);B=imag(Y); %区分节点导纳矩阵的实部和虚部GB%赋初值delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;u(1)=1;u(2)=1;u(3)=1;u(4)=1;p(1)=-0.30;q(1)=-0.18;p(2)=-0.55;q(2)=-0.13;p(3)=0;q(3)=0;p(4)=0.5;q(4)=1.10;p(5)=0.8;q(5)=0.50;k=0;precision=1;Npq=4; %Npq分别是网络中的PQ节点数%Unbalance=-JacobiCorrectionwhile precision0.00001 %设定误差上限，判断是否继续迭代 u(5)=1.06; delt(5)=0; %设定平衡节点电压相角与幅值 k; u; delt; for m=1:Npq for n=1:Nbus pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n)+B(m,n)*sin(delt(m)-delt(n); %由节点电压求得的PQ节点注入有功功率 qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n)-B(m,n)*cos(delt(m)-delt(n); %由节点电压求得的PQ节点注入无功功率 end Unbalance(2*m-1)=p(m)-sum(pt); %计算PQ节点有功功率不平衡量 Unbalance(2*m )=q(m)-sum(qt); %计算PQ节点无功功率不平衡量end %Unbalance是节点不平衡量矩阵 for m=1:Npq for n=1:Nbus h0(n)= u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n)-B(m,n)*cos(delt(m)-delt(n); n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n)+B(m,n)*sin(delt(m)-delt(n); j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n)+B(m,n)*sin(delt(m)-delt(n); l0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n)-B(m,n)*cos(delt(m)-delt(n); end H(m,m)=sum(h0)-u(m)2*(G(m,m)*sin(delt(m)-delt(m)-B(m,m)*cos(delt(m)-delt(m); N(m,m)=sum(n0)+u(m)2*(G(m,m)*cos(delt(m)-delt(m)+B(m,m)*sin(delt(m)-delt(m)-2*u(m)2*G(m,m); J(m,m)=sum(j0)+u(m)2*(G(m,m)*cos(delt(m)-delt(m)+B(m,m)*sin(delt(m)-delt(m); L(m,m)=sum(l0)+u(m)2*(G(m,m)*sin(delt(m)-delt(m)-B(m,m)*cos(delt(m)-delt(m)+2*u(m)2*B(m,m); Jacobi(2*m-1,2*m-1)=H(m,m); Jacobi(2*m-1,2*m )=N(m,m); Jacobi(2*m ,2*m-1)=J(m,m); Jacobi(2*m ,2*m )=L(m,m); end %计算m=n情况下的Jacobi矩阵中的子矩阵元素 for m=1:Npq for n=1:Npq if m=n else H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n)-B(m,n)*cos(delt(m)-delt(n); J(m,n)= u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n)+B(m,n)*sin(delt(m)-delt(n); N(m,n)=-J(m,n); L(m,n)= H(m,n); Jacobi(2*m-1,2*n-1)=H(m,n); Jacobi(2*m-1,2*n )=N(m,n); Jacobi(2*m ,2*n-1)=J(m,n); Jacobi(2*m ,2*n )=L(m,n); end end end %计算mn情况下的Jacobi矩阵中的子矩阵元素 Correction=-Jacobi(Unbalance); %计算电压相角和幅值的修正量 precision = max(abs(Correction); %取误差最大值 for m=1:Npq delt(m)=delt(m)+Correction(2*m-1); %修正PQ节点电压相角 u(m)= u(m)+Correction(2*m ); %修正PQ节点电压幅值 end k=k+1; %迭代轮数+1endk, u, delt, Jacobi, precisionfor m=1:Nbus U(m)=u(m)*(cos(delt(m)+j*sin(delt(m); %采用直角坐标系表示电压 I(m)=Y(Nbus,m)*U(m); %计算注入平衡节点的电流endSslack=U(Nbus)*sum(conj(I) %计算注入平衡节点的功率for m=1:Nbus for n=1:Nbus S(m,n)=U(m)*(conj(U(m)-conj(U(n)*conj(-