Fixing bug on PTDF calculation with distributed slack

src/common.jl

@@ -285,7 +285,7 @@ function make_entries_zero!(vector::Vector{Float64}, tol::Float64)
             vector[i] = 0.0
         end
     end
-    return vector
+    return
 end
 
 """

src/ptdf_calculations.jl

@@ -140,6 +140,7 @@ function _calculate_PTDF_matrix_KLU(
     # inizialize matrices for evaluation
     valid_ix = setdiff(1:buscount, ref_bus_positions)
     PTDFm_t = zeros(buscount, linecount)
+
     if !isempty(dist_slack) && length(ref_bus_positions) != 1
         error(
             "Distibuted slack is not supported for systems with multiple reference buses.",
@@ -153,13 +154,15 @@ function _calculate_PTDF_matrix_KLU(
         @info "Distributed bus"
         copyto!(PTDFm_t, BA)
         PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :])
-        PTDFm_t[ref_bus_positions, :] .= 0.0
+        PTDFm_t[collect(ref_bus_positions), :] .= 0.0
         slack_array = dist_slack / sum(dist_slack)
         slack_array = reshape(slack_array, 1, buscount)
-        return PTDFm_t - (PTDFm_t * slack_array) * ones(1, buscount)
+        return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)
     else
         error("Distributed bus specification doesn't match the number of buses.")
     end
+
+    return
 end
 
 """
