i have this code that formd Z bus matrix for power systems

Hello, hope you are doing well,

i have this code that form the Z bus matrix its work well but when i tried to add lines that has transformer its gave me the wrong Z bus matrix so can i have any suggestion.

and thank you.

her is the code

nbus=3; %Nombre total de noeuds
nbr=2;  %Nombre total de branches
linedata=[1   2   0.02   0.06   0.030   1
    1   3   0.03   0.09   0.025   1.2];
%%
Z=zeros(nbus,nbus);
C=zeros(nbr+1,nbus);
A=zeros(1,nbus);
V=0;
l5=0;
B4=2;
a11=0;
c47=0;
g1=linedata(nbr,1);
Q=zeros(nbr+1,nbus-1);
%%
for P=1:nbr
    for O=1:2
        if linedata(P,O)==g1
            if (V==1)
                if O/2 ~= 1
                    X=linedata(P,2);
                    
                elseif O/2 == 1
                    X=linedata(P,1);
                end
                
                C(P,g1(1,1))=X;
                if linedata(P,6)==1
                    R=(linedata(P,5)*1i)^-1;
                    Z(:,X)= Z(:,g1(1,1));
                    Z(X,:)= transpose(Z(:,g1(1,1)));
                    Z(X,X)=Z(g1(1,1),g1(1,1)) +linedata(P,3)+linedata(P,4)*1i;
                    Z=Z-(Z(:,linedata(P,1))*transpose(Z(:,linedata(P,1))))/(Z(linedata(P,1),linedata(P,1))+R);
                    Z=Z-(Z(:,linedata(P,2))*transpose(Z(:,linedata(P,2))))/(Z(linedata(P,2),linedata(P,2))+R);
                else
                    linedata(P,3)=2*linedata(P,3);
                    linedata(P,4)=2*linedata(P,4);
                    R1=(1i*((1-linedata(P,6))/linedata(P,6)^2)*linedata(P,5))^-1;
                    R2=(1i*((linedata(P,6)-1)/linedata(P,6))*linedata(P,5))^-1;
                    Z(:,X)= Z(:,g1(1,1));
                    Z(X,:)= transpose(Z(:,g1(1,1)));
                    Z(X,X)=Z(g1(1,1),g1(1,1)) +linedata(P,3)+linedata(P,4)*1i;
                    Z=Z-(Z(:,linedata(P,1))*transpose(Z(:,linedata(P,1))))/(Z(linedata(P,1),linedata(P,1))+R1);
                    Z=Z-(Z(:,linedata(P,2))*transpose(Z(:,linedata(P,2))))/(Z(linedata(P,2),linedata(P,2))+R2);
                    
                end
                linedata(P,3:5)=0;
            else
                %  First Z Matix
                if O/2 ~= 1
                    X=linedata(P,2);
                elseif O/2 == 1
                    X=linedata(P,1);
                end
                C(P,g1)=X;
                C(nbr+1,g1)=g1;
                V=V+1;
                if linedata(P,6)==1
                    Z(linedata(P,1),linedata(P,1))=(linedata(P,5)*1i)^-1;
                    Z(linedata(P,2),linedata(P,2))=(linedata(P,5)*1i)^-1;
                    K= Z(:,linedata(P,1))- Z(:,linedata(P,2));
                    Z=Z-(K*transpose(K))*((Z(linedata(P,1),linedata(P,1))+Z(linedata(P,2),linedata(P,2))-2*Z(linedata(P,2),linedata(P,1))+linedata(P,3)+linedata(P,4)*1i))^-1;
                else
                    linedata(P,3)=2*linedata(P,3);
                    linedata(P,4)=2*linedata(P,4);
                    Z(linedata(P,1),linedata(P,1))=(1i*((1-linedata(P,6))/linedata(P,6)^2)*linedata(P,5))^-1;
                    Z(linedata(P,2),linedata(P,2))=(1i*((linedata(P,6)-1)/linedata(P,6))*linedata(P,5))^-1;
                    K= Z(:,linedata(P,1))- Z(:,linedata(P,2));
                    Z=Z-(K*transpose(K))*((Z(linedata(P,1),linedata(P,1))+Z(linedata(P,2),linedata(P,2))-2*Z(linedata(P,2),linedata(P,1))+linedata(P,3)+linedata(P,4)*1i))^-1;
                end
                linedata(P,3:5)=0;
            end
        end
    end
end
A(1,1)=g1;
for x1=1:nbus
    if x1~=g1
        for P=1:nbr
            for O=1:2
                if linedata(P,O)==x1
                    if O/2 ~= 1
                        X=linedata(P,2);
                    elseif O/2 == 1
                        X=linedata(P,1);
                    end
                    C(P,x1)=X;
                    C(nbr+1,x1)=x1;
                end
                
            end
        end
    end
end
for p=1:nbus
    m=A(1,p);
    for j5=1:nbr
        v7=C(j5,m);
        if v7~=0
            for z2=1:nbus
                if v7==A(1,z2)
                    l5=1;
                    break;
                end
            end
            if l5~=1
                A(1,B4)=v7;
                B4=B4+1;
            end
            if l5==1
                l5=0;
            end
            if B4==6
                break;
            end
        end
    end
