diff --git a/src/lodf_calculations.jl b/src/lodf_calculations.jl index b69fb6bd..9a96a9ab 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 @@ -147,7 +147,50 @@ function _calculate_LODF_matrix_DENSE( push!(m_V, 1.0 - ptdf_denominator_t[iline, iline]) end end - lodf_t = LinearAlgebra.diagm(m_V) \ (a * ptdf) + 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_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 + 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 + + # 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 @@ -160,6 +203,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_lodf.jl b/test/test_lodf.jl index 5a26cc2a..8b6bea15 100644 --- a/test/test_lodf.jl +++ b/test/test_lodf.jl @@ -30,12 +30,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) @@ -63,11 +67,9 @@ end @test test_value - # 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) + L5NS_from_ba_aba = LODF(A, P5; linear_solver = "XXX") catch err if err isa ErrorException test_value = true @@ -75,10 +77,11 @@ end @test test_value - # test if error is thrown in case `MKLPardiso` is chosen as a linera solver + # test if error is thrown in case `tol` is defined in PTDF + P5 = PTDF(sys5; tol = 1e-3) test_value = false try - lodf = LODF(sys5; linear_solver = "MKLPardiso") + L5NS_from_ptdf = LODF(A, P5) catch err if err isa ErrorException test_value = true