AC Sensitivity Analysis

Motivation

A power flow answers what is the state of the network? Most operational questions are one step past that: what happens if something changes?

• This generator is about to ramp 50 MW. Which line picks it up?
• This corridor is at 98% of its limit. Which injection relieves it fastest, and by how much per MW?
• The voltage at this bus is sagging. Would a tap change help, or does it need reactive support?

Each is a partial derivative of the solved operating point with respect to an input. You could answer any of them by re-solving — perturb, solve again, subtract — and for one question that is perfectly reasonable. It stops being reasonable at a hundred questions, and it is unnecessary: the derivative is available in closed form from a factorization the solve already computed.

The formulation

A converged power flow satisfies

\[ g(x, p) = s_{calc}(x) - s_{spec}(p) = 0 \]

where \(x = [\theta_{non\text{-}slack};\ |V|_{PQ}]\) is exactly the unknown vector Newton iterates on, and \(p\) is whatever input we are differentiating with respect to. Differentiating the identity and rearranging:

\[ \frac{dx}{dp} = -J^{-1} \frac{\partial g}{\partial p} \]

with \(J = \partial g / \partial x\) — the same Jacobian Newton already builds, filled at the converged state. Nothing new is assembled.

For an injection the right-hand side is especially simple. A specified injection enters exactly one equation, so \(\partial g/\partial p = -e_i\) and

\[ \frac{dx}{dp_i} = J^{-1} e_i \]

— a single triangular solve with a unit right-hand side. That is the entire method for the common case.

Any quantity of interest \(f\) then follows by the chain rule:

\[ \frac{df}{dp} = \left(\frac{\partial f}{\partial x}\right)^{\mathsf T} \frac{dx}{dp} + \left.\frac{\partial f}{\partial p}\right|_{x} \]

For a branch flow, \(\partial f/\partial x\) is already computed by branch_flow::terminal_flow_derivs — written originally for the state estimator, reused here unchanged. The second term is zero unless \(p\) is the branch's own tap; see below.

Two directions

The expression above can be bracketed either way, and the two groupings cost very differently:

Forward is "this generator ramps — what happens to all 5,000 branches?" Adjoint is "this line is overloaded — which of 500 injections relieves it?" Neither is derivable cheaply from the other, so both are offered: state_response/branch_response and function_row.

Both run against one factorization, taken once in AcSensitivity::new. sparse::RealFactorization::solve_transpose reuses the same LU for the transposed direction rather than factorizing \(J^{\mathsf T}\) separately — an LU = PAQ decomposition solves the transposed system by running the same triangular factors in the opposite order.

This mirrors the ptdf_column/ptdf_row pair in DC power flow, which is the same duality on the DC side.

Variables

The zeros are answers, not gaps. A slack bus's injection is not an input at all — it is what the solve determines — so perturbing it moves nothing, exactly as a DC PTDF column is zero at its reference bus. Likewise reactive power at a PV bus is an output of the voltage controller.

Branch indices are the crate's flat order: lines first, then transformers.

Taps and phase shifters

A tap is a different shape of variable from an injection, and gets two things wrong if treated casually.

The sign flips. An injection enters through \(s_{spec}\), so \(\partial g/\partial p = -\partial s_{spec}/\partial p\). A tap changes the branch's admittance and therefore the calculated injections at both its ends, so \(\partial g/\partial p = +\partial s_{calc}/\partial p\). The response picks up the opposite sign.

The direct term. A tapped branch's own flow depends on the tap directly, not only through the state — that \(\partial f/\partial p|_{x}\) term above. It is zero for every branch except the tapped one, where it is large.

This is not a rounding-level concern. On distribution-case, dropping it changes \(dP_8/dk\) from \(-0.276\) to \(+0.286\): a sign flip, on the branch whose tap is being adjusted, which is the branch a user is most likely to read. Every other branch stays correct, which is exactly what makes the omission dangerous — the answer looks entirely plausible.

