www.pudn.com > powerflow-genetical.gorithm.rar > chaoliu.asv, change:2009-02-25,size:6063b


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=; 
%第一列存储节点号 
%第二列存储节点的功率 
%第三列存储节点的电压 
%第四列存储节点的类型  ;其中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));