Unravelling the mystery of #861

[odow]

#861 needed to reduce some tolerances for a det check. That shouldn't have been necessary because computing the matrix of a 5x5 matrix with 17 non-zeros should have 1e-5 error!

This issue is just for me to track my notes while debugging.

Julia 1.6

julia> data = PowerModels.parse_file("/Users/Oscar/.julia/dev/PowerModels/test/data/matpower/case5_ext.m")
[warn | PowerModels]: active branch 8 is connected to inactive bus 11
[warn | PowerModels]: reversing the orientation of branch 6 (4, 3) to be consistent with other parallel branches
[warn | PowerModels]: active generators found at bus 4, updating to bus type from 1 to 2
[warn | PowerModels]: the voltage setpoint on generator 4 does not match the value at bus 4
[warn | PowerModels]: the voltage setpoint on generator 1 does not match the value at bus 1
[warn | PowerModels]: the voltage setpoint on generator 5 does not match the value at bus 10
[warn | PowerModels]: the voltage setpoint on generator 2 does not match the value at bus 1
[warn | PowerModels]: the voltage setpoint on generator 3 does not match the value at bus 3
[info | PowerModels]: removing 1 cost terms from generator 4: [4000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 1: [1400.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 5: [1000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 2: [1500.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 3: [3000.0, 0.0]
Dict{String, Any} with 13 entries:
  "bus"            => Dict{String, Any}("4"=>Dict{String, Any}("zone"=>1, "bus_i"=>4, "bus_type"=>2,…
  "source_type"    => "matpower"
  "name"           => "case5"
  "dcline"         => Dict{String, Any}()
  "source_version" => "2"
  "gen"            => Dict{String, Any}("4"=>Dict{String, Any}("pg"=>0.0, "model"=>2, "shutdown"=>0.…
  "branch"         => Dict{String, Any}("8"=>Dict{String, Any}("br_r"=>0.00297, "rate_a"=>4.26, "shi…
  "storage"        => Dict{String, Any}()
  "switch"         => Dict{String, Any}()
  "baseMVA"        => 100.0
  "per_unit"       => true
  "shunt"          => Dict{String, Any}("1"=>Dict{String, Any}("source_id"=>Any["bus", 11], "shunt_b…
  "load"           => Dict{String, Any}("1"=>Dict{String, Any}("source_id"=>Any["bus", 2], "load_bus…

julia> sm = PowerModels.calc_susceptance_matrix(data)
AdmittanceMatrix(5 buses, 17 entries)

julia> sm.matrix
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 17 stored entries:
 -222.507     35.2348     ⋅        32.569     154.703
   35.2348  -126.911    91.6758      ⋅           ⋅ 
     ⋅        91.6758  -94.7752     3.09933      ⋅ 
   32.569       ⋅        3.09933  -69.005      33.3367
  154.703       ⋅         ⋅        33.3367   -188.04

julia> LinearAlgebra.det(sm.matrix)
0.0

julia> sm.matrix.nzval
17-element Vector{Float64}:
 -222.5068568853514
   35.234840209999646
   32.56904637832204
  154.7029702970297
   35.234840209999646
 -126.91067446009131
   91.67583425009167
   91.67583425009167
  -94.77516156755509
    3.099327317463416
   32.56904637832204
    3.099327317463416
  -69.00504069581879
   33.33666700003334
  154.7029702970297
   33.33666700003334
 -188.03963729706305

Julia 1.9

julia> data = PowerModels.parse_file("/Users/Oscar/.julia/dev/PowerModels/test/data/matpower/case5_ext.m")
[warn | PowerModels]: active branch 8 is connected to inactive bus 11
[warn | PowerModels]: reversing the orientation of branch 6 (4, 3) to be consistent with other parallel branches
[warn | PowerModels]: active generators found at bus 4, updating to bus type from 1 to 2
[warn | PowerModels]: the voltage setpoint on generator 4 does not match the value at bus 4
[warn | PowerModels]: the voltage setpoint on generator 1 does not match the value at bus 1
[warn | PowerModels]: the voltage setpoint on generator 5 does not match the value at bus 10
[warn | PowerModels]: the voltage setpoint on generator 2 does not match the value at bus 1
[warn | PowerModels]: the voltage setpoint on generator 3 does not match the value at bus 3
[info | PowerModels]: removing 1 cost terms from generator 4: [4000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 1: [1400.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 5: [1000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 2: [1500.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 3: [3000.0, 0.0]
Dict{String, Any} with 13 entries:
  "bus"            => Dict{String, Any}("4"=>Dict{String, Any}("zone"=>1, "bus_i"=>4, "bus_type"=>2,…
  "source_type"    => "matpower"
  "name"           => "case5"
  "dcline"         => Dict{String, Any}()
  "source_version" => "2"
  "gen"            => Dict{String, Any}("4"=>Dict{String, Any}("vg"=>1.06414, "mbase"=>100.0, "sourc…
  "branch"         => Dict{String, Any}("8"=>Dict{String, Any}("br_r"=>0.00297, "rate_a"=>4.26, "shi…
  "storage"        => Dict{String, Any}()
  "switch"         => Dict{String, Any}()
  "baseMVA"        => 100.0
  "per_unit"       => true
  "shunt"          => Dict{String, Any}("1"=>Dict{String, Any}("source_id"=>Any["bus", 11], "shunt_b…
  "load"           => Dict{String, Any}("1"=>Dict{String, Any}("source_id"=>Any["bus", 2], "load_bus…

julia> sm = PowerModels.calc_susceptance_matrix(data)
AdmittanceMatrix(5 buses, 17 entries)

julia> sm.matrix
5×5 SparseArrays.SparseMatrixCSC{Float64, Int64} with 17 stored entries:
 -222.507     35.2348     ⋅        32.569     154.703
   35.2348  -126.911    91.6758      ⋅           ⋅ 
     ⋅        91.6758  -94.7752     3.09933      ⋅ 
   32.569       ⋅        3.09933  -69.005      33.3367
  154.703       ⋅         ⋅        33.3367   -188.04

julia> LinearAlgebra.det(sm.matrix)
-1.3498291219486257e-6

julia> sm.matrix.nzval
17-element Vector{Float64}:
 -222.50685688535137
   35.234840209999646
   32.56904637832204
  154.70297029702968
   35.234840209999646
 -126.91067446009131
   91.67583425009167
   91.67583425009167
  -94.7751615675551
    3.0993273174634197
   32.56904637832204
    3.0993273174634197
  -69.0050406958188
   33.33666700003334
  154.70297029702968
   33.33666700003334
 -188.03963729706302

Diffing the two vectors shows

julia> nzval_16 = [-222.5068568853514
          35.234840209999646
          32.56904637832204
         154.7029702970297
          35.234840209999646
        -126.91067446009131
          91.67583425009167
          91.67583425009167
         -94.77516156755509
           3.099327317463416
          32.56904637832204
           3.099327317463416
         -69.00504069581879
          33.33666700003334
         154.7029702970297
          33.33666700003334
        -188.03963729706305
       ]
17-element Vector{Float64}:
 -222.5068568853514
   35.234840209999646
   32.56904637832204
  154.7029702970297
   35.234840209999646
 -126.91067446009131
   91.67583425009167
   91.67583425009167
  -94.77516156755509
    3.099327317463416
   32.56904637832204
    3.099327317463416
  -69.00504069581879
   33.33666700003334
  154.7029702970297
   33.33666700003334
 -188.03963729706305

julia> nzval_16 .- sm.matrix.nzval
17-element Vector{Float64}:
 -2.842170943040401e-14
  0.0
  0.0
  2.842170943040401e-14
  0.0
  0.0
  0.0
  0.0
  1.4210854715202004e-14
 -3.552713678800501e-15
  0.0
 -3.552713678800501e-15
  1.4210854715202004e-14
  0.0
  2.842170943040401e-14
  0.0
 -2.842170943040401e-14

[odow]

So Julia 1.9 is noticeably less-accurate at competing inv of a complex number. And this causes problems with the susceptance matrix. It's probably also the cause of the numeric issue on Windows.

Julia 1.9

julia> imag(inv(0.00064 + 0.0064im))
-154.70297029702968

julia> x = imag(inv(big"0.00064" + big"0.0064" * im))
-154.7029702970297029702970297029702970297029702970297029702970297029702970297017

julia> x - imag(inv(0.00064 + 0.0064im))
-2.335645428439141205041715414217202970297029702970297029702970170069474548472957e-14

Julia 1.6

julia> imag(inv(0.00064 + 0.0064im))
-154.7029702970297

julia> x = imag(inv(big"0.00064" + big"0.0064" * im))
-154.7029702970297029702970297029702970297029702970297029702970297029702970297017

julia> x - imag(inv(0.00064 + 0.0064im))
5.065255146012595384427816560952970297029702970297029702970298299305254515270434e-15

Diff

julia> -154.70297029702968 - -154.7029702970297
2.842170943040401e-14

[odow]

Upstream change seems to be JuliaLang/julia#47255

[odow]

The two Float64s are neighbors

julia> prevfloat(-154.70297029702968)
-154.7029702970297

so it seems unlikely that they'll revert the upstream change. I wonder if we should compute the matrix in better precision on the PowerModels side.