GitHub

How to Build Multiple Matrices Without Repeating Work

To follow along, you can download this tutorial as a Julia script (.jl) or Jupyter notebook (.ipynb).

Every matrix constructor that takes a System rebuilds the same intermediates from scratch — the Ybus, the incidence matrix A (IncidenceMatrix), and the susceptance-weighted BA (BA_Matrix). This guide shows how to compute the shared pieces once and feed them to the constructors that accept pre-built matrices.

using PowerNetworkMatrices
import PowerSystemCaseBuilder

Load an example test system with PowerSystemCaseBuilder.build_system:

sys = PowerSystemCaseBuilder.build_system(
    PowerSystemCaseBuilder.PSITestSystems,
    "c_sys5",
);
[ Info: Loaded time series from storage file existing=/home/runner/.julia/packages/PowerSystemCaseBuilder/zuCHR/data/serialized_system/8584b9e729c8aa68ee5405660c6258cde1f67ed3b68f114823707912e9a0d16c/c_sys5_time_series_storage.h5 new=/tmp/jl_32fi7I compression=InfrastructureSystems.CompressionSettings(false, InfrastructureSystems.CompressionTypesModule.CompressionTypes.DEFLATE = 1, 3, true)

Build the shared intermediates once

The construction dependency chain is

Ybus → IncidenceMatrix, BA_Matrix → ABA_Matrix / PTDF, and PTDF → LODF.

The Ybus is the expensive shared root. Build it — and the IncidenceMatrix and BA_Matrix derived from it — exactly once:

ybus = Ybus(sys)
A = IncidenceMatrix(ybus)
BA = BA_Matrix(ybus)
DC_BA_Matrix
    Dimension 1, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
    Dimension 2, [1, 2, 3, 4, 5]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (5, 6):
  35.5872   156.25     ⋅       ⋅         ⋅     32.8947
 -35.5872      ⋅       ⋅     92.5926     ⋅       ⋅ 
    ⋅          ⋅     33.67  -92.5926     ⋅       ⋅ 
    ⋅          ⋅    -33.67     ⋅       33.67  -32.8947
    ⋅      -156.25     ⋅       ⋅      -33.67     ⋅

Reuse them across constructors

PTDF accepts the incidence and BA matrices directly, skipping its own Ybus build:

ptdf = PTDF(A, BA)
DC_PTDF_Matrix
    Dimension 1, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
    Dimension 2, [1, 2, 3, 4, 5]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (5, 6):
  0.193917   0.368495  0.193917   0.193917  -0.368495  0.437588
 -0.475895   0.217552  0.524105   0.524105  -0.217552  0.258343
 -0.348989   0.159538  0.651011  -0.348989  -0.159538  0.189451
  0.0        0.0       0.0        0.0        0.0       0.0
  0.159538  -0.519548  0.159538   0.159538  -0.480452  0.36001

LODF can be built straight from a PTDF you already have, reusing that work too — no second factorization of the network:

lodf = LODF(A, ptdf)
LODF{Tuple{Vector{Tuple{Int64, Int64}}, Vector{Tuple{Int64, Int64}}}, Tuple{Dict{Tuple{Int64, Int64}, Int64}, Dict{Tuple{Int64, Int64}, Int64}}, Matrix{Float64}}
    Dimension 1, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
    Dimension 2, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (6, 6):
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
  0.307071  -1.0        0.307071   0.307071   1.0        0.692929
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
 -0.307071   1.0       -0.307071  -0.307071  -1.0       -0.692929
  0.344795   0.655205   0.344795   0.344795  -0.655205  -1.0

Alternatively, the factorized ABA_Matrix route builds LODF from the same A and BA. The factorize keyword is what pre-factorizes it; is_factorized reports the result. All three inputs must share the same network reduction — which they do here, because they all descend from one ybus:

aba = ABA_Matrix(ybus; factorize = true)
is_factorized(aba)
true

That factorized ABA_Matrix, the incidence matrix, and the BA matrix are all LODF needs:

lodf_via_aba = LODF(A, aba, BA)
LODF{Tuple{Vector{Tuple{Int64, Int64}}, Vector{Tuple{Int64, Int64}}}, Tuple{Dict{Tuple{Int64, Int64}, Int64}, Dict{Tuple{Int64, Int64}, Int64}}, Matrix{Float64}}
    Dimension 1, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
    Dimension 2, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
Note!! The data shown below corresponds to the transposed matrix.
And data with size (6, 6):
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
  0.307071  -1.0        0.307071   0.307071   1.0        0.692929
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
 -1.0        0.457143  -1.0       -1.0       -0.457143   0.542857
 -0.307071   1.0       -0.307071  -0.307071  -1.0       -0.692929
  0.344795   0.655205   0.344795   0.344795  -0.655205  -1.0

Virtual matrices likewise accept a pre-built Ybus, so VirtualPTDF and its siblings reuse the same root:

vptdf = VirtualPTDF(ybus)
DC_vPTDF_Matrix{PowerNetworkMatrices.KLUWrapper.KLULinSolveCache{Float64, Int64}}
    Dimension 1, [(1, 2), (1, 5), (3, 4), (2, 3), (4, 5), (1, 4)]
    Dimension 2, [1, 2, 3, 4, 5]
And data with size (5, 6)
Keep reductions consistent

Constructors that combine pre-built matrices (e.g. LODF(A, ABA, BA)) require every input to have been built with the same network_reductions. Because they all derive from a single Ybus here, they are automatically consistent. Pass network_reductions once, to the Ybus call, and everything downstream inherits it. See the NetworkReduction docstring for the keyword and its rules.

See also