1、MATLAB 潮流程序（ IEEE14 直角坐标 牛拉法）clearbaseMVA=100; %功率基值%读 Data1中数据 load Data1.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中数据load Data2.txtII=Da
2、ta2(:,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); %支路数for x=1:linswitch Ltype(x)case 1y1(II(x),JJ(x)=1/(R(x)+i*X(x);y1(JJ(x),II(x)=y1(II(x),JJ(x);
3、y3(II(x),JJ(x)=i*B(x);y3(JJ(x),II(x)=i*B(x);case 2y1(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);end endclear xY=zeros(14);for x=1:14Y(x,x)=sum(y1(x,:)+sum(y2(x,:)+sum(y3(x,:)+i*Qsh(x);endclear x;Y
4、=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);for x=1:14Ps(x)=(Pgen(x)-Pload(x)/baseMVA;Qs(x)=(Qgen(x)-Qload(x)/baseMVA;endclear x;N=0;Jacbi=zeros(26);while max(abs(D)0.000001for x=2:14 %节点功率及电压不平衡量switch Vtype(x)case 1 %PQ节点D(2*x-3)=Ps(x)-e(x
5、)*(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); case 2 %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); endendclear mfor I=2:14 %求雅克比矩阵for J=2:14if I=J %非对角元素Jacbi(2*I-3)
6、,(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);switch Vtype(I)case 1 %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);case 2 %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(
7、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);switch Vtype(I)case 1 %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);case 2 %PV节点Jacbi(2*I-2,2*J-3)=-2*f(I);Jacbi(2*I-2
8、,2*J-2)=-2*e(I);endendendendclear I J;Deta=-inv(Jacbi)*D; %修正方程for x=2:14 %新电压初值f(x)=f(x)+Deta(2*x-3),1);e(x)=e(x)+Deta(2*x-2),1);endclear x;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); %末端功率for x=1:20S1(x)=U(II(x)*(conj(U(II(x)*(conj(y2(II(x)+y3(
9、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(Bus Vabs Angle);disp(平衡节点功率);disp(S0);disp(线路功率);disp( II JJ 始端功率 末端功率 线路损耗功率);disp(II JJ S1 S2 detaS); Data1及Data2数据文件也已上传 文件名就是 Data1和Data2
