Skip to content

Commit

Permalink
Merge pull request #60 from NREL-Sienna/jd/debug_LODF
Browse files Browse the repository at this point in the history 
Jd/debug lodf
  • Loading branch information
@jd-lara
jd-lara authored Oct 7, 2023
2 parents fea6218 + 74e093f commit a13f879
Show file tree
Hide file tree
Showing 5 changed files with 142 additions and 36 deletions.
2 changes: 2 additions & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
KLU = "ef3ab10e-7fda-4108-b977-705223b18434"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Expand All @@ -17,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
DocStringExtensions = "~0.8, ~0.9"
HDF5 = "0.17"
InfrastructureSystems = "^1.20"
MKL = "0.6"
KLU = "~0.4"
Pardiso = "~0.5.4"
PowerSystems = "3"
Expand Down
6 changes: 4 additions & 2 deletions src/PowerNetworkMatrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ export VirtualPTDF
export Ybus

using DocStringExtensions
import MKL
import InfrastructureSystems
import PowerSystems
import PowerSystems: ACBusTypes
Expand All @@ -31,9 +32,10 @@ import SparseArrays: rowvals, nzrange
import HDF5
import KLU: klu
import KLU
import LinearAlgebra: LAPACK.getri!, LAPACK.getrf!, BLAS.gemm
import LinearAlgebra: ldiv!, mul!, I, dot
import LinearAlgebra
import LinearAlgebra: BLAS.gemm
import LinearAlgebra: ldiv!, mul!, I, dot
import LinearAlgebra: LAPACK.getrf!, LAPACK.getrs!
import Pardiso