end
for h11=2:nbus
    g11=A(1,h11);
    for u11=1:nbr
        for f4=1:2
            if g11== linedata(u11,f4)
                if f4/2 ~= 1
                    X=linedata(u11,2);
                elseif f4/2 == 1
                    X=linedata(u11,1);
                end
                M=C;
                Q(:,A(1,h11-1))=M(:,A(1,h11-1));
                for v47=1:nbus-1
                    for b47=1:nbr+1
                        if X==Q(b47,v47)
                            c47=1;
                            break;
                        end
                    end
                end
                
                if c47~=1
                    if linedata(P,6)==1
                        if linedata(u11,3)+linedata(u11,4)*1i~=0
                            Z(:,X)= Z(:,g11);
                            Z(X,:)= transpose(Z(:,g11));
                            Z(X,X)=Z(g11,g11) +linedata(u11,3)+linedata(u11,4)*1i;
                        end
                        R=(linedata(u11,5)*1i)^-1;
                        if (R~=inf)
                            Z=Z-(Z(:,linedata(u11,1))*transpose(Z(:,linedata(u11,1))))/(Z(linedata(u11,1),linedata(u11,1))+R);
                            Z=Z-(Z(:,linedata(u11,2))*transpose(Z(:,linedata(u11,2))))/(Z(linedata(u11,2),linedata(u11,2))+R);
                        end
                    else
                        linedata(u11,3)=2*linedata(u11,3);
                        linedata(u11,4)=2*linedata(u11,4);
                        R1=(1i*((1-linedata(u11,6))/linedata(u11,6)^2)*linedata(u11,5))^-1;
                        R2=(1i*((linedata(u11,6)-1)/linedata(u11,6))*linedata(u11,5))^-1;
                        if linedata(u11,3)+linedata(u11,4)*1i~=0
                            Z(:,X)= Z(:,g11);
                            Z(X,:)= transpose(Z(:,g11));
                            Z(X,X)=Z(g11,g11) +linedata(u11,3)+linedata(u11,4)*1i;
                        end
                        
                        if (R1~=inf) && (R2~=inf)
                            Z=Z-(Z(:,linedata(u11,1))*transpose(Z(:,linedata(u11,1))))/(Z(linedata(u11,1),linedata(u11,1))+R1);
                            Z=Z-(Z(:,linedata(u11,2))*transpose(Z(:,linedata(u11,2))))/(Z(linedata(u11,2),linedata(u11,2))+R2);
                        end
                    end
                    linedata(u11,3:5)=0;
                end
                if c47==1
                    if linedata(P,6)==1
                        if (linedata(u11,3)+linedata(u11,4)*1i)~=0
                            K=Z(:,X)-Z(:,g11);
                            Z=Z-(K*transpose(K))*(Z(g11,g11)+Z(X,X)-2*Z(X,g11)+linedata(u11,3)+linedata(u11,4)*1i)^-1;
                        end
                        R=(linedata(u11,5)*1i)^-1;
                        if (R~=inf)
                            Z=Z-(Z(:,linedata(u11,1))*transpose(Z(:,linedata(u11,1))))/(Z(linedata(u11,1),linedata(u11,1))+R);
                            Z=Z-(Z(:,linedata(u11,2))*transpose(Z(:,linedata(u11,2))))/(Z(linedata(u11,2),linedata(u11,2))+R);
                        end
                    else
                        linedata(u11,3)=2*linedata(u11,3);
                        linedata(u11,4)=2*linedata(u11,4);
                        R1=(1i*((1-linedata(u11,6))/linedata(u11,6)^2)*linedata(u11,5))^-1;
                        R2=(1i*((linedata(u11,6)-1)/linedata(u11,6))*linedata(u11,5))^-1;
                        if (linedata(u11,3)+linedata(u11,4)*1i)~=0
                            K=Z(:,X)-Z(:,g11);
                            Z=Z-(K*transpose(K))*(Z(g11,g11)+Z(X,X)-2*Z(X,g11)+linedata(u11,3)+linedata(u11,4)*1i)^-1;
                        end
                        
                        if (R1~=inf) && (R2~=inf)
                            Z=Z-(Z(:,linedata(u11,1))*transpose(Z(:,linedata(u11,1))))/(Z(linedata(u11,1),linedata(u11,1))+R1);
                            Z=Z-(Z(:,linedata(u11,2))*transpose(Z(:,linedata(u11,2))))/(Z(linedata(u11,2),linedata(u11,2))+R2);
                        end
                    end
                    linedata(u11,3:5)=0;
                    c47=0;
                end
                
            end
        end
    end
end
disp(Z)

   0.0052 -16.4605i  -0.0047 -16.4901i   0.0011 -16.4728i
  -0.0047 -16.4901i   0.0054 -16.4598i  -0.0089 -16.5025i
   0.0011 -16.4728i  -0.0089 -16.5025i   0.0570 -16.3051i

Answer 1

[Anushka]

Hi @Abdelkrim,

In my understanding, your code aims to construct a Z-bus matrix for a power sytem, including the complexities introduced by transformers. To ensure accurate handling of transforners in your Z-bus matrix code, consider the following steps:

1. Ensure the transformer's impedance and turns ratio are correctly represented in your calculations.
2. Verify the impedance is properly adjusted according to the turns ratio and check the indexing for bus number to avoid errors.
3. Use debug statements to trace matrix updates and validate results agianst known test cases. Insert print statements at key points in code to monitor changes:

disp('Updating Z-bus matrix for branch:');
disp(['From bus: ', num2str(linedata(P,1)), ' To bus: ', num2str(linedata(P,2))]);
disp('Current Z-matrix:');
disp(Z);

4. Compare results with known test cases or use simulation tools to ensure the Z-bus matrix is accurate.

These steps will help you track Z-bus matrix updates and ensure each step is executed as expected.

Hope this helps!