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/ptdf_calculations.jl b/src/ptdf_calculations.jl index 83074943..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.", @@ -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 """ 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/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_ptdf.jl b/test/test_ptdf.jl index ccb24643..e5d33a1e 100644 --- a/test/test_ptdf.jl +++ b/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 diff --git a/test/test_virtual_ptdf.jl b/test/test_virtual_ptdf.jl index e0c8612d..6ba5d0e6 100644 --- a/test/test_virtual_ptdf.jl +++ b/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