The admittance derivatives themselves are unusually tidy. With \(tap = k,e^{j\alpha}\), and given that \(y_{ff}\) always carries \(1/k^2\), \(y_{ft}\) and \(y_{tf}\) always carry \(1/k\), and \(y_{tt}\) never depends on the tap at all:

\[ \frac{\partial y_{ff}}{\partial k} = \frac{-2,y_{ff}}{k}, \quad \frac{\partial y_{ft}}{\partial k} = \frac{-y_{ft}}{k}, \quad \frac{\partial y_{tf}}{\partial k} = \frac{-y_{tf}}{k}, \quad \frac{\partial y_{tt}}{\partial k} = 0 \]

\[ \frac{\partial y_{ft}}{\partial \alpha} = j,y_{ft}, \quad \frac{\partial y_{tf}}{\partial \alpha} = -j,y_{tf}, \quad \frac{\partial y_{ff}}{\partial \alpha} = \frac{\partial y_{tt}}{\partial \alpha} = 0 \]

Every entry is a scaling of itself, so nothing has to be rebuilt from nameplate data — and the relations hold unchanged for the half-open terminal states, where the affected entries are already zero.

Taps come out of the adjoint at no extra solve, alongside the injection sensitivities. So "which tap relieves this overload?" costs exactly what "which injection relieves it?" does.

How this differs from the DC factors

PTDF and LODF are exact for the model they describe, because the DC model is linear — a PTDF is not an approximation of DC, it is DC. These are different. The AC problem is nonlinear, so these are first-order derivatives: accurate for small perturbations and degrading as the step grows.

What they buy in exchange is everything DC discards — voltage magnitudes, reactive flows, losses. A DC PTDF cannot tell you that an injection will collapse a voltage; an AC sensitivity can.

The usual division of labour: screen thousands of cases in DC, answer precisely about the few that matter in AC.

Islands

None, and none is needed. Every island carries its own slack, so the assembled Jacobian is block-diagonal across islands and non-singular as a whole. A sourceless island contributes no rows at all, because network::mark_unreferenced_islands has already pinned its buses to Slack — so it is excluded rather than making the system singular.

This is simpler than the DC side, where DcSensitivity factorizes each island's \(B\) separately because the reduced \(B\) is built per island.

What is not here

No contingency sensitivities. An outage is a finite change in topology, not an infinitesimal change in an input, so it is not a derivative at all. The DC side handles outages with a Woodbury update (DcSensitivity::outage_flows) that has no equally cheap AC analogue; batch::BatchSolver::solve_contingencies is the AC answer, and it re-solves.

No sensitivities through the outer loops. These differentiate the power-flow equations as posed. If PV→PQ switching would trip at the perturbed point, or a tap controller would respond, the derivative does not know it — it is the derivative of the current configuration.

Validation

tests/ac_sensitivity_test.rs, against an oracle that shares no code with the sensitivity module: perturb an input, run the ordinary Newton solver to convergence, take a central difference of the solved state.

Because AC is nonlinear the finite difference is not exact, so the two constants are chosen together against each other:

• the solve's residual noise is divided by \(2h\), so a smaller step amplifies it;
• the central difference's own truncation error is \(O(h^2)\), so a larger step grows it.

At \(h = 10^{-4}\) and a solve tolerance of \(10^{-10}\) both land near \(10^{-7}\). Tightening the solve further does not help: the test case's residual floors at \(1.2 \times 10^{-11}\) in double precision, so asking for \(10^{-12}\) simply never converges.

Tap sensitivities are checked the same way, by actually changing a Transformer's tap and re-solving, so the whole downstream chain sees it as a real tap change. Forward and adjoint are additionally checked against each other to \(10^{-9}\) — the same scalar bracketed the other way, so any disagreement is a transpose or indexing error rather than a modelling choice.

Using it

From the CLI — forward, then adjoint:

gridoxide sensitivity network.json --dp 2
gridoxide sensitivity network.json --watch 8

