% fir fin r x b k;
branch=...
[1, 3, 0.082 , 0.427, 0.028, 1;
1, 5, 0.163, 0.754, 0.014, 1;
3, 4, 0.145, 0.581, 0.021, 1;
4, 5, 0.104, 0.518, 0.018, 1;
2, 5, 0.031, 0.248, 0 , 0.95];
% no type p q v a;
bus=...
[1, 0, 1, 1, 1.05, 0;
2, 1, 0.35, 0.35, 1.0 , 0;
3, 2, 0.22, 0.14, 1.0, 0;
4, 2, 0.18, 0.09, 1.0, 0;
5, 2, 0.27, 0.13, 1.0, 0];
Yc=[0,0,0.04*i,0,0];
baseMVA=100; %功率基值%%读Data1中数据
NO=bus(:,1);
Type=bus(:,2); %节点类型
P=bus(:,3); %有功
Q=bus(:,4); %无功
U=bus(:,5); %电压设定点
a=bus(:,6); %并联电容电纳标幺值
II=branch(:,1);
JJ=branch(:,2); %支路端点号
R=branch(:,3); %两点间电阻
X=branch(:,4); %两点间电抗
BD=branch(:,5); %线路对地电纳
K=branch(:,6); %变压器非标准变压比
y1=zeros(5);
y2=zeros(5);
y3=zeros(5);
lin=length(II); %支路数
for x=1:1:lin
if(K(x)==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*BD(x);
y3(JJ(x),II(x))=i*BD(x);
else
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(5);
for x=1:1:5
Y(x,x)=sum(y1(x,:))+sum(y2(x,:))+sum(y3(x,:))+Yc(x);
end
clear x;
Y=Y-y1;
B=imag(Y);
G=real(Y);
ph=find(Type(:,1)==0);
BB=B;
BB(:,ph)=[];
BB(ph,:)=[];
B1=BB;
phpv=find(Type(:,1)<2);
BB=B;
BB(:,phpv)=[];
BB(phpv,:)=[];
B2=BB;
count=0;
E=0.00001;
DPI=zeros(4,1);
DQI=zeros(3,1);
nph=find(Type(:,1)>0);
pq=find(Type(:,1)==2);
PD1=1;
PD2=1;
while((PD1>E)||(PD2>E))
for m=1:1:4
sum1=0;
for n=1:1:5
sum1=sum1+U(nph(m),1)*U(n,1)*(G(nph(m),n)*cos(a(nph(m),1)-a(n,1))+...
B(nph(m),n)*sin(a(nph(m),1)-a(n,1)));
end
DPI(m,1)=P(nph(m),1)-sum1;
end
DAA=(-inv(B1)*(DPI(1:4)./U(2:5)))./U(2:5);
a(2:5)=a(2:5)+DAA(1:4);
PD1=max(abs(DPI(1:4)./U(2:5)));
for m=1:1:3
sum2=0;
for n=1:1:5
sum2=sum2+U(pq(m),1)*U(n,1)*(G(pq(m),n)*sin(abs(a(pq(m),1)-a(n,1)))-...
B(pq(m),n)*cos(abs(a(pq(m),1)-a(n,1))));
end
DQI(m,1)=Q(pq(m),1)-sum2;
end
Uq=U;
Uq(phpv,:)=[];
Upq=Uq;
DUI=(-inv(B2)*(DQI(1:3)./Upq(1:3)));
for m=1:1:3
U(pq(m),1)=U(pq(m),1)+DUI(m,1);
end
PD2=max(abs(DQI(1:3)./Upq(1:3)));
count=count+1;
end
S0=U(1)*(conj(Y(1,:))*conj(U)); %平衡节点功率
S1=zeros(5,1); %始端功率
S2=zeros(5,1); %末端功率
for x=1:1:5
S1(x)=U(II(x))*(conj(U(II(x)))*(conj(y2(II(x)))+y3(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)))+...
(conj(U(JJ(x)))-conj(U(II(x))))*conj(y1(II(x),JJ(x))));
end
DS=S1+S2; %线路损耗功率
Uabs=abs(U); %电压幅值
Angle=a(:,1)*(180/pi); %相角
%%显示数据
disp('迭代次数 conut=');
disp(count);
disp('各节点电压');
disp(' 节点 幅值 相角');
disp([NO Uabs Angle]);
disp('平衡节点功率');
disp(S0);
disp('线路功率');
disp(' II JJ 始端功率 末端功率 线路损耗功率');
disp([II JJ S1 S2 DS]);
aa=1;
