From 59387ee2a3e9dcba5563a9259f821c8e17094acc Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Tue, 16 Jun 2026 04:01:45 -0400 Subject: [PATCH 1/2] Expand and correct math details in the power flow docs Verify each added formula against the sensitivity implementations. - DC power flow: write out the switching flow sensitivity (direct and indirect parts) and expand the angle sensitivity term into explicit form - DC OPF: gate the angle difference limits by sw to match the model, and add the explicit per-branch KKT switching derivative blocks - AC power flow: add an Admittance and Topology Sensitivity section and show the sw gate in the bus admittance matrix - AC OPF: note the implicit function framing of the KKT differentiation - Cite Talkington et al. (arXiv:2510.17071) in the AC docs --- docs/src/math/ac-opf.md | 8 +++++ docs/src/math/ac-power-flow.md | 56 ++++++++++++++++++++++++++++++++-- docs/src/math/dc-opf.md | 33 +++++++++++++++++--- docs/src/math/dc-power-flow.md | 33 +++++++++++++++++--- 4 files changed, 118 insertions(+), 12 deletions(-) diff --git a/docs/src/math/ac-opf.md b/docs/src/math/ac-opf.md index b4ecedf..354e67f 100644 --- a/docs/src/math/ac-opf.md +++ b/docs/src/math/ac-opf.md @@ -105,6 +105,8 @@ with total dimension ``6n + 12m + 6k + n_{\text{ref}}``, where ``n`` is the numb ## Implicit Differentiation +Optimal AC OPF solutions are differentiated by applying the implicit function theorem to the KKT system, which is the constrained optimization counterpart of implicitly differentiating the power flow equations themselves. The same theory guarantees that the derivative exists under mild regularity (a nonsingular KKT Jacobian at a strictly complementary solution) for generic radial or meshed networks. + ### KKT Jacobian The KKT Jacobian ``\partial K / \partial z`` is assembled analytically as a @@ -185,3 +187,9 @@ lmps = calc_sensitivity(ac_prob, :lmp, :d) # dLMP/dd ## Solver The AC OPF uses [Ipopt](https://github.com/coin-or/Ipopt) as the default nonlinear programming solver, accessed via JuMP. + +## References + +The implicit differentiation framework underlying these sensitivities, which applies the implicit function theorem to power flow solutions and their optimization counterparts under admittance and topology changes, follows: + +- S. Talkington, D. Turizo, S. A. Dorado-Rojas, R. K. Gupta & D. K. Molzahn, ["Differentiating Through Power Flow Solutions for Admittance and Topology Control,"](https://arxiv.org/abs/2510.17071) 2025. diff --git a/docs/src/math/ac-power-flow.md b/docs/src/math/ac-power-flow.md index e19ae7a..27cc3a2 100644 --- a/docs/src/math/ac-power-flow.md +++ b/docs/src/math/ac-power-flow.md @@ -2,13 +2,13 @@ ## Admittance Matrix -The bus admittance matrix in rectangular form is: +The bus admittance matrix is assembled from the branch admittances and the incidence matrix: ```math -Y = A^\top \operatorname{diag}(g + jb) \, A + \operatorname{diag}(g_{\text{sh}} + jb_{\text{sh}}) +Y = A^\top \operatorname{diag}\bigl(\mathrm{sw} \circ (g + jb)\bigr) \, A + \operatorname{diag}(g_{\text{sh}} + jb_{\text{sh}}) ``` -where ``g`` and ``b`` are branch conductances and susceptances, and ``g_{\text{sh}}``, ``b_{\text{sh}}`` are shunt admittances. +where ``A`` is the ``m \times n`` incidence matrix (row ``e``, written ``a_e^\top``, has ``+1`` at branch ``e``'s from-bus and ``-1`` at its to-bus), ``g`` and ``b`` are branch conductances and susceptances, ``\mathrm{sw}_e \in [0,1]`` is the switching state that scales branch ``e``'s admittance (``\mathrm{sw}_e = 0`` removes the branch), and ``g_{\text{sh}}``, ``b_{\text{sh}}`` are shunt admittances. An off-nominal transformer replaces each branch's rank-1 contribution ``a_e a_e^\top`` by a ``2 \times 2`` primitive built from the complex tap ratio. ## Power Injection Equations @@ -98,3 +98,53 @@ The power flow formulation uses power injections ``(p, q)`` as native parameters ``` These transforms are applied automatically by the unified interface when using `:d` or `:qd` as the parameter symbol. + +## Admittance and Topology Sensitivity + +The sections above differentiate the power flow with respect to power injections. The same implicit function machinery differentiates it with respect to the **network admittance parameters**: branch conductance ``g_e`` and susceptance ``b_e``. This is the basis for topology and admittance control, predicting how the voltage state moves when a branch's admittance changes, or when a line is switched, *without* solving the power flow equations again. + +At a solved operating point the power balance residual ``F(V, \beta) = 0`` depends on a parameter ``\beta \in \{g_e, b_e\}`` only through the admittance matrix ``Y = A^\top \operatorname{diag}(\mathrm{sw} \circ (g + jb)) \, A``. Implicit differentiation reuses the same Jacobian ``J`` used for voltage sensitivity to power injections (same factorization as above, with a different right side): + +```math +\frac{\partial V}{\partial \beta} = -J^{-1} \frac{\partial F}{\partial \beta}, \qquad +\frac{\partial F}{\partial \beta} = \begin{bmatrix} \partial P / \partial \beta \\ \partial Q / \partial \beta \end{bmatrix} +``` + +Each branch enters ``Y`` as a rank-1 update, so the parameter derivatives are sparse: + +```math +\frac{\partial Y}{\partial g_e} = \mathrm{sw}_e \, a_e a_e^\top, \qquad +\frac{\partial Y}{\partial b_e} = j \, \mathrm{sw}_e \, a_e a_e^\top +``` + +(shown for unit tap ratio; an off-nominal complex tap replaces the outer product by the branch's ``2 \times 2`` primitive). Let ``\Delta V = A V`` be the vector of edge voltage drops (so ``\Delta V_e = V_f - V_t`` across branch ``e``). Still at unit tap, define the matrix ``M`` with entries + +```math +M_{i,e} = \mathrm{sw}_e \, \bar{V}_i \,(a_e)_i \, \Delta V_e , +``` + +which is nonzero only at branch ``e``'s two endpoints. Using ``P = \operatorname{Re}(\bar V \circ YV)`` and ``Q = -\operatorname{Im}(\bar V \circ YV)``, the injection derivatives are: + +```math +\frac{\partial P}{\partial g} = \operatorname{Re}(M), \quad +\frac{\partial Q}{\partial g} = -\operatorname{Im}(M); \qquad +\frac{\partial P}{\partial b} = -\operatorname{Im}(M), \quad +\frac{\partial Q}{\partial b} = -\operatorname{Re}(M). +``` + +The conductance and susceptance columns share the same ``M`` and differ only in which part is taken: the factor ``j`` in ``\partial Y / \partial b_e`` rotates the real and imaginary parts into each other. Solving for ``\partial V / \partial \beta`` and projecting onto magnitude and angle exactly as in the voltage sensitivity case, + +```math +\frac{\partial |V_i|}{\partial \beta} = \operatorname{Re}\!\left(\frac{\partial V_i}{\partial \beta} \cdot \frac{\bar V_i}{|V_i|}\right), \qquad +\frac{\partial \theta_i}{\partial \beta} = \operatorname{Im}\!\left(\frac{\partial V_i}{\partial \beta} \cdot \frac{\bar V_i}{|V_i|^2}\right), +``` + +gives the admittance sensitivities `calc_sensitivity(state, :vm, :g)`, `calc_sensitivity(state, :va, :b)`, and so on. Branch current and flow sensitivities to ``g`` and ``b`` follow by the same chain rule used in the Current Sensitivity and Branch Flow Sensitivity sections above, with ``\partial V / \partial \beta`` substituted for ``\partial V / \partial p_k``. + +In AC OPF the analogous topology parameter is the switching state ``\mathrm{sw}``, which scales the whole branch admittance primitive (``\partial Y / \partial \mathrm{sw}_e = (g_e + jb_e)\, a_e a_e^\top`` at unit tap). Its sensitivities come from implicitly differentiating the KKT system (see [AC Optimal Power Flow](@ref)). + +## References + +The implicit differentiation framework for the power flow equations, including the guarantee that the derivative exists for generic radial or meshed networks and its application to voltage, current, and power flow sensitivities under admittance and topology changes, follows: + +- S. Talkington, D. Turizo, S. A. Dorado-Rojas, R. K. Gupta & D. K. Molzahn, ["Differentiating Through Power Flow Solutions for Admittance and Topology Control,"](https://arxiv.org/abs/2510.17071) 2025. diff --git a/docs/src/math/dc-opf.md b/docs/src/math/dc-opf.md index 084e484..1e0a093 100644 --- a/docs/src/math/dc-opf.md +++ b/docs/src/math/dc-opf.md @@ -17,7 +17,7 @@ f &= W A \theta & (\nu_{\text{flow}}) \\ -f_{\max} \leq f &\leq f_{\max} & (\lambda_{\text{lb}}, \lambda_{\text{ub}}) \\ g_{\min} \leq g &\leq g_{\max} & (\rho_{\text{lb}}, \rho_{\text{ub}}) \\ 0 \leq \text{psh} &\leq d_+ & (\mu_{\text{lb}}, \mu_{\text{ub}}) \\ -\alpha_{\min} \leq A\theta &\leq \alpha_{\max} & (\gamma_{\text{lb}}, \gamma_{\text{ub}}) \\ +\mathrm{sw} \circ \alpha_{\min} \leq \mathrm{sw} \circ A\theta &\leq \mathrm{sw} \circ \alpha_{\max} & (\gamma_{\text{lb}}, \gamma_{\text{ub}}) \\ \theta_{\text{ref}} &= 0 & (\eta_{\text{ref}}) \end{aligned} ``` @@ -32,6 +32,7 @@ where: - ``c_{\text{shed}}`` is the load shedding cost vector - ``d_+ = \max(d, 0)`` is the curtailable portion of signed net demand; negative net demand remains in power balance as an injection - ``\tau`` is a small regularization parameter for numerical conditioning +- the angle difference limits are gated by ``\mathrm{sw}`` so an open branch (``\mathrm{sw}_e = 0``) imposes no limit; the factor cancels for a branch in service (``\mathrm{sw}_e = 1``) ## KKT System for Implicit Differentiation @@ -57,12 +58,12 @@ with total dimension ``5n + 6m + 3k + 1``. The KKT residual ``K(z, p)`` consists of: -1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + e_{\text{ref}} \eta_{\text{ref}} + A^\top (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0`` +1. **Stationarity w.r.t. ``\theta``**: ``B^\top \nu_{\text{bal}} + (WA)^\top \nu_{\text{flow}} + e_{\text{ref}} \eta_{\text{ref}} + A^\top \operatorname{diag}(\mathrm{sw}) (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0`` 2. **Stationarity w.r.t. ``g``**: ``2 C_q g + c_l - G_{\text{inc}}^\top \nu_{\text{bal}} - \rho_{\text{lb}} + \rho_{\text{ub}} = 0`` 3. **Stationarity w.r.t. ``f``**: ``\tau^2 f - \nu_{\text{flow}} - \lambda_{\text{lb}} + \lambda_{\text{ub}} = 0`` 4. **Stationarity w.r.t. psh**: ``c_{\text{shed}} - \nu_{\text{bal}} - \mu_{\text{lb}} + \mu_{\text{ub}} = 0`` 5. **Complementary slackness (flow bounds)**: ``\lambda_{\text{lb}} \circ (f + f_{\max}) = 0``, ``\lambda_{\text{ub}} \circ (f_{\max} - f) = 0`` -5b. **Complementary slackness (angle differences)**: ``\gamma_{\text{lb}} \circ (A\theta - \alpha_{\min}) = 0``, ``\gamma_{\text{ub}} \circ (\alpha_{\max} - A\theta) = 0`` +5b. **Complementary slackness (angle differences)**: ``\gamma_{\text{lb}} \circ \mathrm{sw} \circ (A\theta - \alpha_{\min}) = 0``, ``\gamma_{\text{ub}} \circ \mathrm{sw} \circ (\alpha_{\max} - A\theta) = 0`` 5c. **Complementary slackness (generation/shedding bounds)**: ``\rho \circ (\cdot) = 0``, ``\mu \circ (\cdot) = 0`` 6. **Primal feasibility**: ``G_{\text{inc}} g + \text{psh} - d - B\theta = 0`` 7. **Flow definition**: ``f - WA\theta = 0`` @@ -94,14 +95,36 @@ the collapsed bound ``0 \leq \text{psh} \leq 0``. ### Switching (``\mathrm{sw}``) -Switching affects the Laplacian ``B``, weight matrix ``W``, and flow definition through: +Switching enters the Laplacian ``B``, the weight matrix ``W``, and the angle difference limits gated by ``\mathrm{sw}``, with elementary perturbations: ```math \frac{\partial B}{\partial \mathrm{sw}_e} = -b_e \, a_e a_e^\top, \qquad \frac{\partial W}{\partial \mathrm{sw}_e} = -b_e \, e_e e_e^\top ``` -This propagates into the stationarity, power balance, and flow definition blocks of the KKT system. +where ``a_e^\top`` is row ``e`` of ``A`` (so ``a_e^\top \theta = (A\theta)_e``) and ``e_e`` is the unit vector for branch ``e``. These yield the nonzero parameter derivative blocks of ``\partial K / \partial \mathrm{sw}_e``: + +```math +\begin{aligned} +\frac{\partial K_{\nu_{\text{bal}}}}{\partial \mathrm{sw}_e} + &= -\frac{\partial B}{\partial \mathrm{sw}_e}\,\theta = b_e\,(a_e^\top \theta)\, a_e + && \text{(power balance)} \\ +\frac{\partial K_{\nu_{\text{flow}}}}{\partial \mathrm{sw}_e} + &= -\frac{\partial W}{\partial \mathrm{sw}_e}\,A\theta = b_e\,(a_e^\top \theta)\, e_e + && \text{(flow definition)} \\ +\frac{\partial K_{\theta}}{\partial \mathrm{sw}_e} + &= -b_e\,\bigl(a_e^\top \nu_{\text{bal}} + (\nu_{\text{flow}})_e\bigr)\, a_e + + (\gamma_{\text{ub},e} - \gamma_{\text{lb},e})\, a_e + && (\theta\text{ stationarity}) \\ +\frac{\partial K_{\gamma_{\text{lb}}}}{\partial \mathrm{sw}_e} + &= \gamma_{\text{lb},e}\,\bigl((A\theta)_e - \alpha_{\min,e}\bigr)\, e_e, \quad +\frac{\partial K_{\gamma_{\text{ub}}}}{\partial \mathrm{sw}_e} + = \gamma_{\text{ub},e}\,\bigl(\alpha_{\max,e} - (A\theta)_e\bigr)\, e_e + && \text{(angle limits)} +\end{aligned} +``` + +Each block is rank-1 in the incidence row ``a_e`` or supported on branch ``e`` alone, so every column of ``\partial K / \partial \mathrm{sw}`` has only a handful of nonzeros. The stationarity block combines the Laplacian term (``\partial B / \partial \mathrm{sw}_e``), the flow coupling term (``\partial W / \partial \mathrm{sw}_e``), and the gated angle limit term (``\partial \operatorname{diag}(\mathrm{sw}) / \partial \mathrm{sw}_e = e_e e_e^\top``). ### Cost Coefficients (``c_q``, ``c_l``) diff --git a/docs/src/math/dc-power-flow.md b/docs/src/math/dc-power-flow.md index 62c63ef..f950d13 100644 --- a/docs/src/math/dc-power-flow.md +++ b/docs/src/math/dc-power-flow.md @@ -21,7 +21,7 @@ The susceptance-weighted Laplacian is: B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) \, A ``` -where ``A`` is the ``m \times n`` incidence matrix and ``b`` stores the imaginary part of the inverse impedance (``b_e = \operatorname{Im}(1/z_e) < 0`` for inductive branches, so ``-b > 0``). +where ``A`` is the ``m \times n`` incidence matrix (row ``e``, written ``a_e^\top``, has ``+1`` at branch ``e``'s from-bus and ``-1`` at its to-bus, so ``a_e^\top \theta = (A\theta)_e`` is the angle difference across the branch), and ``b`` stores the imaginary part of the inverse impedance (``b_e = \operatorname{Im}(1/z_e) < 0`` for inductive branches, so ``-b > 0``). ## Switching Sensitivity @@ -32,7 +32,7 @@ Switching sensitivity follows from matrix perturbation theory. For a branch ``e` = -B_r^{-1} \frac{\partial B_r}{\partial \mathrm{sw}_e} \, \theta_r ``` -where the perturbation is a rank-1 update from the incidence column of branch ``e`` restricted to non-reference buses: +where the perturbation is a rank-1 update from row ``a_e^\top`` of branch ``e`` restricted to the non-reference buses (denoted ``a_{e,r}``): ```math \frac{\partial B_r}{\partial \mathrm{sw}_e} = -b_e \, a_{e,r} \, a_{e,r}^\top @@ -40,12 +40,37 @@ where the perturbation is a rank-1 update from the incidence column of branch `` ### Flow Sensitivity to Switching -Branch flows are ``f = W A \theta`` where ``W = \operatorname{diag}(-b \circ \mathrm{sw})``. The flow sensitivity has both indirect (via angle changes) and direct (via the switching coefficient) contributions: +Branch flows are ``f = W A \theta`` with ``W = \operatorname{diag}(-b \circ \mathrm{sw})``, so branch ``e`` carries ``f_e = -b_e \, \mathrm{sw}_e \,(a_e^\top \theta)``. Both the weight ``W`` and the angles ``\theta`` depend on ``\mathrm{sw}_e``, so the product rule splits the sensitivity into an **indirect** part (every branch feels the change in angles) and a **direct** part (only branch ``e``'s own weight changes): ```math -\frac{\partial f}{\partial \mathrm{sw}_e} = W A \frac{\partial \theta}{\partial \mathrm{sw}_e} + \text{direct effect on edge } e +\frac{\partial f}{\partial \mathrm{sw}_e} + = \underbrace{W A \frac{\partial \theta}{\partial \mathrm{sw}_e}}_{\text{indirect (all branches)}} + + \underbrace{\frac{\partial W}{\partial \mathrm{sw}_e} A \theta}_{\text{direct (branch } e \text{ only)}} + = W A \frac{\partial \theta}{\partial \mathrm{sw}_e} \;-\; b_e \,(a_e^\top \theta)\, e_e ``` +The direct term uses ``\partial W / \partial \mathrm{sw}_e = -b_e \, e_e e_e^\top``. Only the diagonal entry of ``W`` for branch ``e`` depends on ``\mathrm{sw}_e``, so it adds ``-b_e \,(a_e^\top \theta) = -b_e \,(A\theta)_e`` to the flow on branch ``e`` alone, where ``e_e \in \mathbb{R}^m`` is the unit vector for branch ``e``. + +Substituting the angle sensitivity ``\partial \theta / \partial \mathrm{sw}_e`` we get + +```math +\frac{\partial \theta_r}{\partial \mathrm{sw}_e} + = -B_r^{-1} \frac{\partial B_r}{\partial \mathrm{sw}_e}\, \theta_r + = b_e \,(a_{e,r}^\top \theta_r)\, B_r^{-1} a_{e,r} +``` + +(embedded into the full bus vector with the reference entry held at zero) which makes the flow sensitivity fully explicit: + +```math +\frac{\partial f}{\partial \mathrm{sw}_e} + = b_e \,(a_e^\top \theta) \Big( + \underbrace{W A_r\, B_r^{-1} a_{e,r}}_{\text{indirect (all branches)}} + \;-\; \underbrace{e_e}_{\text{direct (branch } e \text{ only)}} + \Big) +``` + +where ``A_r`` is the incidence matrix with its reference bus column removed (so ``A \,\partial\theta/\partial \mathrm{sw}_e = A_r \,\partial\theta_r/\partial \mathrm{sw}_e``), and we used ``a_{e,r}^\top \theta_r = a_e^\top \theta`` since ``\theta_{\text{ref}} = 0``. The shared scalar ``b_e \,(a_e^\top \theta) = b_e \,(A\theta)_e`` is exactly the negative of branch ``e``'s flow when it is in service (``f_e = -b_e (A\theta)_e`` at ``\mathrm{sw}_e = 1``), so the whole flow sensitivity scales with how heavily branch ``e`` is loaded. + ## Demand Sensitivity Since ``p = g - d`` and generation is fixed, ``\partial p / \partial d = -I``. The angle sensitivity to demand is: From 8d96f3becb8b342f2e5eb8cfa1d80f732a4ccceb Mon Sep 17 00:00:00 2001 From: ckhanpour3 Date: Tue, 16 Jun 2026 19:21:29 -0400 Subject: [PATCH 2/2] Add the direct admittance term to AC current and flow sensitivities Branch current and flow sensitivities to g and b include a direct term on the perturbed branch (its own admittance acting on its voltage drop), in addition to the term through the voltage change. This matches the implementation (the direct topology current contribution) and the paper's line current and flow sensitivity equations. --- docs/src/math/ac-power-flow.md | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/docs/src/math/ac-power-flow.md b/docs/src/math/ac-power-flow.md index 27cc3a2..441edf2 100644 --- a/docs/src/math/ac-power-flow.md +++ b/docs/src/math/ac-power-flow.md @@ -139,7 +139,15 @@ The conductance and susceptance columns share the same ``M`` and differ only in \frac{\partial \theta_i}{\partial \beta} = \operatorname{Im}\!\left(\frac{\partial V_i}{\partial \beta} \cdot \frac{\bar V_i}{|V_i|^2}\right), ``` -gives the admittance sensitivities `calc_sensitivity(state, :vm, :g)`, `calc_sensitivity(state, :va, :b)`, and so on. Branch current and flow sensitivities to ``g`` and ``b`` follow by the same chain rule used in the Current Sensitivity and Branch Flow Sensitivity sections above, with ``\partial V / \partial \beta`` substituted for ``\partial V / \partial p_k``. +gives the admittance sensitivities `calc_sensitivity(state, :vm, :g)`, `calc_sensitivity(state, :va, :b)`, and so on. Branch current and flow sensitivities to ``g`` and ``b`` carry a **direct** term in addition to this voltage chain rule, because the perturbed branch's own admittance enters its flow. For a shuntless branch ``I_\ell = (g_\ell + j b_\ell)(V_f - V_t)``, + +```math +\frac{\partial I_\ell}{\partial \beta} + = \underbrace{\frac{\partial (g_\ell + j b_\ell)}{\partial \beta}\,(V_f - V_t)}_{\text{direct (branch } \ell)} + + \underbrace{(g_\ell + j b_\ell)\left(\frac{\partial V_f}{\partial \beta} - \frac{\partial V_t}{\partial \beta}\right)}_{\text{indirect (through the voltage change)}}, +``` + +which mirrors the direct and indirect split of the DC switching sensitivity. The direct term is nonzero only when ``\beta`` is branch ``\ell``'s own ``g`` or ``b`` (then ``\partial (g_\ell + j b_\ell)/\partial g_\ell = 1`` or ``\partial (g_\ell + j b_\ell)/\partial b_\ell = j``); for every other branch only the voltage term remains. Power flow sensitivities follow from ``S_\ell = V_f \bar I_\ell`` by the product rule, and the implementation adds this direct contribution for `:im` and `:f` queries. In AC OPF the analogous topology parameter is the switching state ``\mathrm{sw}``, which scales the whole branch admittance primitive (``\partial Y / \partial \mathrm{sw}_e = (g_e + jb_e)\, a_e a_e^\top`` at unit tap). Its sensitivities come from implicitly differentiating the KKT system (see [AC Optimal Power Flow](@ref)).