diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl index 4785383d..79c5c061 100644 --- a/src/lodf_calculations.jl +++ b/src/lodf_calculations.jl @@ -32,8 +32,8 @@ function _buildlodf( lodf_t = _calculate_LODF_matrix_KLU(a, ptdf) elseif linear_solver == "Dense" lodf_t = _calculate_LODF_matrix_DENSE(a, ptdf) - else - error("MKLPardiso still to be implemented") + elseif linear_solver == "MKLPardiso" + lodf_t = _calculate_LODF_matrix_MKLPardiso(a, ptdf) end return lodf_t end @@ -109,50 +109,65 @@ function _calculate_LODF_matrix_KLU( return lodf_t end -# ! temp for evaluation -function _calculate_LODF_matrix_KLU2( +function _calculate_LODF_matrix_DENSE( a::SparseArrays.SparseMatrixCSC{Int8, Int}, ptdf::Matrix{Float64}, ) - ptdf_denominator_t = a * ptdf linecount = size(ptdf, 2) - lodf_t = zeros(linecount, linecount) - for i in 1:linecount - for j in 1:linecount - if i == j - lodf_t[i, i] = -1.0 - else - if (1.0 - ptdf_denominator_t[i, i]) < 1.0E-06 - lodf_t[i, j] = ptdf_denominator_t[i, j] - else - lodf_t[i, j] = ptdf_denominator_t[i, j] / (1 - ptdf_denominator_t[i, i]) - end - end + ptdf_denominator_t = a * ptdf + m_V = Float64[] + for iline in 1:linecount + if (1.0 - ptdf_denominator_t[iline, iline]) < 1.0E-06 + push!(m_V, 1.0) + else + push!(m_V, 1.0 - ptdf_denominator_t[iline, iline]) end end + lodf_t = LinearAlgebra.diagm(m_V) \ ptdf_denominator_t + lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0 return lodf_t end -function _calculate_LODF_matrix_DENSE( +function _calculate_LODF_matrix_MKLPardiso( a::SparseArrays.SparseMatrixCSC{Int8, Int}, ptdf::Matrix{Float64}, ) linecount = size(ptdf, 2) ptdf_denominator_t = a * ptdf + m_I = Int[] + m_V = Float64[] for iline in 1:linecount if (1.0 - ptdf_denominator_t[iline, iline]) < 1.0E-06 - ptdf_denominator_t[iline, iline] = 0.0 + push!(m_I, iline) + push!(m_V, 1.0) + else + push!(m_I, iline) + push!(m_V, 1 - ptdf_denominator_t[iline, iline]) end end - (Dem, dipiv, dinfo) = getrf!( - Matrix{Float64}(1.0I, linecount, linecount) - - Array(LinearAlgebra.Diagonal(ptdf_denominator_t)), - ) - lodf_t = gemm('N', 'N', getri!(Dem, dipiv), ptdf_denominator_t) - lodf_t = - lodf_t - Array(LinearAlgebra.Diagonal(lodf_t)) - - Matrix{Float64}(1.0I, linecount, linecount) + # intialize MKLPardiso + ps = Pardiso.MKLPardisoSolver() + defaults = Pardiso.get_iparms(ps) + Pardiso.set_iparm!(ps, 1, 1) + for (ix, v) in enumerate(defaults[2:end]) + Pardiso.set_iparm!(ps, ix + 1, v) + end + Pardiso.set_iparm!(ps, 2, 2) + Pardiso.set_iparm!(ps, 59, 2) + Pardiso.set_iparm!(ps, 6, 1) + # inizialize matrix for evaluation + lodf_t = zeros(linecount, linecount) + # solve system + Pardiso.pardiso( + ps, + lodf_t, + SparseArrays.sparse(m_I, m_I, m_V), + ptdf_denominator_t, + ) + Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL) + # set diagonal to -1 + lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0 return lodf_t end @@ -164,6 +179,10 @@ Builds the LODF matrix given the data of branches and buses of the system. vector of the System AC branches - `buses::Vector{PSY.ACBus}`: vector of the System buses +- `linear_solver`::String + linear solver to use for matrix evaluation. + Available options: "KLU", "Dense" (OpenBLAS) or "MKLPardiso". + Default solver: "KLU". - `tol::Float64`: Tolerance to eliminate entries in the LODF matrix (default eps()) """ diff --git a/test/test_BA_ABA_matrix.jl b/test/test_BA_ABA_matrix.jl index a7b19911..eb728061 100644 --- a/test/test_BA_ABA_matrix.jl +++ b/test/test_BA_ABA_matrix.jl @@ -67,17 +67,7 @@ end # test if error is correctly thrown when ref bus is called rb = collect(ba.ref_bus_positions)[1] - test_val = false - try - aba[rb, :] - catch err - if err isa ErrorException - test_val = true - else - error("Expected an ErrorException but was not thrown.") - end - end - @test test_val + @test_throws ErrorException aba[rb, :] end @testset "Test show for A, BA and ABA matrix" begin diff --git a/test/test_lodf.jl b/test/test_lodf.jl index 5a26cc2a..0dec2892 100644 --- a/test/test_lodf.jl +++ b/test/test_lodf.jl @@ -13,12 +13,6 @@ @test isapprox(maximum(L5.data), maximum(Lodf_5), atol = 1e-3) @test isapprox(L5[branches_5[1], branches_5[2]], 0.3447946513849093) - # get LODF with second method - a = IncidenceMatrix(sys5) - ptdf = PTDF(sys5) - lodf_t_3 = PNM._calculate_LODF_matrix_KLU2(a.data, ptdf.data) - @test isapprox(lodf_t_3, L5.data, atol = 1e-5) - # get LODF with system L5NS = LODF(sys5) @test getindex(L5NS, "5", "6") - -0.3071 <= 1e-4 @@ -30,12 +24,16 @@ P5 = PTDF(sys5) L5NS_from_ptdf = LODF(A, P5) L5NS_from_ptdf2 = LODF(A, P5; linear_solver = "Dense") + L5NS_from_ptdf3 = LODF(A, P5; linear_solver = "MKLPardiso") @test getindex(L5NS_from_ptdf, "5", "6") - -0.3071 <= 1e-4 @test getindex(L5NS_from_ptdf2, "5", "6") - -0.3071 <= 1e-4 + @test getindex(L5NS_from_ptdf3, "5", "6") - -0.3071 <= 1e-4 total_error = abs.(L5NS_from_ptdf.data' .- Lodf_5) total_error2 = abs.(L5NS_from_ptdf2.data' .- Lodf_5) + total_error3 = abs.(L5NS_from_ptdf3.data' .- Lodf_5) @test isapprox(sum(total_error), 0.0, atol = 1e-3) @test isapprox(sum(total_error2), 0.0, atol = 1e-3) + @test isapprox(sum(total_error3), 0.0, atol = 1e-3) # A, ABA, and BA case ABA = ABA_Matrix(sys5; factorize = true) @@ -53,38 +51,14 @@ end # test if error is thrown in case other linear solvers are called - test_value = false - try - L5NS_from_ba_aba = LODF(A, ABA, BA; linear_solver = "Dense") - catch err - if err isa ErrorException - test_value = true - end - end - @test test_value + @test_throws ErrorException LODF(A, ABA, BA; linear_solver = "Dense") + + @test_throws ErrorException LODF(A, P5; linear_solver = "XXX") # test if error is thrown in case `tol` is defined in PTDF P5 = PTDF(sys5; tol = 1e-3) test_value = false - try - L5NS_from_ptdf = LODF(A, P5) - catch err - if err isa ErrorException - test_value = true - end - end - @test test_value - - # test if error is thrown in case `MKLPardiso` is chosen as a linera solver - test_value = false - try - lodf = LODF(sys5; linear_solver = "MKLPardiso") - catch err - if err isa ErrorException - test_value = true - end - end - @test test_value + @test_throws ErrorException L5NS_from_ptdf = LODF(A, P5) # get 14 bus system buses_14 = nodes14() diff --git a/test/test_ptdf.jl b/test/test_ptdf.jl index 882affa5..0a06ce9f 100644 --- a/test/test_ptdf.jl +++ b/test/test_ptdf.jl @@ -230,37 +230,13 @@ end dist_slack = 1 / buscount * ones(buscount) slack_array = dist_slack / sum(dist_slack) - test_val1 = false - try - ptdf_1 = PTDF(sys; dist_slack = slack_array) - catch err - if err isa ErrorException - test_val1 = true - else - error( - "Error was not thrown for PTDF with distributed slack bus and more than one reference bus.", - ) - end - end + @test_throws ErrorException ptdf_1 = PTDF(sys; dist_slack = slack_array) # incorrect dist_slack arrya length sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") buscount = length(PNM.get_buses(sys5)) + 1 dist_slack = 1 / buscount * ones(buscount) slack_array = dist_slack / sum(dist_slack) - test_val2 = false - try - ptdf_2 = PTDF(sys5; dist_slack = slack_array) - catch err - if err isa ErrorException - test_val2 = true - else - error( - "Error was not thrown for PTDF with distributed slack array of incorrect length.", - ) - end - end - @test test_val1 - @test test_val2 + @test_throws ErrorException ptdf_2 = PTDF(sys5; dist_slack = slack_array) end diff --git a/test/test_virtual_lodf.jl b/test/test_virtual_lodf.jl index c5811f7e..9dd37ed9 100644 --- a/test/test_virtual_lodf.jl +++ b/test/test_virtual_lodf.jl @@ -112,15 +112,7 @@ end length(axes(vlodf)[1]) * length(axes(vlodf)[2]) # check if error is correctly thrown - test_value = false - try - vlodf[1, 1] = 1 - catch err - if err isa ErrorException - test_value = true - end - end - @test test_value + @test_throws ErrorException vlodf[1, 1] = 1 # test show test_value = false diff --git a/test/test_virtual_ptdf.jl b/test/test_virtual_ptdf.jl index 3ec4b509..69cffd98 100644 --- a/test/test_virtual_ptdf.jl +++ b/test/test_virtual_ptdf.jl @@ -101,19 +101,8 @@ end dist_slack = 1 / buscount * ones(buscount) slack_array = dist_slack / sum(dist_slack) - test_val1 = false - try - ptdf_1 = VirtualPTDF(sys; dist_slack = slack_array) - ptdf_1[1, 1] - catch err - if err isa ErrorException - test_val1 = true - else - error( - "Error was not thrown for PTDF with distributed slack bus and more than one reference bus.", - ) - end - end + ptdf_1 = VirtualPTDF(sys; dist_slack = slack_array) + @test_throws ErrorException ptdf_1[1, 1] # incorrect dist_slack array length sys5 = PSB.build_system(PSB.PSITestSystems, "c_sys5") @@ -121,21 +110,8 @@ end dist_slack = 1 / buscount * ones(buscount) slack_array = dist_slack / sum(dist_slack) test_val2 = false - try - ptdf_2 = VirtualPTDF(sys5; dist_slack = slack_array) - ptdf_2[1, 1] - catch err - if err isa ErrorException - test_val2 = true - else - error( - "Error was not thrown for PTDF with distributed slack array of incorrect length.", - ) - end - end - - @test test_val1 - @test test_val2 + ptdf_2 = VirtualPTDF(sys5; dist_slack = slack_array) + @test_throws ErrorException ptdf_2[1, 1] end @testset "Test Virtual PTDF auxiliary functions" begin @@ -150,15 +126,7 @@ end @test isempty(vptdf) == false # check if error is correctly thrown - test_value = false - try - vptdf[1, 1] = 1 - catch err - if err isa ErrorException - test_value = true - end - end - @test test_value + @test_throws ErrorException vptdf[1, 1] = 1 # get the rows and full PTDF matrix, test get_ptdf_data ptdf = PTDF(sys)