loadflow.m

一种前推回代的配电网潮流计算

clc;
clear;
eval nine_node
[ busnum, dump ] = size( Bus );
Vbus = ones(busnum,1);
Vbus(1,1) = 1.02;
Vbus2 = Vbus;
busindex = Bus(:,1);

[ branchnum, dump ] = size( Branch ); 
frombus = Branch(:,1);
tobus = Branch(:,2);

errorr = 1;
kkk = 0;
Pij = zeros( branchnum+1, 1 );
Qij = zeros( branchnum+1, 1 );
Ploss = zeros( branchnum, 1 );
Qloss = zeros( branchnum, 1 );
while errorr > 1e-6;
    [ endbus,branchindex ] = max ( tobus );
    ii = endbus;%子节点 
    kk = Branch( branchindex,1 );%父节点
     
    for loop = 1:1:branchnum;
   
        Pload = Bus( ii,2 );
        Qload = Bus( ii,3 );
        
        R = Branch( branchindex ,3 );
        X = Branch( branchindex, 4 );
        Pij0 = Pij( ii+1 );
        Qij0 = Qij( ii+1 );
        VV = Vbus( ii );
        II = (( Pload + Pij0 )^2 + ( Qload + Qij0 )^2)/(VV*conj(VV));
        Ploss( ii ) = II * R;
        Qloss( ii ) = II * X;
        Pij( kk+1 ) = Pload + Pij0 + Ploss( ii );
        Qij( kk+1 ) = Qload + Qij0 + Qloss( ii );
       
        if kk ~= 0; 
            ii = kk;
            branchindex = find( tobus == ii );
            kk = Branch( branchindex, 1 );
        end
    end
    
    lineindex = find( frombus == kk );
    ii = Branch( branchindex, 2 );
    
    for  loop = 1:1:busnum;
        R = Branch( branchindex ,3 );
        X = Branch( branchindex, 4 );
   
        VV =   Vbus( kk+1 );
        Iki = ( Pij( kk+1 ) - j*Qij( kk+1 ))/conj( VV );
        
        Vbus( ii+1,1 ) = VV - Iki * ( R+j*X );
      
        if ii ~= endbus
           kk = ii;
           branchindex = find( frombus == kk );
           ii = Branch( branchindex, 2 );
        end
    end
    
    errorr = max( abs(Vbus2-Vbus) );
    Vbus2 = Vbus;
    kkk = kkk + 1;
end

Vbusamp = abs(Vbus);
bar( Vbusamp )

Lij = zeros( busnum-1, 1 );
for  ii = 1:(busnum-1);   
    Pj = Pij( ii ) - Ploss( ii );
    Qj = Qij( ii ) - Qloss( ii );
    Rij = Branch( ii,3 );
    Xij = Branch( ii,4 );
    Ui = Vbusamp( ii );
    Lij( ii )=4*( (Pj*Xij-Qj*Rij)^2 + (Pj*Rij+Qj*Xij)*Ui^2 )/Ui^4;    
end

Vbus_V=zeros(branchnum+1,1);
aa=zeros(branchnum,1);
LL=zeros(branchnum,1);
for i=1:branchnum
    Vbus_V(i)=abs(Vbus(i));
    aa(i)=angle(Vbus(i));  
    LL(i)=(2*cos(aa(i))*Vbus_V(i+1)/Vbus_V(i)-1);
end