Different Jacobian Algorithms

[luke-kiernan]

Current data layout

Our Jacobian matrix is stored as a SparseMatrixCSC, of dimensions 2n-by-2n, where n is the number of buses. Each 2x2 block is structurally nonzero if and only if there's a branch between that pair of buses. If any one entry in a 2x2 block is structurally nonzero, then all 4 are structurally nonzero. This means we're storing a few extra zeros, but the tradeoff is that our sparse structure doesn't depend on the bus types (so PV/PQ switching isn't a problem).

When writing code, I frequently find myself referring to the following chart:

Bus Type	Input variables in x	Power Balance Equations in F
PQ	V_m, V_θ	ΔP, ΔQ
PV	Q_inj, V_θ	ΔP, ΔQ
REF	P_inj, Q_inj	ΔP, ΔQ

Rows in J correspond to entries in F, columns to entries in x, so the above also describes the structure of J. If that isn't clear, I wrote down a sample Jacobian matrix for a 3-bus system in the docstring for _create_jacobian_matrix_structure.

Iterating over Neighbors

This is what we're currently using. There's been some refactoring for performance reasons, but the core algorithm amounts to

for bus_from in 1:num_buses
    for bus_to in neighbors(bus_from)
        # set entries in 2x2 block [2*bus_from - 1, 2*bus_from] x [2*bus_to - 1, 2*bus_to]
    end
end

Note that this works across-then-down, when down-then-across cooperates much better with CSC format. The no-frills loop can be found here. Main performance improvements since then:

1. Multiple dispatch on Val of the bus type of the bus_to, for compilation reasons.
2. The diagonal entries (from and to bus identical) are a sum() over the other neighbors. Since we're iterating over the neighbors anyway, total it up as we go.
3. Constructing Val objects at runtime is slow: instead of having Val(bus_type), spell out the 3 possibilities in an if-elseif so that it's Val(constant).

2-pass algorithm

This is described in Algorithm 1 of the paper "Bus Admittance Matrix Revisited: Performance Challenges on Modern Computers": it iterates directly over the sparse structure. However, it stores the angle derivatives and voltage magnitude derivatives in 2 separate CSC matrices, dS_Θ and dS_V. The sparse structure of dS_Θ and dS_V is the same as that of the Ybus matrix: their initialization step of "copy the contents of Ybus to dS_Θ" is linear memory access, hence very fast. With our data layout, we don't have that feature.

Nevertheless, I adapted their general algorithm on the lk/2-pass-jacbian branch. Profiling the Jacobian calls alone, I get that it's about 2x faster than our current implementation. Running the TrustRegion method on a 10k bus system, it's a 15% reduction in total runtime, albeit with a 15% increase in memory usage. I've eliminated all unnecessary allocations in my implementation of their algorithm, but I haven't experimented with further performance tweaks.

Matrix algorithm

Worth mentioning for readibility and mathematical elegance: the matpower docs have a nice derivation using chain rule. Implemented a while back: poor performance, like 10x slower than the current method. However, it's still used in testing and can be found in test/test_utils/legacy_pf.jl.

[luke-kiernan]

I made some revisions and rebased the lk/2-pass-jacobian to a more recent branch of PF. Two things worth mentioning:

1. My implementation isn't actually 2-pass. Where the paper does element-wise operations in a straight loop, I'm broadcasting operations over the whole range. Part of this is for readability reasons; part of this is that I also expect this to be faster (vectorization). However, I haven't written a straight-loop version and done a performance comparison.
2. 80% of the complexity in the code comes from the fact that Jv stores things using real entries, while the Y bus matrix storing things using complex entries. As a result, where Y.rowval has a single index (k if in kth row), Jv.rowval will have a pair of indices (2k-1 and 2k) for its pair of entries. The right abstraction or interface could make this all so much easier. For example, if reinterpret worked on sparse arrays as I suggest here, instead of

J_nonzeros = reinterpret(ComplexF64, SparseArrays.nonzeros(Jv))
for col in 1:num_buses
    Yb_val_range = nzrange(Yb, col)
    J_Vm_start = div(Jv.colptr[2 * col - 1] + 1, 2)
    J_Vm_val_range = Yb_val_range .+ (J_Vm_start - Yb_val_range.start)
    J_Va_start = div(Jv.colptr[2 * col] + 1, 2)
    J_Va_val_range = Yb_val_range .+ (J_Va_start - Yb_val_range.start)
    # do something with the ranges
end

I could write

J_complex = reinterpret(ComplexF64, Jv)
for col in 1:num_buses
    Yb_val_range = nzrange(Yb, col)
    J_Vm_val_range = nzrange(J_complex, 2*col-1)
    J_Va_val_range = nzrange(J_complex, 2*col)
    # do something with the ranges
end

and the reindexing would be handled behind the scenes. [Explicitly: the data width doubled, so J_complex.rowval equals half of every-other element of Jv.rowval.]