Y Bus Assembly (Using Singular Tranformation) for mutually coupled elements

yBusSingular(Mutual).m

clc; clearvars;
% Data from IEEE 5-bus Data
% "http://shodhganga.inflibnet.ac.in/bitstream/10603/26549/14/14_appendix.pdf""
lineData = ...
[
%   line   sending bus   receiving bus     line impedance
     1          1             2              0.02+0.060j
     2          1             3              0.08+0.240j
     3          2             3              0.06+0.250j
     4          2             4              0.06+0.180j
     5          2             5              0.04+0.120j
     6          3             4              0.01+0.030j
     7          4             5              0.08+0.240j
];

mutualCoupling = ...
[
%     line 1  line 2      mutual impedance
        2       3           0.00+0.500j
        5       6           0.00+0.250j
];
numberMutuallyCoupledLines = length(mutualCoupling(:,1));
% Calculate the number of lines in the line data
lines = length(lineData);
% Calculate the number of buses in the line data
buses = max(max(lineData(:,2)),max(lineData(:,3)));
% Calculate the line admittance by taking element by element inverse
lineAdmittance = lineData(:,4).^-1;
% Transform the lineAdmittance array into a diagonal matrix
y = diag(lineAdmittance);
% Fill the primitive y matrix with mutual coupling elements
for x = 1:numberMutuallyCoupledLines
    line1 = mutualCoupling(x,1);    line2 = mutualCoupling(x,2);
    y(line1,line2) = mutualCoupling(x,3);
    y(line2,line1) = mutualCoupling(x,3);
end
% Initialize the admittance matrix with dimensions (lines X buses)
admittanceMatrix = zeros(lines,buses);

% Fill the Admittance Matrix by reading the lineData
for k = 1:lines
    sendingBus = lineData(k,2);
    receivingBus = lineData(k,3);
    % Fill in the sending nodes
    if sendingBus ~= 0               % condition used to remove reference node
        admittanceMatrix(k, sendingBus) = 1;
    end
    % Fill in the receiving nodes
    if receivingBus ~= 0              % condition used to remove reference node
        admittanceMatrix(k, receivingBus) = -1;
    end
end

% Calculate the Y Bus
ybus = admittanceMatrix'*y*admittanceMatrix;

disp('Ybus (Using singular transformation technique)');     disp(ybus);