(home)

Writing DC Optimal Power Flow in Matrix Form

The DC optimal power flow (DC OPF) approximation is a widely used model for power system optimization. In this post, we will write the DC OPF problem in matrix form and differentiate its solution with respect to the network parameters. This post follows the notation from our recent e-Energy paper.

Parameterized DC OPF program formulation

The program parameters are fixed problem data collected in the vector $\boldsymbol{\theta} := (\boldsymbol{d}^\top,\boldsymbol{\alpha}^\top,\boldsymbol{\beta}^\top)^\top$. The parameters $\boldsymbol{d} \in \mathbb{R}^{N}$ are the transmission-level active power demands $d_k = \sum_{i \in \mathcal{ D}_k} d_{k,i}$; $\boldsymbol{\alpha}\in \mathbb{R}^K$ and $\boldsymbol{\beta} \in \mathbb{R}^K$ are the quadratic and linear generator cost coefficients, respectively. The parameterized DC OPF program is then

$$\min_{\boldsymbol{g},\boldsymbol{p}} \ \sum_{i=1}^k \alpha_i g_i^2 + \beta_i g_i,\\ \operatorname{\sf s.t.} \quad \boldsymbol{F}(\boldsymbol{B} \boldsymbol{g} - \boldsymbol{d}) = \boldsymbol{p} \\ \quad \quad \mathbf{ 1}^\top\boldsymbol{B} \boldsymbol{g} = \mathbf{ 1}^\top\boldsymbol{d}, \\ \quad\quad -\boldsymbol{p}^{\sf max} \leq \boldsymbol{p} \leq \boldsymbol{p}^{\sf max} \\ \quad\quad \mathbf{ 0} \leq \boldsymbol{g} \leq \boldsymbol{g}^{\sf max},$$

where $\boldsymbol{g} \in \mathbb{R}^K$ are the power generation set points and $\boldsymbol{p} \in \mathbb{R}^M$ are the line power flows. The matrix $\boldsymbol{B} \in \{0,1\}^{N \times K}$ is a generator-to-node incidence matrix with entries of the form

$$B_{i,j} = \begin{cases} 1 & \mathsf{generator} \ j \ \mathsf{at \ node} \ i\\ 0 & \mathsf{otherwise,} \end{cases}$$

and the matrix $\boldsymbol{F} \in \mathbb{R}^{M \times N}$ is the power transfer distribution factor (PTDF) matrix.

Quadratic program formulation of DC OPF

With some algebra, the parameterized DC OPF program (1) can be written as a compact, general form quadratic program (QP). To achieve this, define the constraint matrices

$$\boldsymbol{A} \triangleq\begin{pmatrix} \boldsymbol{F} \boldsymbol{B} & - \mathbf{ I}_M\\ \mathbf{ 1}^\top \boldsymbol{B} & \mathbf{ 0}_{M}^\top \end{pmatrix} \quad \boldsymbol{G} \triangleq\begin{pmatrix} -\mathbf{ I}_K & \mathbf{ 0}_{K\times M}\\ \mathbf{ I}_K & \mathbf{ 0}_{K\times M}\\ \mathbf{ 0}_{M\times K} & -\mathbf{ I}_M\\ \mathbf{ 0}_{M\times K} & \mathbf{ I}_M \end{pmatrix},$$

and the constraint vectors

$$\boldsymbol{x} \triangleq \begin{pmatrix} \boldsymbol{g}\\ \boldsymbol{p} \end{pmatrix}, \ \ \boldsymbol{y} \triangleq \begin{pmatrix} \boldsymbol{F} \boldsymbol{d}\\ \mathbf{ 1}^\top \boldsymbol{d} \end{pmatrix}, \ \ \boldsymbol{h} \triangleq \begin{pmatrix} \mathbf{ 0}\\ \boldsymbol{g}^{\sf max}\\ \boldsymbol{p}^{\sf max}\\ \boldsymbol{p}^{\sf max} \end{pmatrix}.$$

Furthermore, define the quadratic cost coefficient matrix and linear cost vector as

$$\boldsymbol{Q} \triangleq \begin{pmatrix} \operatorname{\sf diag}(\boldsymbol{\alpha}) & \mathbf{ 0}_{K \times M}\\ \mathbf{ 0}_{M \times K} & \tau \mathbf{ I}_{M \times M} \end{pmatrix}, \ \ \boldsymbol{w} \triangleq \begin{pmatrix} \boldsymbol{\beta}\\ \mathbf{ 0}_M \end{pmatrix},$$

