Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Compute the susceptance matrix in higher precision #863

Closed
wants to merge 10 commits into from
Closed

Compute the susceptance matrix in higher precision #863

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
4936aea
Compute the susceptance matrix in higher precision
odow Jun 1, 2023
ae59ca5
Update test/pf.jl
odow Jun 1, 2023
8b194bc
Another case of inv(::ComplexF64)
odow Jun 1, 2023
8c8928c
Update admittance_matrix.jl
odow Jun 1, 2023
4909adc
Update test/pf.jl
odow Jun 1, 2023
678ba6d
Update test/pf.jl
odow Jun 1, 2023
b10a40a
Update admittance_matrix.jl
odow Jun 4, 2023
716cf2d
Apply suggestions from code review
odow Jun 4, 2023
c7d0473
Apply suggestions from code review
odow Jun 5, 2023
174194a
Update test/pf.jl
odow Jun 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 9 additions & 9 deletions src/core/admittance_matrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,15 @@ function calc_admittance_matrix(data::Dict{String,<:Any})
if branch[pm_component_status["branch"]] != pm_component_status_inactive["branch"] && haskey(bus_to_idx, f_bus) && haskey(bus_to_idx, t_bus)
f_bus = bus_to_idx[f_bus]
t_bus = bus_to_idx[t_bus]
y = inv(branch["br_r"] + branch["br_x"]im)
y = inv(big(branch["br_r"] + branch["br_x"]im))
tr, ti = calc_branch_t(branch)
t = tr + ti*im
lc_fr = branch["g_fr"] + branch["b_fr"]im
lc_to = branch["g_to"] + branch["b_to"]im
push!(I, f_bus); push!(J, t_bus); push!(V, -y/conj(t))
push!(I, t_bus); push!(J, f_bus); push!(V, -(y/t))
push!(I, f_bus); push!(J, f_bus); push!(V, (y + lc_fr)/abs2(t))
push!(I, t_bus); push!(J, t_bus); push!(V, (y + lc_to))
push!(I, f_bus); push!(J, t_bus); push!(V, convert(ComplexF64, -y/conj(t)))
push!(I, t_bus); push!(J, f_bus); push!(V, convert(ComplexF64, -(y/t)))
push!(I, f_bus); push!(J, f_bus); push!(V, convert(ComplexF64, (y + lc_fr)/abs2(t)))
push!(I, t_bus); push!(J, t_bus); push!(V, convert(ComplexF64, (y + lc_to)))
end
end

Expand Down Expand Up @@ -104,7 +104,7 @@ function calc_susceptance_matrix(data::Dict{String,<:Any})
if branch[pm_component_status["branch"]] != pm_component_status_inactive["branch"] && haskey(bus_to_idx, f_bus) && haskey(bus_to_idx, t_bus)
f_bus = bus_to_idx[f_bus]
t_bus = bus_to_idx[t_bus]
b_val = imag(inv(branch["br_r"] + branch["br_x"]im))
b_val = convert(Float64, imag(inv(big(branch["br_r"] + branch["br_x"]im))))
push!(I, f_bus); push!(J, t_bus); push!(V, -b_val)
push!(I, t_bus); push!(J, f_bus); push!(V, -b_val)
push!(I, f_bus); push!(J, f_bus); push!(V, b_val)
Expand Down Expand Up @@ -153,7 +153,7 @@ function calc_susceptance_matrix_inv(data::Dict{String,<:Any})
sm = calc_susceptance_matrix(data)
S = sm.matrix
num_buses = length(sm.idx_to_bus) # this avoids inactive buses

ref_bus = reference_bus(data)
ref_idx = sm.bus_to_idx[ref_bus["index"]]
if !(ref_idx > 0 && ref_idx <= num_buses)
Expand All @@ -162,14 +162,14 @@ function calc_susceptance_matrix_inv(data::Dict{String,<:Any})
S[ref_idx, :] .= 0.0
S[:, ref_idx] .= 0.0
S[ref_idx, ref_idx] = 1.0

F = LinearAlgebra.ldlt(Symmetric(S); check=false)
if !LinearAlgebra.issuccess(F)
Memento.error(_LOGGER, "Failed factorization in calc_susceptance_matrix_inv")
end
M = F \ Matrix(1.0I, num_buses, num_buses)
M[ref_idx, :] .= 0.0 # zero-out the row of the slack bus

return AdmittanceMatrixInverse(sm.idx_to_bus, sm.bus_to_idx, ref_idx, M)
end

Expand Down
4 changes: 2 additions & 2 deletions test/am.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,15 +42,15 @@ end

@test isa(sm, AdmittanceMatrix{Float64})
@test SparseArrays.nnz(sm.matrix) == 17
@test isapprox(LinearAlgebra.det(sm.matrix), 0.0, atol=1e-5)
@test isapprox(LinearAlgebra.det(sm.matrix), 0.0, atol=1e-6)
end
@testset "5-bus ext case" begin
data = PowerModels.parse_file("../test/data/matpower/case5_ext.m")
sm = calc_susceptance_matrix(data)

@test isa(sm, AdmittanceMatrix{Float64})
@test SparseArrays.nnz(sm.matrix) == 17
@test isapprox(LinearAlgebra.det(sm.matrix), 0.0, atol=1e-5)
@test isapprox(LinearAlgebra.det(sm.matrix), 0.0, atol=1e-6)
end
@testset "14-bus pti case" begin
data = PowerModels.parse_file("../test/data/pti/case14.raw")
Expand Down
10 changes: 3 additions & 7 deletions test/pf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,14 +81,10 @@ end

@testset "test ac rect pf" begin
@testset "5-bus asymmetric case" begin
nlp_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-6)
result = run_pf("../test/data/matpower/case5_asym.m", ACRPowerModel, nlp_solver)
if VERSION >= v"1.9" && Sys.iswindows()
# Some numerical issue on Windows with Julia 1.9?
@test result["termination_status"] in (LOCALLY_SOLVED, ITERATION_LIMIT)
else
@test result["termination_status"] == LOCALLY_SOLVED
@test isapprox(result["objective"], 0; atol = 1e-2)
end
@test result["termination_status"] == LOCALLY_SOLVED
@test isapprox(result["objective"], 0; atol = 1e-2)
end
#=
# numerical issues with ipopt, likely div. by zero issue in jacobian
Expand Down