@@ -246,15 +249,21 @@ function _calculate_PTDF_matrix_DENSE(
         return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
+        PTDFm_t[valid_ixs, :] = gemm(
+            'N',
+            'N',
+            getri!(ABA, bipiv),
+            BA[valid_ixs, :],
+        )
         slack_array = dist_slack / sum(dist_slack)
-        slack_array = reshape(slack_array, buscount, 1)
+        slack_array = reshape(slack_array, 1, buscount)
         return PTDFm_t -
-               gemm('N', 'N', gemm('N', 'N', PTDFm, slack_array), ones(1, buscount))
+               gemm('N', 'N', ones(buscount, 1), gemm('N', 'N', slack_array, PTDFm_t))
     else
         error("Distributed bus specification doesn't match the number of buses.")
     end
 
-    return PTDFm_t
+    return
 end
 
 """
@@ -316,28 +325,29 @@ function _calculate_PTDF_matrix_MKLPardiso(
     ps = Pardiso.MKLPardisoSolver()
     Pardiso.set_iparm!(ps, 59, 2)
     Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
-
     Pardiso.solve!(ps, ABA_inv, ABA, Ix)
+    # inizialize matrices for evaluation
+    valid_ix = setdiff(1:buscount, ref_bus_positions)
     PTDFm_t = zeros(buscount, linecount)
-    row_idx = setdiff(1:size(PTDFm_t, 1), ref_bus_positions)
 
     if !isempty(dist_slack) && length(ref_bus_positions) != 1
         error(
             "Distibuted slack is not supported for systems with multiple reference buses.",
         )
-    elseif isempty(dist_slack) && length(ref_bus_positions) < buscount
-        PTDFm_t[row_idx, :] = ABA_inv * @view BA[row_idx, :]
+    elseif isempty(dist_slack) && length(ref_bus_positions) != buscount
+        PTDFm_t[valid_ix, :] = ABA_inv * @view BA[valid_ix, :]
+        return PTDFm_t
     elseif length(dist_slack) == buscount
         @info "Distributed bus"
-        PTDFm_t[row_idx, :] = ABA_inv * @view BA[row_idx, :]
+        PTDFm_t[valid_ix, :] = ABA_inv * @view BA[valid_ix, :]
         slack_array = dist_slack / sum(dist_slack)
-        slack_array = reshape(slack_array, buscount, 1)
-        PTDFm_t = PTDFm_t - (PTDFm_t * slack_array) * ones(1, buscount)
+        slack_array = reshape(slack_array, 1, buscount)
+        return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)
     else
         error("Distributed bus specification doesn't match the number of buses.")
     end
 
-    return PTDFm_t
+    return
 end
 
 """

src/virtual_lodf_calculations.jl

@@ -242,10 +242,9 @@ function _getindex(
         lodf_row[row] = -1.0
         # add slack bus value (zero) and make copy of temp into the cache
         if get_tol(vlodf) > eps()
-            vlodf.cache[row] = make_entries_zero!(deepcopy(lodf_row), get_tol(vlodf))
-        else
-            vlodf.cache[row] = deepcopy(lodf_row)
+            make_entries_zero!(lodf_row, get_tol(vlodf))
         end
+        vlodf.cache[row] = deepcopy(lodf_row)
         return vlodf.cache[row][column]
     end
 end

src/virtual_ptdf_calculations.jl

@@ -71,6 +71,9 @@ function VirtualPTDF(
     tol::Float64 = eps(),
     max_cache_size::Int = MAX_CACHE_SIZE_MiB,
     persistent_lines::Vector{String} = String[])
+    if length(dist_slack) != 0
+        @info "Distributed bus"
+    end
 
     #Get axis names
     line_ax = [PSY.get_name(branch) for branch in branches]
@@ -171,19 +174,35 @@ function _getindex(
         # Needs improvement
         valid_ix = vptdf.valid_ix
         lin_solve = KLU.solve!(vptdf.K, Vector(vptdf.BA[valid_ix, row]))
+        buscount = size(vptdf, 1)
 
-        # ! missing dist_slack case
-
-        for i in eachindex(valid_ix)
-            vptdf.temp_data[valid_ix[i]] = lin_solve[i]
+        if !isempty(vptdf.dist_slack) && length(vptdf.ref_bus_positions) != 1
+            error(
+                "Distibuted slack is not supported for systems with multiple reference buses.",
+            )
+        elseif isempty(vptdf.dist_slack) && length(vptdf.ref_bus_positions) < buscount
+            for i in eachindex(valid_ix)
+                vptdf.temp_data[valid_ix[i]] = lin_solve[i]
+            end
+            vptdf.cache[row] = deepcopy(vptdf.temp_data)
+        elseif length(vptdf.dist_slack) == buscount
+            for i in eachindex(valid_ix)
+                vptdf.temp_data[valid_ix[i]] = lin_solve[i]
+            end
+            slack_array = vptdf.dist_slack / sum(vptdf.dist_slack)
+            slack_array = reshape(slack_array, buscount)
+            vptdf.cache[row] =
+                deepcopy(vptdf.temp_data .- dot(vptdf.temp_data, slack_array))
+        else
+            error("Distributed bus specification doesn't match the number of buses.")
         end
 
         # add slack bus value (zero) and make copy of temp into the cache
         if get_tol(vptdf) > eps()
-            vptdf.cache[row] = make_entries_zero!(deepcopy(vptdf.temp_data), get_tol(vptdf))
-        else
-            vptdf.cache[row] = deepcopy(vptdf.temp_data)
+            # vptdf.cache[row] = make_entries_zero!(deepcopy(vptdf.temp_data), get_tol(vptdf))
+            make_entries_zero!(vptdf.cache[row], get_tol(vptdf))
         end
+
         return vptdf.cache[row][column]
     end
 end

test/test_ptdf.jl

@@ -192,3 +192,20 @@ end
     # Test Throw error when isolated buses are connected to an available branch
     @test_throws IS.ConflictingInputsError ptdf_2 = PTDF(sys_2)
 end
+
+@testset "PTDF matrices with distributed slack" begin
+    sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5")
+
+    buscount = length(PNM.get_buses(sys5))
+
+    dist_slack = 1 / buscount * ones(buscount)
+    slack_array = dist_slack / sum(dist_slack)
+
+    P5_1 = PTDF(sys5; dist_slack = slack_array, linear_solver = "KLU")
+    P5_2 = PTDF(sys5; dist_slack = slack_array, linear_solver = "Dense")
+    P5_3 = PTDF(sys5; dist_slack = slack_array, linear_solver = "MKLPardiso")
+
+    @assert isapprox(P5_1.data, P5_2.data, atol = 1e-5)
+    @assert isapprox(P5_1.data, P5_3.data, atol = 1e-5)
+    @assert isapprox(P5_2.data, P5_3.data, atol = 1e-5)
+end

test/test_virtual_ptdf.jl

@@ -67,3 +67,19 @@ end
         @test vptdf.lookup[1][l] ∈ keys(vptdf.cache.temp_cache)
     end
 end
+
+@testset "Test virtual PTDF with distributed slack" begin
+    # get 5 bus system
+    sys = PSB.build_system(PSB.PSITestSystems, "c_sys5")
+    bus_number = length(PNM.get_buses(sys))
+    dist_slack = 1 / bus_number * ones(bus_number)
+    # compute full PTDF
+    ptdf = PTDF(sys; dist_slack = dist_slack)
+    # compute each row of the virtual PTDF and compare values
+    vptdf = VirtualPTDF(sys; dist_slack = dist_slack)
+    for row in 1:size(ptdf.data, 2)
+        # evaluate the column (just needs one element)
+        vptdf[row, 1]
+        @assert isapprox(vptdf.cache[row], ptdf[row, :], atol = 1e-5)
+    end
+end