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.