MATLAB牛拉法,MATLAB潮流程序(IEEE14直角坐标牛拉法).doc

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

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

clear

baseMVA=100; %功率基值

%%读Data1中数据

load Data1.txt

Bus=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.txt

II=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: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));

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

end

clear x

Y=zeros(14);

for x=1:14

Y(x,x)=sum(y1(x,:))+sum(y2(x,:))+sum(y3(x,:))+i*Qsh(x);

end

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;

end

clear 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)=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(