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

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=Data2(

2、:,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:lin switch Ltype(x) case 1 y1(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 2 y1(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:14 Y(x,x)=sum(y1(x,:)+sum(y2(x,:)+sum(y3(x,:)+i*Qsh(x);end

4、clear x;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);for x=1:14 Ps(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.000001 for x=2:14 %节点功率及电压不平衡量 switch Vtype(x) case 1 %PQ节点

5、 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); 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); end end clear m for I=2:14 %求雅克比矩阵 for J=2:14

6、 if I=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); 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; end else

7、%对角元素 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); 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

8、 %PV节点 Jacbi(2*I-2,2*J-3)=-2*f(I); Jacbi(2*I-2,2*J-2)=-2*e(I); end end end endclear 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

9、x=1:20 S1(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(Bus Vabs Angle);disp(&#