www.pudn.com > powerflow-genetical.gorithm.rar > chaoliu.m, change:2009-05-20,size:6176b

```function [U,sunhao]=chaoliu(B2)
%本程序的功能是用牛顿-拉夫逊法进行潮流计算
n=5;
nl=5;
isb=5;
pr=0.00001;
X=[1,0;2,0;3,0;4,0;5,0];
%第一列存储节点号
%第二列存储对地导纳，若无则取0
B1=[1,2,0.04+0.25j,0.50j;1,3,0.1+0.35j,0;2,3,0.08+0.30j,0.50j;2,-4,0.015j,1.05;3,-5,0.03j,1.05];
%第一列存储支路的首端点
% 第二列存储支路的末端点
% 第三列存储支路的阻抗
% 第四列存储支路的对地导纳或变压器支路的变比
%B2=;
B2=[1,-1.6-0.8j,1,1;2,-2-4j,1,1;3,-3.7-1.3j,1,1;4,5,1.05,2;5,0,1.00,3];
%第一列存储节点号
%第二列存储节点的功率
%第三列存储节点的电压
%第四列存储节点的类型  ;其中1为PQ节点，2为PV节点，3为平衡节点
Y=zeros(n);U=zeros(1,n);cta=zeros(1,n);S=zeros(1,n);
for i=1:n
if X(i,2)~=0;
p=X(i,1);
Y(p,p)=X(i,2);
end
end
for i=1:nl
if B1(i,1)>0&B1(i,2)>0
p=B1(i,1);q=B1(i,2);
Y(p,q)=Y(p,q)-1./B1(i,3);
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)/2;
Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)/2;
else p=abs(B1(i,1));q=abs(B1(i,2));
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,4));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+1./(B1(i,3)*B1(i,4)^2);
Y(q,q)=Y(q,q)+1./B1(i,3);
end
end
disp(Y);            %形成节点导纳矩阵
G=real(Y);B=imag(Y);
for i=1:n
cta(i)=angle(B2(i,3));%i节点的相角
U(i)=abs(B2(i,3));%i节点的电压幅值
end
for i=1:n
S(i)=B2(i,2);
end
P=real(S);Q=imag(S);
ICT1=0;IT2=1;
while IT2~=0
IT2=0;t1=1;t2=1;
for i=1:n
if i~=isb%isb已经定义为5，即平衡节点不参与运算，程序这样写和给出的B2的型式是相关的，B2中其平衡节点在最后一列
C(i)=0;
D(i)=0;
for j1=1:n
C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));%Pi
D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));%Qi
end
DP(t1)=P(i)-C(i);%\$Pi
t1=t1+1;
if B2(i,4)==1           %若为pq节点,也就是若为pv节点，就不必算\$Qi了
DQ(t2)=Q(i)-D(i);%\$Qi
t2=t2+1;
end
end
end
t1=t1-1;t2=t2-1;
DPQ=[DP';DQ'];             %求DP,DQ
for i=1:t1+t2
if abs(DPQ(i))>pr%pr已经定义了为0.00001，功率不平衡量要小于这个值，则迭代完成
IT2=IT2+1;
end
end
H=zeros(t1,t1);N=zeros(t1,t2);J=zeros(t2,t1);L=zeros(t2,t2);
for i=1:t1
for j1=1:t1
if j1~=isb&j1~=i
H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
elseif j1~=isb&j1==i
H(i,j1)=U(i)^2*B(i,j1)+D(i);%？
end
end
end
for i=1:t1
for j1=1:t2
if j1~=isb&j1~=i
N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
elseif j1~=isb&j1==i
N(i,j1)=0-U(i)^2*G(i,j1)-C(i);%？
end
end
end
for  i=1:t2
for j1=1:t1
if j1~=isb&j1~=i
J(i,j1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
elseif j1~=isb&j1==i
J(i,j1)=U(i)^2*G(i,j1)-C(i);%？
end
end
end
for i=1:t2
for j1=1:t2
if j1~=isb&j1~=i
L(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
elseif j1~=isb&j1==i
L(i,j1)=U(i)^2*B(i,j1)-D(i);%？
end
end
end
K=[H,N;J,L];%求雅可比矩阵
modify=-K\DPQ;
Dcta=modify([1:t1],:);%相角的偏差量
t3=U(:,[1:t2]);
DU=diag(t3,0)*modify([t1+1:t1+t2],:);%diag是对角阵
t4=1;
for i=1:t1
if B2(i,4)~=3
cta(1,i)=cta(1,i)+Dcta(t4,1);%相角修正
t4=t4+1;
end
end
t5=1;
for i=1:t2
if B2(i,4)==1
U(1,i)=U(1,i)+DU(t5,1);%电压修正
t5=t5+1;
end
end
ICT1=ICT1+1;
if ICT1>15
%disp('不收敛');
sunhao=200;

break;
end
end       %修正原值
for i=1:n
UU(i)=U(i)*cos(cta(i))+1i*U(i)*sin(cta(i));
end
for p=1:n
c(p)=0;
for q=1:n
c(p)=c(p)+conj(Y(p,q))*conj(UU(q));
end
s(p)=UU(p)*c(p);%
end
%disp('--------------------------------------------------------------------------------');
%disp('各节点电压U为（节点从小到大排列）:');
%disp(UU);
%disp('--------------------------------------------------------------------------------');
%disp('各节点电压幅值为（节点从小到大排列）:');
%disp(abs(UU));
%disp('--------------------------------------------------------------------------------');
%disp('各节点电压相角为（节点从小到大排列）:');
%disp(180*angle(UU)/pi);
%disp('--------------------------------------------------------------------------------');
%disp('按公式计算全部线路功率，结果如下:');
for i=1:nl
if B1(i,1)>0&B1(i,2)>0
p=B1(i,1);q=B1(i,2);
Si(p,q)=UU(p)*(conj(UU(p))*conj(B1(i,4)./2)+(conj(UU(p))-conj(UU(q)))*conj(1./(B1(i,3))));
else p=abs(B1(i,1));q=abs(B1(i,2));
Si(p,q)=UU(p)*((conj(UU(p)*B1(i,4))-conj(UU(q)))*conj(1./(B1(i,3)*B1(i,4))));%各条支路首端功率Si
end
f=[p,q,Si(p,q)];
%disp(f);
end
for i=1:nl
if B1(i,1)>0&B1(i,2)>0
p=B1(i,1);q=B1(i,2);
Sj(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q))-conj(UU(p)))*conj(1./(B1(i,3))));
else p=abs(B1(i,1));q=abs(B1(i,2));
Sj(q,p)=UU(q)*((conj(UU(q)*B1(i,4))-conj(UU(p)))*conj(1./(B1(i,3)*B1(i,4))));%各条支路末端功率Sj
end
f=[q,p,Sj(q,p)];
%disp(f);
end
%disp('--------------------------------------------------------------------------------');
%disp('各条支路的功率损耗DS为(顺序同您输入B1时一样):');
for i=1:nl
if B1(i,1)>0&B1(i,2)>0
p=B1(i,1);q=B1(i,2);
else p=abs(B1(i,1));q=abs(B1(i,2));
end
DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS
%disp(DS(i));
end
Sp=0;
for i=1:n
Sp=Sp+UU(isb)*conj(Y(isb,i))*conj(UU(i));%平衡节点功率
end
%disp('平衡节点的功率:');
%disp(Sp);下面是用于遗传算法的部分程序
q(1)=DS(1);
for i=2:5
q(i)=q(i-1)+DS(i);            %累加个体适应度形成赌轮
end
sunhao=real(q(5));```