基于MATLAB潮流计算软件程序1

基于MATLAB的潮流计算软件程序%求解节点导纳矩阵程序clcclearN=input('请输入节点总数为:N=');nl=input('请输入支路数为：nl=');B1=input('请输入由支路参数形成的矩阵为:B1=');Y=zeros(N);fori=1:nlp=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1./B1(i,3);Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./B1(i,3);Y(p,p)=Y(p,p)+1./B1(i,3);enddisp('导纳矩阵为:Y=');%求出导纳矩阵disp(Y)%开始clccleardisp('节点总数为：');N=5disp('平衡节点为：');1disp('PQ节点为：');JD=[2,3,4,5]Y=[6.25-18.75i-5+15i-1.25+3.75i00;-5+15i10.834-32.5i-1.667+5i-1.667+5i-2.5+7.5i;-1.25+3.75i-1.667+5i12.917-38.75i-10+30i0;0-1.667+5i-10+30i12.917-38.75i-1.25+3.75i;0-2.5+7.5i0-1.25+3.75i3.75-11.25i]disp('导纳矩阵为:Y=');%求出导纳矩阵ei=[1.06;1;1;1;1;];fi=[0;0;0;0;0];P=[0;,0.2;-0.45;-0.40;-0.60];Q=[0;0.20;-0.15;-0.05;-0.10];disp(Y);wucha=0.0001;k=0;G=real(Y);B=imag(Y);N1=4;whilewucha0.00001form=2:5pp=0;qq=0;forn=1:5pp=pp+ei(m)*(G(m,n)*ei(n)-B(m,n)*fi(n))+fi(m)*(G(m,n)*fi(n)+B(m,n)*ei(n));%Pi（0）qq=qq+fi(m)*(G(m,n)*ei(n)-B(m,n)*fi(n))-ei(m)*(G(m,n)*fi(n)+B(m,n)*ei(n));I(m)=(pp-qq*i)/conj(ei(m)+fi(m)*i);%节点注入电流Ii（0）endDP(m)=P(m)-pp;%有功功率的不平衡量DQ(m)=Q(m)-qq;%无功功率的不平衡量H(m,m)=-B(m,m)*ei(m)+G(m,m)*fi(m)+imag(I(m));N(m,m)=G(m,m)*ei(m)+B(m,m)*fi(m)+real(I(m));J(m,m)=-G(m,m)*ei(m)-B(m,m)*fi(m)+real(I(m));L(m,m)=-B(m,m)*ei(m)+G(m,m)*fi(m)-imag(I(m));%公式P131JJ(2*m-3,2*m-3)=H(m,m);JJ(2*m-3,2*m-2)=N(m,m);JJ(2*m-2,2*m-3)=J(m,m);JJ(2*m-2,2*m-2)=L(m,m);endform=2:Nforn=2:Nifn==melseH(m,n)=-B(m,n)*ei(m)+G(m,n)*fi(m);N(m,n)=G(m,n)*ei(m)+B(m,n)*fi(m);J(m,n)=-B(m,n)*fi(m)-G(m,n)*ei(m);L(m,n)=G(m,n)*fi(m)-B(m,n)*ei(m);JJ(2*m-3,2*n-3)=H(m,n);JJ(2*m-3,2*n-2)=N(m,n);JJ(2*m-2,2*n-3)=J(m,n);JJ(2*m-2,2*n-2)=L(m,n);%公式P130endendenddisp('雅克比矩阵JJ:');disp(JJ);form=2:NDW(2*m-3)=DP(m);DW(2*m-2)=DQ(m);%形成DP和DQ不平衡量的矩阵endDY=JJ\DW';%解修正方程式wucha=max(abs(DY));disp('第M次修正方程的解DY:');disp(DY);forn=2:Nfi(n)=fi(n)+DY(2*n-3);ei(n)=ei(n)+DY(2*n-2);%计算各节点电压新值endu=ei+fi*i;disp('节点电压的第C(k)次近似值:');disp(u);disp('各点的电压实部ei(单位:V)为(节点号从小到大排列):');disp(ei);disp('各点的电压虚部fi(单位:V)为(节点号从小到大排列):');disp(fi);k=k+1;disp('迭代次数:');disp(k);endforn=1:N1+1I(n)=Y(1,n)*u(n);%(公式P139)enddisp('平衡节点的功率');S1=u(1)*sum(conj(I))%计算平衡节点的功率form=1:N1+1forn=1:N1+1S(m,n)=u(m)*(conj(u(m))-conj(u(n)))*conj(-Y(m,n));%计算m-n各支路功率S2(n,m)=u(n)*(conj(u(n))-conj(u(m)))*conj(-Y(n,m));%计算n-m各支路功率DS(m,n)=S(m,n)+S2(n,m);%计算各支路功率的损耗(公式P1354-51a)endenddisp('各支路m-n功率');disp(S)disp('各支路n-m功率');disp(S2)disp('各支路功率损耗');disp(DS)%结束