diff --git a/src/PowerNetworkMatrices.jl b/src/PowerNetworkMatrices.jl
index 1546cd00..0d6a36e4 100644
--- a/src/PowerNetworkMatrices.jl
+++ b/src/PowerNetworkMatrices.jl
@@ -16,6 +16,7 @@ export to_hdf5
 export validate_connectivity
 export VirtualLODF
 export VirtualPTDF
+export Ward_data
 export Ybus
 
 using DocStringExtensions
@@ -58,5 +59,6 @@ include("lodf_calculations.jl")
 include("virtual_lodf_calculations.jl")
 include("system_utils.jl")
 include("serialization.jl")
+include("ward_calculations.jl")
 
-end
+end
\ No newline at end of file
diff --git a/src/ward_calculations.jl b/src/ward_calculations.jl
new file mode 100644
index 00000000..3997689f
--- /dev/null
+++ b/src/ward_calculations.jl
@@ -0,0 +1,49 @@
+"""
+structure and functions for power systems analysis based on Ward's decomposition
+"""
+
+"""
+Data related to the considered Ward decomposition is here stored.
+Given the set of buses under study, the original Ybus can be divided into 
+sub matrices belongin to study (S) and external (E) buses.
+       |Y_ss  Y_se|
+Ybus = |          |
+       |Y_es  Y_ee|
+
+# Arguments
+- `A::Matrix{ComplexF64}`:
+        the product of Y_se * inv(Y_ee) * Y_es
+- `B::Matrix{ComplexF64}`:
+        the transposed product of Y_se * inv(Y_ee)
+- `study_buses::Vector{Int64}`:
+        vector of the study buses (all other buses will be aggregated according 
+        to the Ward's decomposition)
+"""
+struct Ward_data
+    A::Matrix{ComplexF64}
+    B::Matrix{ComplexF64}
+    study_buses::Vector{Int64}
+end
+
+function Ward_data(
+    ybus::Ybus{Tuple{Vector{Int64}, Vector{Int64}}, Tuple{Dict{Int64, Int64}, Dict{Int64, Int64}}},
+    study_buses::Vector{Int64}
+)
+    # get the indeces of the study and external buses
+    select_x = [ybus.lookup[1][i] for i in study_buses];
+    select_y = [ybus.lookup[2][i] for i in study_buses];
+    unselect_x = [ybus.lookup[1][i] for i in keys(ybus.lookup[1]) if i ∉ study_buses];
+    unselect_y = [ybus.lookup[2][i] for i in keys(ybus.lookup[1]) if i ∉ study_buses];
+    # first reorder the Ybus to get the intermediate sub matrices
+    Y_se = ybus.data[select_x, unselect_y];
+    Y_es = ybus.data[unselect_x, select_y];
+    Y_ee = ybus.data[unselect_x, unselect_y];
+    # get the matrices
+    A = -Y_se*KLU.solve!(KLU.klu(Y_ee), Matrix(Y_es));
+    B = KLU.solve!(KLU.klu(Y_ee), Matrix(transpose(Y_se)));
+    return Ward_data(
+        A,
+        B,
+        study_buses 
+    )
+end
\ No newline at end of file
diff --git a/test/try_YBus_inversion.jl b/test/try_YBus_inversion.jl
new file mode 100644
index 00000000..be951e69
--- /dev/null
+++ b/test/try_YBus_inversion.jl
@@ -0,0 +1,66 @@
+using Revise
+using Test
+import Logging
+import LinearAlgebra: I
+using PowerNetworkMatrices
+using PowerSystems
+using InfrastructureSystems
+using PowerSystemCaseBuilder
+using TimeSeries
+using KLU
+using SparseArrays
+using PowerFlows
+
+const IS = InfrastructureSystems
+const PSY = PowerSystems
+const PSB = PowerSystemCaseBuilder
+const PNM = PowerNetworkMatrices
+
+# ============================================================================
+# ============================================================================
+
+# try the code with a RTS system
+sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys");
+
+# get network matrices
+a_mat = PNM.IncidenceMatrix(sys)
+ba_mat = PNM.BA_Matrix(sys)
+ptdf_mat = PNM.PTDF(sys)
+
+# get the 2T HVDC data, then the buses connecting
+hvdc_line = PSY.get_component(PSY.TwoTerminalHVDCLine, sys, "DC1");
+study_buses = [
+    PSY.get_from(PSY.get_arc(hvdc_line)),
+    PSY.get_to(PSY.get_arc(hvdc_line))
+]
+
+# get all other buses
+other_buses = PSY.get_components(x -> !(x in study_buses), PSY.ACBus, sys)
+
+# get flows and angles
+data = PowerFlowData(DCPowerFlow(), sys)
+power_injection =
+    deepcopy(data.bus_activepower_injection - data.bus_activepower_withdrawals)
+matrix_data = deepcopy(data.power_network_matrix.K)       # LU factorization of ABA
+aux_network_matrix = deepcopy(data.aux_network_matrix)    # BA matrix
+
+valid_ix = setdiff(1:length(power_injection), data.aux_network_matrix.ref_bus_positions)
+ref_bus_angles = deepcopy(data.bus_angles)
+ref_flow_values = deepcopy(data.branch_flow_values)
+
+ref_bus_angles[valid_ix] = matrix_data \ power_injection[valid_ix]
+ref_flow_values = transpose(aux_network_matrix.data) * ref_bus_angles
+
+# get the angles from the PSY data
+thetas = zeros((length(ba_mat.axes[1]),))
+for b in PSY.get_components(PSY.ACBus, sys)
+    idx = ba_mat.lookup[1][PSY.get_number(b)]
+    thetas[idx] = PSY.get_angle(b)
+end
+
+# ! angles are different after computation
+
+
+# now evaluate the Ward decomposition
+ybus = PNM.Ybus(sys)
+ward_data = Ward_data(ybus, PSY.get_number.(study_buses))