DC Optimal Power Flow

B-theta Formulation

The DC OPF minimizes generation cost subject to linearized power flow constraints using the susceptance-weighted Laplacian:

\[\min_{g, \theta, f, \text{psh}} \quad g^\top C_q g + c_l^\top g + c_{\text{shed}}^\top \text{psh} + \frac{\tau^2}{2} \|f\|^2\]

subject to:

\[\begin{aligned} G_{\text{inc}} g + \text{psh} - d &= B \theta & (\nu_{\text{bal}}) \\ 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}}) \\ \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{refs}} &= 0 & (\eta_{\text{ref}}) \end{aligned}\]

where:

• $B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) \, A$ is the susceptance-weighted Laplacian
• $W = \operatorname{diag}(-b \circ \mathrm{sw})$ is the branch weight matrix
• $A$ is the $m \times n$ incidence matrix (branches × buses)
• $G_{\text{inc}}$ is the $n \times k$ generator-bus incidence matrix
• $C_q = \operatorname{diag}(c_q)$ contains quadratic cost coefficients
• $c_l$ contains linear cost coefficients
• $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$)
• $\text{refs}$ contains one reference bus per energized island, including isolated buses. The choice is deterministic: the configured $\text{ref_bus}$ is the reference for its island, and every other island uses its lowest sequential bus index

The built OPF model stores one reference constraint for each entry of $\text{refs}$. Consequently, topology changes that can alter the energized island partition must rebuild the model: use update_switching! for switch changes, and rebuild the DCOPFProblem after direct susceptance edits that move a branch across zero.

KKT System for Implicit Differentiation

OPF sensitivities are computed via the implicit function theorem applied to the KKT conditions. At an optimal solution $z^*$, the KKT residual satisfies $K(z^*, p) = 0$ where $z$ collects all primal and dual variables and $p$ is a parameter.

By the implicit function theorem:

\[\frac{dz}{dp} = -\left(\frac{\partial K}{\partial z}\right)^{-1} \frac{\partial K}{\partial p}\]

KKT Variable Ordering

The KKT variable vector $z$ is structured as:

\[z = [\theta, g, f, \text{psh}, \lambda_{\text{lb}}, \lambda_{\text{ub}}, \gamma_{\text{lb}}, \gamma_{\text{ub}}, \rho_{\text{lb}}, \rho_{\text{ub}}, \mu_{\text{lb}}, \mu_{\text{ub}}, \nu_{\text{bal}}, \nu_{\text{flow}}, \eta_{\text{ref}}]\]

with total dimension $5n + 6m + 3k + n_{\text{ref}}$.

KKT Conditions

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 \operatorname{diag}(\mathrm{sw}) (\gamma_{\text{ub}} - \gamma_{\text{lb}}) = 0$, where $E_{\text{ref}} \in \mathbb{R}^{n \times n_{\text{ref}}}$ selects the per-island reference buses
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 \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$

1. Primal feasibility: $G_{\text{inc}} g + \text{psh} - d - B\theta = 0$
2. Flow definition: $f - WA\theta = 0$
3. Reference buses: $\theta_{\text{refs}} = 0$

Analytical Sparse KKT Jacobian

The KKT Jacobian $\partial K / \partial z$ is computed analytically as a sparse matrix, which is more efficient than automatic differentiation for large problems.

Parameter Derivatives

For each parameter type, we compute $\partial K / \partial p$ and then the full derivative via the IFT formula.

Demand ($d$)

Demand enters the power balance and the load shedding upper bound through $d_+ = \max(d, 0)$:

\[\frac{\partial K_{\nu_{\text{bal}}}}{\partial d} = -I, \qquad \frac{\partial K_{\mu_{\text{ub}}}}{\partial d} = \operatorname{diag}\left(\mu_{\text{ub}} \circ \frac{\partial d_+}{\partial d}\right)\]

For strictly positive demand, $\partial d_+ / \partial d = 1$. For negative demand, it is $0$. At zero demand, the positive part function is nonsmooth; the implementation uses the fixed zero shedding convention already required by the collapsed bound $0 \leq \text{psh} \leq 0$.

Switching ($\mathrm{sw}$)

Switching enters the Laplacian $B$, the weight matrix $W$, and the angle difference limits gated by $\mathrm{sw}$, with elementary perturbations:

\[\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\]

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$:

\[\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$).

Sensitivities are defined while the energized island partition is fixed. Opening or closing a bridge splits or merges islands and changes the reference set, so the derivative is nonsmooth at that topology boundary.

Cost Coefficients ($c_q$, $c_l$)

Quadratic cost $c_q$ enters stationarity w.r.t. $g$:

\[\frac{\partial K_g}{\partial c_{q,i}} = 2 g_i e_i\]

Linear cost $c_l$ enters stationarity w.r.t. $g$:

\[\frac{\partial K_g}{\partial c_l} = I_k\]

Flow Limits ($f_{\max}$)

Flow limits enter the complementary slackness conditions:

\[\frac{\partial K_{\lambda_{\text{lb}}}}{\partial f_{\max}} = -\operatorname{diag}(\lambda_{\text{lb}}), \qquad \frac{\partial K_{\lambda_{\text{ub}}}}{\partial f_{\max}} = \operatorname{diag}(\lambda_{\text{ub}})\]

Susceptances ($b$)

Susceptances affect the same blocks as switching (through $B$ and $W$), but with different partial derivatives since $B = A^\top \operatorname{diag}(-b \circ \mathrm{sw}) A$.

LMP Decomposition

Locational marginal prices are the power balance duals $\nu_{\text{bal}}$, decomposed as:

\[\text{LMP} = \underbrace{\text{energy}}_{\text{per island uniform component}} + \underbrace{\text{congestion}}_{\text{flow limit component}}\]

The congestion component is extracted by solving:

\[\text{congestion}[\text{non-ref}] = B_r^{-1} \left(A_r^\top W (\lambda_{\text{ub}}^{\text{std}} - \lambda_{\text{lb}}^{\text{std}}) + A_r^\top \operatorname{diag}(\mathrm{sw}) (\gamma_{\text{ub}}^{\text{std}} - \gamma_{\text{lb}}^{\text{std}})\right)\]

where $\lambda^{\text{std}}$, $\gamma^{\text{std}}$ use the standard sign convention (non-negative for binding constraints). The $\gamma$ terms capture congestion from binding phase angle difference limits. The energy component is uniform within each energized island and reflects that island's marginal cost of generation.