Power Flow in Bipolar DC Distribution Networks

Power Flow in Bipolar DC Distribution Networks Considering Current Limits 10.1109/TPWRS.2022.3181851

%Created for the paper titled 
% "Power Flow in Bipolar DC Distribution Networks Considering Current Limits" 
% This paper will be published in IEEE Transactions on Power Systems
%% POWER FLOW FOR BIPOLAR DC GRIDS
clc
clear
% Take the data
load('Feeder.mat')
% Modifications: switches, voltage regulators and system transformers are eliminated

NumN = Feeder.NumN;NumL = Feeder.NumL;
NODES = zeros(NumN,6);
LINES = zeros(NumL,3);

for k = 1:NumL
    N1 = Feeder.Topology(k,1);
    N2 = Feeder.Topology(k,2);
    Rkm = real(sum(diag(Feeder.ZLIN(:,:,1)))/3)*Feeder.Topology(k,3);
    LINES(k,:) = [N1 N2 Rkm];
end

rand('twister',12345)
Loadpropose1 = 0.8 + (0.4)*rand(NumN,3);
Loadpropose2 = randi([0,1],NumN,3);
LOAD = Loadpropose1.*Loadpropose2; 

LOAD = 2*3490*LOAD/(sum(sum(LOAD)));
Imax = 1.1*LOAD/4.16;
Nodes_ID = unique(LINES(:,[1 2]));
NODES= [Nodes_ID LOAD Imax];


% SYSTEM DATA
Solid = 0;    % 0 for solidly grounded, 1 non-grounded
Vary = 1;     % 0 for Varying inital voltage, 1 for Varying inital load
Vb = 4.16;    % kV rated Voltage
Pb = 1000;    % kW rated Power
Rb = (Vb^2)/(Pb/1000); % only part resistive
Ibase = Pb/Vb;
Vref = 1;
Damping = [20 30];
Nodog = [80 120]; % Controlled Node 
tic
NE = 200;
% Per-Unit transformation
    LINES(:,3) = LINES(:,3)/Rb;
    NODES(:,2:4) = NODES(:,2:4)/Pb;
    NODES(:,5:end) = NODES(:,5:end)/Ibase;
for kk = 1:NE
    if Vary == 1
        NODESin(:,2:4) = NODES(:,2:4)*(1.25*rand+0.75);
    else
        NODESin(:,2:4) = NODES(:,2:4);
    end    
    %% INCIDENCE MATRIX
    N = size(NODES,1);
    L = size(LINES,1);
    A = zeros(N,L);
    for k = 1:L
        n = LINES(k,1); m = LINES(k,2);
        A(n,k) = 1; A(m,k) = -1;
    end
    A3p = zeros(3*N,3*L);
    for m = 1:L
        for k = 1:N
            if A(k,m) == 1
                A3p(3*k-2:3*k,3*m-2:3*m) = eye(3);
            elseif  A(k,m) == -1
                A3p(3*k-2:3*k,3*m-2:3*m) = -eye(3);
            end
        end
    end
    %% PRIMITIVE ADMITANCE MATRIX
    G = diag(1./LINES(:,3));
    G3p = zeros(3*L,3*L);
    for k = 1:L
        G3p(3*k-2:3*k,3*k-2:3*k) = G(k,k)*eye(3);
    end
    G3bus = A3p*G3p*A3p.';
    %% MATRIZ SEPARATION
    G3bds = G3bus(4:end,1:3);
    G3bdd = G3bus(4:end,4:end);
    Z3bdd = inv(G3bdd);
    Vdpn = ones(3*N-3,1);
    Vdpn0 = ones(3*N-3,1);
    PdpnY = ones(2*N-2,1);
    PdpnD = ones(2*N-2,1);
    Vspn = [1;0;-1];
    for k = 1:N-1
        Vdpn0(3*k-2:3*k,1) = Vspn;
        if Vary == 0
            Vdpn0(3*k-2:3*k,1) = (1.1*rand+0.2)*Vspn;
        end
        PdpnY(2*k-1:2*k,1) = [NODESin(k+1,2);NODESin(k+1,3)];
        PdpnD(2*k-1:2*k,1) = [NODESin(k+1,4);NODESin(k+1,4)];
    end
    %% Checking convergency
    z = zeros(3*(N-1),1); k2 = 1;
    for k = 2:N
        for k1 = 1:3
            if NODES(k,k1+1) == 0
            else
                z(k2) = NODES(k,k1+4)^2/abs(NODESin(k,k1+1));
            end
            k2 = k2 +1;
        end
    end
    L = norm(z)*norm(Z3bdd);
    Lkk(kk) = L;
    %% SUCESSIVE APPROXIMATION POWER FLOW METHOD
    itermax = 1000;
    e = 1e-10;
    IdpnY = zeros(size(Vdpn0));
    IdpnD = zeros(size(Vdpn0));
    
    for t = 1:itermax
        ll = 0;
        for k = 2:N
            if any(k == Nodog)
                ll = 1+ ll;
                IdpnY((3*k-5),1) = (PdpnY(2*k-3)-Damping(ll)*(Vref-Vdpn0(3*k-5,1)))/(Vdpn0(3*k-5,1)-Vdpn0(3*k-4,1));
                IdpnY((3*k-4),1) = PdpnY(2*k-3)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-5,1)) + PdpnY(2*k-2)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-3,1));
                IdpnY((3*k-3),1) = (PdpnY(2*k-2)+Damping(ll)*(-Vref-Vdpn0(3*k-3,1)))/(Vdpn0(3*k-3,1)- Vdpn0(3*k-4,1));                
            else
                %% GROUND LOADS
                IdpnY((3*k-5),1) = PdpnY(2*k-3)/(Vdpn0(3*k-5,1)-Vdpn0(3*k-4,1));
                IdpnY((3*k-4),1) = PdpnY(2*k-3)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-5,1)) + PdpnY(2*k-2)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-3,1));
                IdpnY((3*k-3),1) = PdpnY(2*k-2)/(Vdpn0(3*k-3,1)- Vdpn0(3*k-4,1));
            end
            %% LINE-TO-LINE LOADS
            IdpnD((3*k-5),1) = PdpnD(2*k-3)/(Vdpn0(3*k-5,1) - Vdpn0(3*k-3,1));
            IdpnD((3*k-4),1) = 0;
            IdpnD((3*k-3),1) = PdpnD(2*k-2)/(Vdpn0(3*k-3,1) - Vdpn0(3*k-5,1));
            
        end
        Idpn = IdpnY + IdpnD;
        Vdpn = -Z3bdd*(G3bds*Vspn + Idpn);
        if Solid == 0
            for m = 1:N-1
                Vdpn(3*m-1,1) = 0;
            end
        end
        error(t,kk) = max(abs(abs(Vdpn)-abs(Vdpn0)));
        if error(t,kk) < e
            Vpn = [Vspn;Vdpn];
            break;
        else
            Vdpn0 = Vdpn;
        end
    end
end
t1 = toc;
Vall = zeros(NumN,4);
for k = 2:NumN+1    
    Vall(k-1,:) = [k-1, Vpn(3*k-5), Vpn(3*k-4), Vpn(3*k-3)] ;
end