where $\tau>0$ is a small regularization parameter to yield a strongly convex objective. The program (1) can now be written as

$$\min_{\boldsymbol{x}}\quad \frac{1}{2}\boldsymbol{x}^\top \boldsymbol{Q} \boldsymbol{x} + \boldsymbol{w}^\top \boldsymbol{x},\\%\boldsymbol{w}^\top \boldsymbol{x},\\ \mathsf{s.t.}\quad \boldsymbol{G} \boldsymbol{x} \preceq \boldsymbol{h}\\ \quad \boldsymbol{A} \boldsymbol{x} = \boldsymbol{y}.$$

Let $\boldsymbol{\mu} \in \mathbb{R}^{2K + 2M}$ and $\boldsymbol{\nu} \in \mathbb{R}^{M+1}$ be the dual variables associated with the inequality constraints $\boldsymbol{G} \boldsymbol{x} \preceq \boldsymbol{h}$ and equality constraints $\boldsymbol{A} \boldsymbol{x} = \boldsymbol{y}$, respectively. Let $\boldsymbol{z} := (\boldsymbol{x}^\top,\boldsymbol{\mu}^\top,\boldsymbol{\nu}^\top)^\top$ be the concatenation of the primal and dual variables of the program (6). The Lagrangian of the program is

$$l(\boldsymbol{z};\boldsymbol{\theta}) \triangleq \frac{1}{2}\boldsymbol{x}^\top \boldsymbol{Q} \boldsymbol{x} + \boldsymbol{w}^\top \boldsymbol{x} + \boldsymbol{\mu}^\top(\boldsymbol{G}\boldsymbol{x} - \boldsymbol{h}) + \boldsymbol{\nu}^\top(\boldsymbol{A}\boldsymbol{x} - \boldsymbol{y}).$$

Next, we construct a linear system of equations that describes the Karush-Kuhn-Tucker conditions for optimality of the program (6). Let $\boldsymbol{z}^*$ be a solution to the program (6). The stationarity, complimentary slackness, and equality constraint feasibility conditions form a system of equations $\boldsymbol{k}(\cdot)$ whose root is $\boldsymbol{z}^*$. The KKT operator $\boldsymbol{k}(\cdot)$ is then

$$\boldsymbol{k}(\boldsymbol{z};\boldsymbol{\theta}) = \begin{pmatrix} \nabla_{\boldsymbol{x}}\mathcal{L}(\boldsymbol{z};\boldsymbol{\theta})\\ \operatorname{\sf diag}(\boldsymbol{\mu})(\boldsymbol{G} \boldsymbol{x} - \boldsymbol{h}) \\ \boldsymbol{A} \boldsymbol{x} - \boldsymbol{y} \end{pmatrix} = \begin{pmatrix} \boldsymbol{Q} \boldsymbol{x} + \boldsymbol{w} + \boldsymbol{G}^\top \boldsymbol{\mu} + \boldsymbol{A}^\top \boldsymbol{\nu}\\ \operatorname{\sf diag}(\boldsymbol{\mu})(\boldsymbol{G} \boldsymbol{x} - \boldsymbol{h})\\ \boldsymbol{A} \boldsymbol{x} - \boldsymbol{y} \end{pmatrix},$$

which satisfies $\boldsymbol{k}(\boldsymbol{z}^*;\boldsymbol{\theta}) = \mathbf{ 0}$.

Linear Programming Representation of DC OPF

An even simpler model can be obtained by eliminating the quadratic terms in the cost function. This is achieved by introducing a new variable $\boldsymbol{z} \in \mathbb{R}^{K+M}$ and rewriting the program (6) as a linear program (LP). The LP is then

$$\min_{\boldsymbol{x},\boldsymbol{z}} \quad \boldsymbol{w}^\top \boldsymbol{x},\\ \mathsf{s.t.} \quad \boldsymbol{Q} \boldsymbol{x} + \boldsymbol{G}^\top \boldsymbol{z} = -\boldsymbol{w},\\ \boldsymbol{A} \boldsymbol{x} = \boldsymbol{y},\\ \boldsymbol{G} \boldsymbol{x} \preceq \boldsymbol{h}.$$