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

2、=Data2(:,1);JJ=Data2(:,2);% 支路端点号Ltype=Data2(:,5);% 线路类型R=Data2(:,6);% 两点间电阻X=Data2(:,7);% 两点间电抗B=Data2(:,8)/2;% 线路对地电纳K=Data2(:,14);% 变压器非标准变压比.% 求导纳矩阵 Y y1=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)

3、,JJ(x);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);endendclear xY=zeros(14);for x=1:14.Y(x,x)=sum(y1(x,:)+sum(y2(x,:)+sum(y3(x,:)+i*Qsh(x);endcl

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

5、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);endendclear mfor I=2:14% 求雅克比矩阵for J=2:14if I=J% 非对角元素Jacbi(

6、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;endelse% 对角元素Jacbi(2*I-3,2*J-3)=-(G(I,

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

8、i(2*I-2,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

9、(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