--dp/--dq take a bus, --dk/--dalpha a branch, --watch a branch to watch; --terminal from|to picks the end. It refuses to differentiate a power flow that did not converge, rather than printing derivatives of nothing.

From Rust:

#![allow(unused)]
fn main() {
use gridoxide::ac_sensitivity::{AcSensitivity, Function, Variable};
use gridoxide::branch_flow::Terminal;

let sens = AcSensitivity::new(&solved_buses, &ybus, &lines, &transformers)
    .expect("nonsingular at a converged point");

// Forward: bus 2 ramps — what does every branch do?
let flows = sens.branch_response(Variable::ActiveInjection(2), Terminal::From).unwrap();

// Adjoint: branch 8 is loaded — what moves it?
let row = sens
    .function_row(Function::BranchActivePower { branch: 8, terminal: Terminal::From })
    .unwrap();
let (best_bus, _) = row
    .d_active
    .iter()
    .enumerate()
    .max_by(|a, b| a.1.abs().total_cmp(&b.1.abs()))
    .unwrap();
}

From Python:

import gridoxide

s = gridoxide.AcSensitivityModel("network.json")

col = s.column(active_injection=2)   # bus 2 ramps — what responds?
row = s.row(branch=8)                # branch 8 is loaded — what moves it?

# Which tap has the most leverage on branch 8?
best = max(range(s.n_branches), key=lambda b: abs(row.d_transformer_ratio[b]))

AcSensitivityModel is a class rather than a function because there is genuinely something worth keeping between calls: one factorization, arbitrarily many questions.

Second derivatives

ac_sensitivity differentiates the operating point once. injection_hessian differentiates it twice, and exists for one reason: any interior-point method for AC-OPF needs the Hessian of the Lagrangian,

\[ \nabla^2_{xx} L = \nabla^2 f + \sum_i \lambda_i \nabla^2 P_i + \sum_i \mu_i \nabla^2 Q_i \]

whose objective term is trivial and whose constraint terms are these.

Three design choices are worth stating, because each is the opposite of what the power-flow code does and each has a reason.

The weighted sum, never one \(\nabla^2 P_i\) at a time. Each individual Hessian is a sparse \(2n \times 2n\) matrix and there are \(2n\) of them, so materializing them separately would cost \(O(n)\) matrices to build one. Accumulating as we go makes it a single pass over the Y-bus — the same cost as assembling the Jacobian.

The full \(2n\) state, not the reduced one. Newton drops the slack angle and every PV magnitude, because it solves with bus types fixed. AC-OPF makes generator voltages decision variables, so the rows Newton discards are exactly the ones the optimizer needs. Emitting everything and letting the caller select is both more general and easier to check.

No fill beyond the Jacobian. This is not obvious and it is what keeps AC-OPF tractable: \(\nabla^2 P_i\) has no entry coupling two distinct neighbours of bus \(i\), even though both appear in \(P_i\). Each term of \(P_i\) involves exactly one neighbour, so differentiating twice can never bring two together. A second-order method therefore needs no more sparsity budget than a first-order one.

How it is checked

Hessians are hard to test — hand-deriving them twice tends to reproduce the same mistake twice — so the validation is deliberately two-staged, with each stage grounded outside the code it tests.

1. A second Jacobian, checked against the first. injection_jacobian computes the same first derivatives as jacobian::JacobianPattern by a different route: summing over neighbours directly, where the power-flow assembler reads its diagonal blocks out of precomputed injections (H_ii = −Q_i − V_i²B_ii). They share no code, so the mature assembler serves as an independent oracle. On case118_ieee that is over a thousand entries agreeing to 1e-9.
2. The Hessian, central-differenced against that verified Jacobian. Since stage one established the Jacobian independently, this grounds the second derivatives in something outside their own implementation.

Doing it in one step — second-order differences of the injections themselves — would work in principle and be much weaker in practice. A second difference has an error floor of \(O(h^2 + \epsilon/h^2)\), leaving roughly six usable digits where a first difference of an exact Jacobian gives ten.