add Yft, Ytf to power_network_matrix Ybus

src/ybus_calculations.jl

@@ -2,55 +2,85 @@
 Nodal admittance matrix (Ybus) is an N x N matrix describing a power system with N buses. It represents the nodal admittance of the buses in a power system.
 
 The Ybus Struct is indexed using the Bus Numbers, no need for them to be sequential
+
+The fields yft and ytf are the branch admittance matrices for the from-to and to-from branch admittances respectively. They are indexed with the branch numbers and the bus numbers.
+The branch numbers are sequential. The bus numbers are represented by fb, tb arrays of sequential bus indices (internal bus indices). 
+Using yft, ytf, and the voltage vector, the branch currents and power flows can be calculated.
 """
 struct Ybus{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF64}
     data::SparseArrays.SparseMatrixCSC{ComplexF64, Int}
     axes::Ax
     lookup::L
+    yft::Union{SparseArrays.SparseMatrixCSC{ComplexF64, Int}, Nothing}
+    ytf::Union{SparseArrays.SparseMatrixCSC{ComplexF64, Int}, Nothing}
+    fb::Union{Vector{Int64}, Nothing}
+    tb::Union{Vector{Int64}, Nothing}
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    y11::Vector{ComplexF64},
+    y12::Vector{ComplexF64},
+    y21::Vector{ComplexF64},
+    y22::Vector{ComplexF64},
     br::PSY.ACBranch,
     num_bus::Dict{Int, Int},
+    branch_ix::Int64,
+    fb::Vector{Int64},
+    tb::Vector{Int64},
 )
     arc = PSY.get_arc(br)
     bus_from_no = num_bus[arc.from.number]
     bus_to_no = num_bus[arc.to.number]
+    fb[branch_ix] = bus_from_no
+    tb[branch_ix] = bus_to_no
     Y_l = (1 / (PSY.get_r(br) + PSY.get_x(br) * 1im))
     Y11 = Y_l + (1im * PSY.get_b(br).from)
     if !isfinite(Y11) || !isfinite(Y_l)
         error(
             "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br)), x = $(PSY.get_x(br))",
         )
     end
-    ybus[bus_from_no, bus_from_no] += Y11
+    y11[branch_ix] = Y11
     Y12 = -Y_l
-    ybus[bus_from_no, bus_to_no] += Y12
-    #Y21 = Y12
-    ybus[bus_to_no, bus_from_no] += Y12
+    y12[branch_ix] = Y12
+    Y21 = Y12
+    y21[branch_ix] = Y21
     Y22 = Y_l + (1im * PSY.get_b(br).to)
-    ybus[bus_to_no, bus_to_no] += Y22
+    y22[branch_ix] = Y22
     return
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    y11::Vector{ComplexF64},
+    y12::Vector{ComplexF64},
+    y21::Vector{ComplexF64},
+    y22::Vector{ComplexF64},
     br::PSY.DynamicBranch,
     num_bus::Dict{Int, Int},
+    branch_ix::Int64,
+    fb::Vector{Int64},
+    tb::Vector{Int64},
 )
-    _ybus!(ybus, br.branch, num_bus)
+    _ybus!(y11, y12, y21, y22, br.branch, num_bus, branch_ix, fb, tb)
     return
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    y11::Vector{ComplexF64},
+    y12::Vector{ComplexF64},
+    y21::Vector{ComplexF64},
+    y22::Vector{ComplexF64},
     br::PSY.Transformer2W,
     num_bus::Dict{Int, Int},
+    branch_ix::Int64,
+    fb::Vector{Int64},
+    tb::Vector{Int64},
 )
     arc = PSY.get_arc(br)
     bus_from_no = num_bus[arc.from.number]
     bus_to_no = num_bus[arc.to.number]
+    fb[branch_ix] = bus_from_no
+    tb[branch_ix] = bus_to_no
     Y_t = 1 / (PSY.get_r(br) + PSY.get_x(br) * 1im)
     Y11 = Y_t
     b = PSY.get_primary_shunt(br)
@@ -60,50 +90,70 @@
         )
     end
 
-    ybus[bus_from_no, bus_from_no] += Y11 - (1im * b)
-    ybus[bus_from_no, bus_to_no] += -Y_t
-    ybus[bus_to_no, bus_from_no] += -Y_t
-    ybus[bus_to_no, bus_to_no] += Y_t
+    y11[branch_ix] = Y11 - (1im * b)
+    Y12 = -Y_t
+    y12[branch_ix] = Y12
+    Y21 = Y12
+    y21[branch_ix] = Y21
+    Y22 = Y_t
+    y22[branch_ix] = Y22
     return
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    y11::Vector{ComplexF64},
+    y12::Vector{ComplexF64},
+    y21::Vector{ComplexF64},
+    y22::Vector{ComplexF64},
     br::PSY.TapTransformer,
     num_bus::Dict{Int, Int},
+    branch_ix::Int64,
+    fb::Vector{Int64},
+    tb::Vector{Int64},
 )
     arc = PSY.get_arc(br)
     bus_from_no = num_bus[arc.from.number]
     bus_to_no = num_bus[arc.to.number]
+    fb[branch_ix] = bus_from_no
+    tb[branch_ix] = bus_to_no
 
     Y_t = 1 / (PSY.get_r(br) + PSY.get_x(br) * 1im)
     c = 1 / PSY.get_tap(br)
     b = PSY.get_primary_shunt(br)
 
     Y11 = (Y_t * c^2)
-    ybus[bus_from_no, bus_from_no] += Y11 - (1im * b)
+    y11[branch_ix] = Y11
     Y12 = (-Y_t * c)
     if !isfinite(Y11) || !isfinite(Y12) || !isfinite(b)
         error(
             "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br)), x = $(PSY.get_x(br))",
         )
     end
-    ybus[bus_from_no, bus_to_no] += Y12
-    #Y21 = Y12
-    ybus[bus_to_no, bus_from_no] += Y12
+
+    y12[branch_ix] = Y12
+    Y21 = Y12
+    y21[branch_ix] = Y21
     Y22 = Y_t
-    ybus[bus_to_no, bus_to_no] += Y22
+    y22[branch_ix] = Y22
     return
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    y11::Vector{ComplexF64},
+    y12::Vector{ComplexF64},
+    y21::Vector{ComplexF64},
+    y22::Vector{ComplexF64},
     br::PSY.PhaseShiftingTransformer,
     num_bus::Dict{Int, Int},
+    branch_ix::Int64,
+    fb::Vector{Int64},
+    tb::Vector{Int64},
 )
     arc = PSY.get_arc(br)
     bus_from_no = num_bus[arc.from.number]
     bus_to_no = num_bus[arc.to.number]
+    fb[branch_ix] = bus_from_no
+    tb[branch_ix] = bus_to_no
     Y_t = 1 / (PSY.get_r(br) + PSY.get_x(br) * 1im)
     tap = (PSY.get_tap(br) * exp(PSY.get_α(br) * 1im))
     c_tap = (PSY.get_tap(br) * exp(-1 * PSY.get_α(br) * 1im))
@@ -114,29 +164,33 @@
             "Data in $(PSY.get_name(br)) is incorrect. r = $(PSY.get_r(br)), x = $(PSY.get_x(br))",
         )
     end
-    ybus[bus_from_no, bus_from_no] += Y11 - (1im * b)
+
+    y11[branch_ix] = Y11 - (1im * b)
     Y12 = (-Y_t / c_tap)
-    ybus[bus_from_no, bus_to_no] += Y12
+    y12[branch_ix] = Y12
     Y21 = (-Y_t / tap)
-    ybus[bus_to_no, bus_from_no] += Y21
+    y21[branch_ix] = Y21
     Y22 = Y_t
-    ybus[bus_to_no, bus_to_no] += Y22
+    y22[branch_ix] = Y22
     return
 end
 
 function _ybus!(
-    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
+    ysh::Vector{ComplexF64},
     fa::PSY.FixedAdmittance,
     num_bus::Dict{Int, Int},
+    fa_ix::Int64,
+    sb::Vector{Int64},
 )
     bus = PSY.get_bus(fa)
     bus_no = num_bus[PSY.get_number(bus)]
+    sb[fa_ix] = bus_no
     if !isfinite(fa.Y)
         error(
             "Data in $(PSY.get_name(fa)) is incorrect. Y = $(fa.Y)",
         )
     end
-    ybus[bus_no, bus_no] += fa.Y
+    ysh[fa_ix] = fa.Y
     return
 end
 
@@ -145,23 +199,43 @@
     buses::Vector{PSY.ACBus},
     fixed_admittances::Vector{PSY.FixedAdmittance},
 )
-    buscount = length(buses)
     num_bus = Dict{Int, Int}()
 
+    branchcount = length(branches)
+    fa_count = length(fixed_admittances)
+    fb = zeros(Int64, branchcount)
+    tb = zeros(Int64, branchcount)
+    sb = zeros(Int64, fa_count)
+
     for (ix, b) in enumerate(buses)
         num_bus[PSY.get_number(b)] = ix
     end
-    ybus = SparseArrays.spzeros(ComplexF64, buscount, buscount)
-    for b in branches
+
+    y11 = zeros(ComplexF64, branchcount)
+    y12 = zeros(ComplexF64, branchcount)
+    y21 = zeros(ComplexF64, branchcount)
+    y22 = zeros(ComplexF64, branchcount)
+    ysh = zeros(ComplexF64, fa_count)
+
+    for (ix, b) in enumerate(branches)
         if PSY.get_name(b) == "init"
             throw(DataFormatError("The data in Branch is invalid"))
         end
-        PSY.get_available(b) && _ybus!(ybus, b, num_bus)
+        PSY.get_available(b) && _ybus!(y11, y12, y21, y22, b, num_bus, ix, fb, tb)
     end
-    for fa in fixed_admittances
-        PSY.get_available(fa) && _ybus!(ybus, fa, num_bus)
+    for (ix, fa) in enumerate(fixed_admittances)
+        PSY.get_available(fa) && _ybus!(ysh, fa, num_bus, ix, sb)
     end
-    return SparseArrays.dropzeros!(ybus)
+    return (
+        y11,
+        y12,
+        y21,
+        y22,
+        ysh,
+        fb,
+        tb,
+        sb,
+    )
 end
 
 """
