MATLAB潮流程序(IEEE14-直角坐标-牛拉法)

MATLAB潮流程序（IEEE14直角坐标牛拉法）clearbaseMVA=100;%功率基值%%读Data1中数据loadData1.txtBus=Data1(:,1);%节点号Vtype=Data1(:,5);%节点类型Pload=Data1(:,8);%负载有功Qload=Data1(:,9);%负载无功Pgen=Data1(:,10);%发电机发出有功Qgen=Data1(:,11);%发电机发出无功Vset=Data1(:,13);%电压设定点Qsh=Data1(:,17);%并联电容电纳标幺值%%读Data2中数据loadData2.txtII=Data2(:,1);JJ=Data2(:,2);%支路端点号Ltype=Data2(:,5);%线路类型R=Data2(:,6);%两点间电阻X=Data2(:,7);%两点间电抗B=Data2(:,8)/2;%线路对地电纳K=Data2(:,14);%变压器非标准变压比%%求导纳矩阵Yy1=zeros(14);y2=zeros(14);y3=zeros(14);lin=length(II);%支路数forx=1:linswitchLtype(x)case1y1(II(x),JJ(x))=1/(R(x)+i*X(x));y1(JJ(x),II(x))=y1(II(x),JJ(x));y3(II(x),JJ(x))=i*B(x);y3(JJ(x),II(x))=i*B(x);case2y1(II(x),JJ(x))=1/((R(x)+i*X(x))*K(x));y1(JJ(x),II(x))=y1(II(x),JJ(x));y2(II(x),JJ(x))=(1-K(x))/((R(x)+i*X(x))*K(x)^2);y2(JJ(x),II(x))=(K(x)-1)/((R(x)+i*X(x))*K(x));endendclearxY=zeros(14);forx=1:14Y(x,x)=sum(y1(x,:))+sum(y2(x,:))+sum(y3(x,:))+i*Qsh(x);endclearx;Y=Y-y1;G=real(Y);B=imag(Y);%%设电压初值U=Vset;e=real(U);f=imag(U);%%Ps=zeros(1,14);Qs=zeros(1,14);D=ones(26,1);forx=1:14Ps(x)=(Pgen(x)-Pload(x))/baseMVA;Qs(x)=(Qgen(x)-Qload(x))/baseMVA;endclearx;N=0;Jacbi=zeros(26);whilemax(abs(D))0.000001forx=2:14%节点功率及电压不平衡量switchVtype(x)case1%PQ节点D(2*x-3)=Ps(x)-e(x)*(G(x,:)*e-B(x,:)*f)-f(x)*(G(x,:)*f+B(x,:)*e);D(2*x-2)=Qs(x)-f(x)*(G(x,:)*e-B(x,:)*f)+e(x)*(G(x,:)*f+B(x,:)*e);case2%PV节点D(2*x-3)=Ps(x)-e(x)*(G(x,:)*e-B(x,:)*f)-f(x)*(G(x,:)*f+B(x,:)*e);D(2*x-2)=Vset(x).*Vset(x)-(e(x).^2+f(x).^2);endendclearmforI=2:14%求雅克比矩阵forJ=2:14ifI~=J%非对角元素Jacbi((2*I-3),(2*J-3))=B(I,J)*e(I)-G(I,J)*f(I);Jacbi((2*I-3),(2*J-2))=-(G(I,J)*e(I)+B(I,J)*f(I));switchVtype(I)case1%PQ节点Jacbi((2*I-2),(2*J-3))=G(I,J)*e(I)+B(I,J)*f(I);Jacbi((2*I-2),(2*J-2))=B(I,J)*e(I)-G(I,J)*f(I);case2%PV节点Jacbi((2*I-2),(2*J-3))=0;Jacbi((2*I-2),(2*J-2))=0;endelse%对角元素Jacbi(2*I-3,2*J-3)=-(G(I,:)*f+B(I,:)*e)+B(I,I)*e(I)-G(I,I)*f(I);Jacbi(2*I-3,2*J-2)=-(G(I,:)*e-B(I,:)*f)-G(I,I)*e(I)-B(I,I)*f(I);switchVtype(I)case1%PQ节点Jacbi(2*I-2,2*J-3)=-(G(I,:)*e-B(I,:)*f)+G(I,I)*e(I)+B(I,I)*f(I);Jacbi(2*I-2,2*J-2)=(G(I,:)*f+B(I,:)*e)+B(I,I)*e(I)-G(I,I)*f(I);case2%PV节点Jacbi(2*I-2,2*J-3)=-2*f(I);Jacbi(2*I-2,2*J-2)=-2*e(I);endendendendclearIJ;Deta=-inv(Jacbi)*D;%修正方程forx=2:14%新电压初值f(x)=f(x)+Deta((2*x-3),1);e(x)=e(x)+Deta((2*x-2),1);endclearx;U=e+i*f;N=N+1;endN=N-1;S0=U(1)*(conj(Y(1,:))*conj(U));%平衡节点功率S1=zeros(20,1);%始端功率S2=zeros(20,1);%末端功率forx=1:20S1(x)=U(II(x))*(conj(U(II(x)))*(conj(y2(II(x)))+y3(II(x))+i*Qsh(II(x)))+(conj(U(II(x)))-conj(U(JJ(x))))*conj(y1(II(x),JJ(x))));S2(x)=U(JJ(x))*(conj(U(JJ(x)))*(conj(y2(JJ(x)))+y3(JJ(x))+i*Qsh(JJ(x)))+(conj(U(JJ(x)))-conj(U(II(x))))*conj(y1(II(x),JJ(x))));enddetaS=S1+S2;%线路损耗功率Vabs=abs(U);%电压幅值Angle=atan(f./e)*180/pi;%相角%%显示数据disp('迭代次数N=');disp(N);disp('各节点电压');disp('节点幅值相角');disp([BusVabsAngle]);disp('平衡节点功率');disp(S0);disp('线路功率');disp('IIJJ始端功率末端功率线路损耗功率');disp([IIJJS1S2detaS]);Data1及Data2数据文件也已上传文件名就是Data1和Data2