From 3ea3d1216876725f87fa618a638c3220c72f0a2f Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Mon, 31 Jul 2023 15:17:09 -0600 Subject: [PATCH 01/11] added fix to distributed slack and new test --- src/ptdf_calculations.jl | 14 ++++++++++---- test/test_ptdf.jl | 19 +++++++++++++++++++ 2 files changed, 29 insertions(+), 4 deletions(-) diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl index 83074943..9f1f5312 100644 --- a/src/ptdf_calculations.jl +++ b/src/ptdf_calculations.jl @@ -153,10 +153,10 @@ 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 @@ -246,10 +246,16 @@ 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 diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index ccb24643..e2ed5467 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -192,3 +192,22 @@ end # Test Throw error when isolated buses are connected to an available branch @test_throws IS.ConflictingInputsError ptdf_2 = PTDF(sys_2) end + +@test "PTDF matrices with distributed slack" begin + sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") + + bus_number = length(PNM.get_buses(sys)) + + dist_slack = 1/bus_number*ones(bus_number) + slack_array = dist_slack / sum(dist_slack) + slack_array = reshape(slack_array, 1, buscount) + + P5_1 = PTDF(sys5; slack_array=slack_array, linear_solver = "KLU") + P5_2 = PTDF(sys5; slack_array=slack_array, linear_solver = "Dense") + P5_3 = PTDF(sys5; slack_array=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 \ No newline at end of file From 11679e602602924671399902f08500266a7f50af Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Mon, 31 Jul 2023 15:34:50 -0600 Subject: [PATCH 02/11] found error --- test/test_ptdf.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index e2ed5467..51e5c5a5 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -193,7 +193,7 @@ end @test_throws IS.ConflictingInputsError ptdf_2 = PTDF(sys_2) end -@test "PTDF matrices with distributed slack" begin +@testset "PTDF matrices with distributed slack" begin sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") bus_number = length(PNM.get_buses(sys)) From ebc8b445bea081d80b10154a8e6b7d812b5791be Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Mon, 31 Jul 2023 15:38:45 -0600 Subject: [PATCH 03/11] formatter --- test/test_ptdf.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index 51e5c5a5..aeacd386 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -198,16 +198,15 @@ end bus_number = length(PNM.get_buses(sys)) - dist_slack = 1/bus_number*ones(bus_number) + dist_slack = 1 / bus_number * ones(bus_number) slack_array = dist_slack / sum(dist_slack) slack_array = reshape(slack_array, 1, buscount) - P5_1 = PTDF(sys5; slack_array=slack_array, linear_solver = "KLU") - P5_2 = PTDF(sys5; slack_array=slack_array, linear_solver = "Dense") - P5_3 = PTDF(sys5; slack_array=slack_array, linear_solver = "MKLPardiso") + P5_1 = PTDF(sys5; slack_array = slack_array, linear_solver = "KLU") + P5_2 = PTDF(sys5; slack_array = slack_array, linear_solver = "Dense") + P5_3 = PTDF(sys5; slack_array = 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 \ No newline at end of file + @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 From 65bb10ba3067bc2431bfba3a804773b03bc11eb7 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Mon, 31 Jul 2023 15:48:29 -0600 Subject: [PATCH 04/11] other error found --- test/test_ptdf.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index aeacd386..bcd96eb4 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -196,7 +196,7 @@ end @testset "PTDF matrices with distributed slack" begin sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") - bus_number = length(PNM.get_buses(sys)) + bus_number = length(PNM.get_buses(sys5)) dist_slack = 1 / bus_number * ones(bus_number) slack_array = dist_slack / sum(dist_slack) From 1b89faf3a7873fe6165cfaee8a7c5cc25e22adb9 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Mon, 31 Jul 2023 16:10:23 -0600 Subject: [PATCH 05/11] other error found 2 --- test/test_ptdf.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index bcd96eb4..937c2b6c 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -196,9 +196,9 @@ end @testset "PTDF matrices with distributed slack" begin sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") - bus_number = length(PNM.get_buses(sys5)) + buscount = length(PNM.get_buses(sys5)) - dist_slack = 1 / bus_number * ones(bus_number) + dist_slack = 1 / buscount * ones(buscount) slack_array = dist_slack / sum(dist_slack) slack_array = reshape(slack_array, 1, buscount) From 7e8ebee0e0ef08f66ff6e5268116aa9514970b7b Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 09:43:25 -0600 Subject: [PATCH 06/11] other error found 3 --- test/test_ptdf.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index 937c2b6c..8ae729ca 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -194,19 +194,20 @@ end 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) - slack_array = reshape(slack_array, 1, buscount) - P5_1 = PTDF(sys5; slack_array = slack_array, linear_solver = "KLU") - P5_2 = PTDF(sys5; slack_array = slack_array, linear_solver = "Dense") - P5_3 = PTDF(sys5; slack_array = slack_array, linear_solver = "MKLPardiso") + 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 From 4a07e6b18fe7c4ad20407c25dbada40aa8e19ee2 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 09:47:06 -0600 Subject: [PATCH 07/11] formatter --- test/test_ptdf.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index 8ae729ca..e5d33a1e 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -194,7 +194,6 @@ end end @testset "PTDF matrices with distributed slack" begin - sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") buscount = length(PNM.get_buses(sys5)) @@ -209,5 +208,4 @@ end @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 From 9a9255abbc89daa60cfb4de334fbf558c62ed390 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 11:27:44 -0600 Subject: [PATCH 08/11] other error found 4 --- src/ptdf_calculations.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl index 9f1f5312..e352159f 100644 --- a/src/ptdf_calculations.jl +++ b/src/ptdf_calculations.jl @@ -160,6 +160,8 @@ function _calculate_PTDF_matrix_KLU( else error("Distributed bus specification doesn't match the number of buses.") end + + return end """ @@ -260,7 +262,7 @@ function _calculate_PTDF_matrix_DENSE( error("Distributed bus specification doesn't match the number of buses.") end - return PTDFm_t + return end """ @@ -333,17 +335,18 @@ function _calculate_PTDF_matrix_MKLPardiso( ) elseif isempty(dist_slack) && length(ref_bus_positions) < buscount PTDFm_t[row_idx, :] = ABA_inv * @view BA[row_idx, :] + return PTDFm_t elseif length(dist_slack) == buscount @info "Distributed bus" PTDFm_t[row_idx, :] = ABA_inv * @view BA[row_idx, :] 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) + 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 """ From c4dde3f6c98462bd6fdd33561b3135fe929b31b9 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 12:09:57 -0600 Subject: [PATCH 09/11] other error found 5 --- src/ptdf_calculations.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/ptdf_calculations.jl b/src/ptdf_calculations.jl index e352159f..0e6b55d2 100644 --- a/src/ptdf_calculations.jl +++ b/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.", @@ -324,23 +325,23 @@ 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) + 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.") From 00c8b203784c4b2a3bba083b84fdffc480305842 Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 21:55:29 -0600 Subject: [PATCH 10/11] updated functions and added distributed slakc bus calculation for vPTDF --- src/common.jl | 2 +- src/virtual_ptdf_calculations.jl | 33 +++++++++++++++++++++++++------- test/test_virtual_ptdf.jl | 14 ++++++++++++++ 3 files changed, 41 insertions(+), 8 deletions(-) diff --git a/src/common.jl b/src/common.jl index c88f7f90..a31c392b 100644 --- a/src/common.jl +++ b/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 """ diff --git a/src/virtual_ptdf_calculations.jl b/src/virtual_ptdf_calculations.jl index f499ab7a..516bc1b4 100644 --- a/src/virtual_ptdf_calculations.jl +++ b/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 diff --git a/test/test_virtual_ptdf.jl b/test/test_virtual_ptdf.jl index e0c8612d..fb2d6c57 100644 --- a/test/test_virtual_ptdf.jl +++ b/test/test_virtual_ptdf.jl @@ -67,3 +67,17 @@ 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") + # compute full PTDF + ptdf = PTDF(sys; dist_slack = slack_array) + # compute each row of the virtual PTDF and compare values + vptdf = VirtualPTDF(sys; dist_slack = slack_array) + 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 From 48412248ee9915000eb990bd1b66c7509a958fbb Mon Sep 17 00:00:00 2001 From: alefcastelli Date: Tue, 1 Aug 2023 22:30:23 -0600 Subject: [PATCH 11/11] found some bugs --- src/virtual_lodf_calculations.jl | 5 ++--- test/test_virtual_ptdf.jl | 6 ++++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl index 5ae1ced7..314a24f5 100644 --- a/src/virtual_lodf_calculations.jl +++ b/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 diff --git a/test/test_virtual_ptdf.jl b/test/test_virtual_ptdf.jl index fb2d6c57..6ba5d0e6 100644 --- a/test/test_virtual_ptdf.jl +++ b/test/test_virtual_ptdf.jl @@ -71,10 +71,12 @@ 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 = slack_array) + ptdf = PTDF(sys; dist_slack = dist_slack) # compute each row of the virtual PTDF and compare values - vptdf = VirtualPTDF(sys; dist_slack = slack_array) + 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]