@@ -175,17 +249,48 @@
     buses::Vector{PSY.ACBus},
     fixed_admittances::Vector{PSY.FixedAdmittance} = Vector{PSY.FixedAdmittance}();
     check_connectivity::Bool = true,
+    make_branch_admittance_matrices::Bool = false,
 )
     bus_ax = PSY.get_number.(buses)
     axes = (bus_ax, bus_ax)
     bus_lookup = make_ax_ref(bus_ax)
+    busnumber = length(buses)
     look_up = (bus_lookup, bus_lookup)
-    ybus = _buildybus(branches, buses, fixed_admittances)
+    y11, y12, y21, y22, ysh, fb, tb, sb = _buildybus(branches, buses, fixed_admittances)
+    ybus = SparseArrays.sparse(
+        [fb; fb; tb; tb; sb],  # row indices
+        [fb; tb; fb; tb; sb],  # column indices
+        [y11; y12; y21; y22; ysh],  # values
+        busnumber,  # size (rows) - setting this explicitly is necessary for the case there are no branches 
+        busnumber,  # size (columns) - setting this explicitly is necessary for the case there are no branches
+    )
+    SparseArrays.dropzeros!(ybus)
     if check_connectivity && length(buses) > 1
         islands = find_subnetworks(ybus, bus_ax)
         length(islands) > 1 && throw(IS.DataFormatError("Network not connected"))
     end