@template DEFAULT = """
Expand Down
2 changes: 2 additions & 0 deletions src/definitions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,5 @@ const MiB = KiB * KiB
const GiB = MiB * KiB
const MAX_CACHE_SIZE_MiB = 100
const ROW_PERSISTENT_CACHE_WARN = 1 * GiB

DEFAULT_LODF_CHUNK_SIZE = 18_000
131 changes: 104 additions & 27 deletions src/lodf_calculations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,9 +80,9 @@ function _calculate_LODF_matrix_KLU(
end
end
Dem_LU = klu(SparseArrays.sparse(m_I, m_I, m_V))
lodf_t = Dem_LU \ ptdf_denominator
lodf_t[SparseArrays.diagind(lodf_t)] .= -1.0
return lodf_t
KLU.solve!(Dem_LU, ptdf_denominator)
ptdf_denominator[SparseArrays.diagind(ptdf_denominator)] .= -1.0
return ptdf_denominator
end

function _calculate_LODF_matrix_KLU(
Expand Down Expand Up @@ -123,9 +123,101 @@ function _calculate_LODF_matrix_DENSE(
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
(mV, bipiv, binfo) = getrf!(Matrix(LinearAlgebra.diagm(m_V)))
_binfo_check(binfo)
getrs!('N', mV, bipiv, ptdf_denominator_t)
ptdf_denominator_t[LinearAlgebra.diagind(ptdf_denominator_t)] .= -1.0
return ptdf_denominator_t
end

function _pardiso_sequential_LODF!(
lodf_t::Matrix{Float64},
A::SparseArrays.SparseMatrixCSC{Float64, Int},
ptdf_denominator_t::Matrix{Float64},
chunk_size::Int = DEFAULT_LODF_CHUNK_SIZE,
)
@info "Line Count too large for single compute using Pardiso. Employing Sequential Calulations using a chunk_size=$(chunk_size)"
linecount = size(lodf_t, 1)
@assert LinearAlgebra.ishermitian(A)
ps = Pardiso.MKLPardisoSolver()
Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
Pardiso.pardisoinit(ps)
# Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
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, 12, 1)
Pardiso.set_iparm!(ps, 11, 0)
Pardiso.set_iparm!(ps, 13, 0)
Pardiso.set_iparm!(ps, 32, 1)
#Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
Pardiso.set_phase!(ps, Pardiso.ANALYSIS)
Pardiso.pardiso(
ps,
lodf_t,
A,
ptdf_denominator_t,
)

Pardiso.set_phase!(ps, Pardiso.NUM_FACT)
Pardiso.pardiso(
ps,
A,
Float64[],
)
Pardiso.set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
i_count = 1
tmp = zeros(Float64, linecount, chunk_size)
while i_count <= linecount
edge = min(i_count + chunk_size - 1, linecount)
if linecount - edge <= 0
tmp = tmp[:, 1:(edge - i_count + 1)]
end
Pardiso.pardiso(
ps,
tmp,
A,
ptdf_denominator_t[:, i_count:edge],
)
lodf_t[:, i_count:edge] .= tmp
i_count = edge + 1
end
Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
Pardiso.pardiso(ps)
return
end

function _pardiso_single_LODF!(
lodf_t::Matrix{Float64},
A::SparseArrays.SparseMatrixCSC{Float64, Int},
ptdf_denominator_t::Matrix{Float64},
)
@assert LinearAlgebra.ishermitian(A)
ps = Pardiso.MKLPardisoSolver()
Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM_POSDEF)
Pardiso.pardisoinit(ps)
Pardiso.set_iparm!(ps, 1, 1)
defaults = Pardiso.get_iparms(ps)
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, 12, 1)
#Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
Pardiso.pardiso(
ps,
lodf_t,
A,
ptdf_denominator_t,
)
Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
Pardiso.pardiso(ps)
return
end

function _calculate_LODF_matrix_MKLPardiso(
Expand All @@ -145,28 +237,13 @@ function _calculate_LODF_matrix_MKLPardiso(
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
A = SparseArrays.sparse(m_I, m_I, m_V)
if linecount > DEFAULT_LODF_CHUNK_SIZE
_pardiso_sequential_LODF!(lodf_t, A, ptdf_denominator_t)
else
_pardiso_single_LODF!(lodf_t, A, ptdf_denominator_t)
end
lodf_t[LinearAlgebra.diagind(lodf_t)] .= -1.0
return lodf_t
end
Expand Down
37 changes: 30 additions & 7 deletions src/ptdf_calculations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,19 +123,17 @@ function _calculate_PTDF_matrix_KLU(
# inizialize matrices for evaluation
valid_ix = setdiff(1:buscount, ref_bus_positions)
PTDFm_t = zeros(buscount, linecount)

copyto!(PTDFm_t, BA)
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
copyto!(PTDFm_t, BA)
PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :])
PTDFm_t[collect(ref_bus_positions), :] .= 0.0
return PTDFm_t
elseif length(dist_slack) == buscount
@info "Distributed bus"
copyto!(PTDFm_t, BA)
PTDFm_t[valid_ix, :] = KLU.solve!(K, PTDFm_t[valid_ix, :])
PTDFm_t[collect(ref_bus_positions), :] .= 0.0
slack_array = dist_slack / sum(dist_slack)
Expand Down Expand Up @@ -172,6 +170,19 @@ function calculate_PTDF_matrix_KLU(
return PTDFm, A
end

function _binfo_check(binfo::Int)
if binfo != 0
if binfo < 0
error("Illegal Argument in Inputs")
elseif binfo > 0
error("Singular value in factorization. Possibly there is an islanded bus")
else
@assert false
end
end
return
end

"""
Funciton for internal use only.
Expand All @@ -198,17 +209,21 @@ function _calculate_PTDF_matrix_DENSE(
valid_ixs = setdiff(1:buscount, ref_bus_positions)
ABA = Matrix(calculate_ABA_matrix(A, BA, ref_bus_positions))
PTDFm_t = zeros(buscount, linecount)

(ABA, bipiv, binfo) = getrf!(ABA)
_binfo_check(binfo)
BA = Matrix(BA[valid_ixs, :])
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[valid_ixs, :] = ABA \ BA[valid_ixs, :]
getrs!('N', ABA, bipiv, BA)
PTDFm_t[valid_ixs, :] = BA
return PTDFm_t
elseif length(dist_slack) == buscount
@info "Distributed bus"
PTDFm_t[valid_ixs, :] = ABA \ BA[valid_ixs, :]
getrs!('N', ABA, bipiv, BA)
PTDFm_t[valid_ixs, :] = BA
slack_array = dist_slack / sum(dist_slack)
slack_array = reshape(slack_array, 1, buscount)
return PTDFm_t -
Expand Down Expand Up @@ -268,8 +283,10 @@ function _calculate_PTDF_matrix_MKLPardiso(
buscount = size(BA, 1)

ABA = calculate_ABA_matrix(A, BA, ref_bus_positions)
# Here add the subnetwork detection
@assert LinearAlgebra.issymmetric(ABA)
ps = Pardiso.MKLPardisoSolver()
Pardiso.set_matrixtype!(ps, Pardiso.REAL_SYM)
Pardiso.pardisoinit(ps)
# Pardiso.set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
defaults = Pardiso.get_iparms(ps)
Pardiso.set_iparm!(ps, 1, 1)
Expand All @@ -279,6 +296,10 @@ function _calculate_PTDF_matrix_MKLPardiso(
Pardiso.set_iparm!(ps, 2, 2)
Pardiso.set_iparm!(ps, 59, 2)
Pardiso.set_iparm!(ps, 6, 1)
Pardiso.set_iparm!(ps, 12, 1)
Pardiso.set_iparm!(ps, 11, 0)
Pardiso.set_iparm!(ps, 13, 0)
Pardiso.set_iparm!(ps, 32, 1)

# inizialize matrices for evaluation
valid_ix = setdiff(1:buscount, ref_bus_positions)
Expand All @@ -293,12 +314,14 @@ function _calculate_PTDF_matrix_MKLPardiso(
Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
PTDFm_t[valid_ix, :] = full_BA
Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
Pardiso.pardiso(ps)
return PTDFm_t
elseif length(dist_slack) == buscount
@info "Distributed bus"
Pardiso.pardiso(ps, PTDFm_t[valid_ix, :], ABA, full_BA)
PTDFm_t[valid_ix, :] = full_BA
Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
Pardiso.pardiso(ps)
slack_array = dist_slack / sum(dist_slack)
slack_array = reshape(slack_array, 1, buscount)
return PTDFm_t - ones(buscount, 1) * (slack_array * PTDFm_t)
Expand Down

0 comments on commit a13f879

Please sign in to comment.