Faster PTDF matrix computation (#849)

* Faster PTDF matrix computation
* Move factorization to calc_susceptance_matrix_inv; use matrix-matrix product in PTDF
* update changelog
* Add mtanneau to contributors
* Update documentation

Co-authored-by: mtanneau3 <[email protected]>

CHANGELOG.md

@@ -3,6 +3,7 @@ PowerModels.jl Change Log
 
 ### Staged
 - Improved network data model documentation
+- Faster PTDF matrix computation (#849)
 
 ### v0.19.7
 - Improve linear algebra test robustness (#827)

README.md

@@ -65,6 +65,7 @@ The primary developer is Carleton Coffrin (@ccoffrin) with support from the foll
 - Miles Lubin (@mlubin) MIT, Julia/JuMP advise
 - Yeesian Ng (@yeesian) MIT, Documenter.jl setup
 - Kaarthik Sundar (@kaarthiksundar) LANL, OBBT utility
+- Mathieu Tanneau (@mtanneau) Georgia Tech, PTDF matrix computation
 - Byron Tasseff (@tasseff) LANL, multi-infrastructure updates
 
 

src/core/admittance_matrix.jl

@@ -136,16 +136,41 @@ end
 Base.show(io::IO, x::AdmittanceMatrixInverse{<:Number}) = print(io, "AdmittanceMatrixInverse($(length(x.idx_to_bus)) buses, $(length(x.matrix)) entries)")
 
 
-"note, data should be a PowerModels network data model; only supports networks with exactly one refrence bus"
+"""
+    calc_susceptance_matrix_inv(data)
+
+Compute the inverse of the network's susceptance matrix.
+
+Note: `data`` should be a PowerModels network data model; only supports networks with exactly one refrence bus.
+
+While the susceptance matrix is sparse, its inverse it typically quite dense.
+This implementation first computes a sparse factorization, then recovers the (dense)
+    matrix inverse via backward substitution. This is more efficient
+    than directly computing a dense inverse with `LinearAlgebra.inv`.
+"""
 function calc_susceptance_matrix_inv(data::Dict{String,<:Any})
     #TODO check single connected component
-
     sm = calc_susceptance_matrix(data)
-
+    S  = sm.matrix
+    num_buses = length(sm.idx_to_bus)  # this avoids inactive buses
+
     ref_bus = reference_bus(data)
-    sm_inv = calc_admittance_matrix_inv(sm, sm.bus_to_idx[ref_bus["index"]])
-
-    return sm_inv
+    ref_idx = sm.bus_to_idx[ref_bus["index"]]
+    if !(ref_idx > 0 && ref_idx <= num_buses)
+        Memento.error(_LOGGER, "invalid ref_idx in calc_susceptance_matrix_inv")
+    end
+    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
 
 "calculates the inverse of the susceptance matrix"

src/core/data_basic.jl

@@ -322,8 +322,6 @@ function compute_basic_dc_pf(data::Dict{String,<:Any})
     return theta
 end
 
-
-
 """
 given a basic network data dict, returns a real valued ptdf matrix with one
 row for each branch and one column for each bus in the network.
@@ -335,25 +333,15 @@ function calc_basic_ptdf_matrix(data::Dict{String,<:Any})
         Memento.warn(_LOGGER, "calc_basic_ptdf_matrix requires basic network data and given data may be incompatible. make_basic_network can be used to transform data into the appropriate form.")
     end
 
-    num_bus = length(data["bus"])
-    num_branch = length(data["branch"])
-
     b_inv = calc_susceptance_matrix_inv(data).matrix
 
-    ptdf = zeros(num_branch, num_bus)
-    for (i,branch) in data["branch"]
-        branch_idx = branch["index"]
-        bus_fr = branch["f_bus"]
-        bus_to = branch["t_bus"]
-        g,b = calc_branch_y(branch)
-        for n in 1:num_bus
-            ptdf[branch_idx, n] = b*(b_inv[bus_fr, n] - b_inv[bus_to, n])
-        end
-    end
+    B = calc_basic_branch_susceptance_matrix(data)
+    ptdf = B * b_inv
 
     return ptdf
 end
 
+
 """
 given a basic network data dict and a branch index returns a row of the ptdf
 matrix reflecting that branch.