-    return Ybus(ybus, axes, look_up)
+    if make_branch_admittance_matrices
+        yft = SparseArrays.sparse(
+            [1:length(fb); 1:length(fb)],
+            [fb; tb],
+            [y11; y12],
+            length(fb),
+            length(buses),
+        )
+        ytf = SparseArrays.sparse(
+            [1:length(tb); 1:length(tb)],
+            [tb; fb],
+            [y22; y21],
+            length(tb),
+            length(buses),
+        )
+    else
+        yft = nothing
+        ytf = nothing
+        fb = nothing
+        tb = nothing
+    end
+    return Ybus(ybus, axes, look_up, yft, ytf, fb, tb)
 end
 
 """

test/test_ybus.jl

@@ -28,8 +28,10 @@
     @test isapprox(Y5NS[10, 4], -3.3336667 + 33.336667im, atol = 1e-4)
 
     Y5NS = Ybus([branches_5[b] for b in Br5NS_ids], [buses_5[b] for b in Bu5NS_ids])
-    for buf in Bu5NS_ids, but in Bu5NS_ids
-        @test isapprox(Y5NS[buf, but], Ybus5_matpower[buf, but], atol = 1e-3)
+    for buf in Bu5NS_ids
+        for but in Bu5NS_ids
+            @test isapprox(Y5NS[buf, but], Ybus5_matpower[buf, but], atol = 1e-3)
+        end
     end
 
     @test Ybus5[buses_5[1], buses_5[2]] == (-3.5234840209999647 + 35.234840209999646im)