A Gentle Introduction to Optimal Power Flow

In an earlier blog post, we discussed the power flow problem, which serves as the key component of a much more challenging task: the optimal power flow (OPF). OPF is an umbrella term that covers a wide range of constrained optimization problems, the most important ingredients of which are: variables that optimize an objective function, some equality constraints, including the power balance and power flow equations, and inequality constraints, including bounds on the variables. Depending on the actual type of OPF, the form of the objective, as well as the sets of variables and constraints, will vary.

In the following, first, we derive the most widely used OPF variant, the economic dispatch problem, which aims to find a minimum cost solution of power demand and supply in the electricity grid, and show how to solve it. Then, we will briefly review other basic types of OPF. We will see that—unlike the power flow problem—OPF is in general rather expensive to solve, so we will finally discuss approximations that simplify the original model (known as relaxations) as well as computational intelligence and machine learning techniques that try to reduce real-time computational costs.

The power flow problem

In a previous blog post, we introduced the concept of the power flow problem and derived the equations for a simple and for an improved model. Since OPF relies on the power flow problem, we present here the final equations of the improved power flow model.

Let $\mathcal{N}$ denote the set of buses (represented by nodes in the corresponding graph of the power grid), and $\mathcal{E}, \mathcal{E}^{R} \subseteq \mathcal{N} \times \mathcal{N}$ the sets of forward and reverse orientations of transmission lines or branches (represented by directed edges in the corresponding graph of the power grid). Further, let $\mathcal{G}_{i}$, $\mathcal{L}_{i}$ and $\mathcal{S}_{i}$ denote the sets of generators, loads and shunt elements, respectively, belonging to bus $i$. Finally, the sets of all generators, loads and shunt elements are designated by $\mathcal{G} = \bigcup \limits_{i \in \mathcal{N}} \mathcal{G}_{i}$, $\mathcal{L} = \bigcup \limits_{i \in \mathcal{N}} \mathcal{L}_{i}$ and $\mathcal{S} = \bigcup \limits_{i \in \mathcal{N}} \mathcal{S}_{i}$, respectively.

Power flow models explicitly describe the power balance at each bus, based on Tellegen's theorem: \begin{equation} S_{i}^{\mathrm{gen}} - S_{i}^{\mathrm{load}} - S_{i}^{\mathrm{shunt}} = S_{i}^{\mathrm{trans}} \qquad \forall i \in \mathcal{N}. \label{balance} \end{equation} Eq. \ref{balance} states that the injected power by the connected generators ($S_{i}^{\mathrm{gen}}$) into the bus, reduced by the connected loads ($S_{i}^{\mathrm{load}}$) and shunt elements ($S_{i}^{\mathrm{shunt}}$), is equal to the transmitted power from and to the adjacent buses ($S_{i}^{\mathrm{trans}}$). There are two basic power flow models: the bus injection model (BIM) and the branch flow model (BFM) that are based on the undirected and directed graph representation of the power grid, respectively. For convenience, in this discussion we consider the BFM but note that the BIM can be also used to derive an equivalent OPF. For an improved power flow model, the BFM has four sets of equations: the first set expresses the power balance for each bus, the second set is the series current of the $\pi$-section model for each transmission line and the last two sets are the (asymmetric) power flows for the two directions of each transmission line: \begin{equation} \begin{aligned} & \sum \limits_{g \in \mathcal{G}_{i}} S_{g}^{G} - \sum \limits_{l \in \mathcal{L}_{i}} S_{l}^{L} - \sum \limits_{s \in \mathcal{S}_{i}} \left( Y_{s}^{S} \right)^{*} \lvert V_{i} \rvert^{2} \\ & \qquad \quad \ = \sum \limits_{(i, j) \in \mathcal{E}} S_{ij} + \sum \limits_{(k, i) \in \mathcal{E}^{R}} S_{ki} \hspace{3.0em} \forall i \in \mathcal{N}, \\ & I_{ij}^{s} = Y_{ij} \left( \frac{V_{i}}{T_{ij}} - V_{j} \right) \hspace{6.0em} \forall (i, j) \in \mathcal{E}, \\ & S_{ij} = V_{i} \left( \frac{I_{ij}^{s}}{T_{ij}^{*}} + Y_{ij}^{C} \frac{V_{i}}{\lvert T_{ij} \rvert^{2}}\right)^{*} \hspace{2.5em} \forall (i, j) \in \mathcal{E}, \\ & S_{ji} = V_{j} \left( -I_{ij}^{s} + Y_{ji}^{C}V_{j} \right)^{*} \hspace{4.5em} \forall (j, i) \in \mathcal{E}^{R}, \\ \end{aligned} \label{bfm_concise} \end{equation} where $S_{g}^{G}$ and $S_{l}^{L}$ represent the complex powers of generator $g$ and load $l$, respectively, $Y_{s}^{S}$ is the admittance of shunt element $s$, and $V_{i}$ is the complex voltage of bus $i$. Also, $I_{ij}^{s}$ and $S_{ij}$ denote the series current, power flow, respectively; $Y_{ij}$ and $Y_{ij}^{C}$ represent the series admittance, shunt admittance, respectively and $T_{ij}$ designates the complex tap ratio of the $\pi$-section model of branch $i \to j$.

Table 1. Comparison of improved BIM and BFM formulations of power flow: complex variables, number of complex variables and equations. Corresponding real variables, their numbers, as well as the number of real equations are also shown in parentheses. $N = \lvert \mathcal{N} \rvert$, $E = \lvert \mathcal{E} \rvert$, $G = \lvert \mathcal{G} \rvert$ and $L = \lvert \mathcal{L} \rvert$ denote the total number of buses, transmission lines, generators and loads, respectively. $v_{i}$, $\delta_{i}$, $i_{ij}^{s}$ and $\gamma_{ij}^{s}$ designate the voltage magnitude and angle of bus $i$ and the magnitude and angle of series current flowing from bus $i$ to bus $j$, respectively. $p$ and $q$ denote the active and reactive components of the corresponding complex powers, respectively.
Formulation Variables Number of variables Number of equations
BIM $V_{i} \ (v_{i}, \delta_{i})$
$S_{g}^{G} \ (p_{g}^{G}, q_{g}^{G})$
$S_{l}^{L} \ (p_{l}^{L}, q_{l}^{L})$
$N + G + L$
$(2N + 2G + 2L)$
$N$
$(2N)$
BFM $V_{i} \ (v_{i}, \delta_{i})$
$S_{g}^{G} \ (p_{g}^{G}, q_{g}^{G})$
$S_{l}^{L} \ (p_{l}^{L}, q_{l}^{L})$
$I_{ij}^{s} \ (i_{ij}^{s}, \gamma_{ij}^{s})$
$S_{ij} \ (p_{ij}, q_{ij})$
$S_{ji} \ (p_{ji}, q_{ji})$
$N + G + L + 3E$
$(2N + 2G + 2L + 6E)$
$N + 3E$
$(2N + 6E)$

The economic dispatch problem

Our starting OPF model is a specific variant of the economic dispatch (ED) model, which tries to find the economically optimal (i.e. lowest possible cost) output of the generators that meets the load and physical constraints of the power system. First, we discuss the elements of an ED formulation using the improved BMF model. Then we discuss the interior point method, which is one of the most widely used techniques to solve ED problems.

Variables

In power flow problems, the number of equations determines the maximum number of unknown variables for which we can solve the non-linear system. For instance, based on Table 1 it means that $G+L$ number of complex-valued (or $2G + 2L$ real-valued) variables need to be specified for both BIM and BFM formulations.

In ED, the variables are basically the same as in the power flow problems, but, as an optimization problem, ED has a higher number of unknown variables. For generality, in this discussion we treat all variables in ED formally as unknown. For variables whose value is actually specified, we can simply define an equality constraint or, more typically, the corresponding lower and upper bounds with the same specific value. This latter formalism makes the model construction more straightforward, and it also follows general practice [Coffrin18]. For the following sections, let $X$ denote all variables in ED.

Objective function

The objective (or cost) function of the most widely used economic dispatch OPF is the cost of the total real power generation. Let $C_{g} = C_{g}(p_{g}^{G})$ denote the individual cost curve of generator, $g \in \mathcal{G}$, that is a function solely depending on the (active) power generated in generator $g$. $C_{g}$ is a monotonically increasing function of $p_{g}^{G}$, usually modeled by a quadratic or piecewise linear function (Figure 1). The objective function, $C(X)$, then can be written as: \begin{equation} C(X) = \sum \limits_{g \in \mathcal{G}} C_{g}(p_{g}^{G}). \end{equation} We note that there can be alternative objective functions in OPF that can be used to investigate different features of the power grid [AlRashidi09]. For instance, one can minimize the line loss (i.e. power loss on transmission lines) over the network. Also, the objective function can express additional environmental and economical costs other than pure generation cost. As an example, with the increasing concerns over climate change, multiple studies have suggested that reducing overall emissions from generation should be also taken into account in the objective function [Gholami14].

Figure 1. Nonlinear generation cost curve (solid line) and its piecewise linear approximation (dashed lines) by three straight-line segments. $p^{\mathrm{min}}$ and $p^{\mathrm{max}}$ specify the operation limits.

Equality constraints

Equality constraints are expressions of the optimization variables characterized by equations. Equivalently, one can think of an equality constraint as a function of the optimization variables whose function value is fixed.

Since power flow equations express fundamental physical laws, they must naturally be contained in any OPF. The power balance equations together with the power flow equations (eq. $\ref{bfm_concise}$) are treated as equality constraints in OPF.

In the equations of power flow, the voltage angle differences play a key role, rather than the individual angles. Therefore, without loss of generality we can select a specific value of the voltage angle of a bus as reference and apply an additional equality constraint to this reference bus by setting its voltage angle to zero. We call this the slack or reference bus. In a more general setting, there can be multiple slack buses (i.e. buses used to make up any generation and demand mismatch caused by line losses) and corresponding constraints on their voltage angle. Denoting the set of all reference/slack buses by $\mathcal{R}$ we can write this set of constraints for all $r \in \mathcal{R}$ as: $\delta_{r} = 0$.

Inequality constraints

Finally, OPF also includes a set of inequality constraints. Similarly to equality constraints, an inequality constraint can be described by a function of the optimization variables, whose function value is allowed to change within a specified interval. Inequality constraints usually express some physical limitations of the elements of the system. For example, the most typical constraints for the ED problem are:

  • lower and upper bounds on voltage magnitude of the buses: $\underline{v}_{i} \le v_{i} \le \overline{v}_{i} \ \ \ \ \forall i \in \mathcal{N}$,
  • lower and upper bounds on voltage angle difference between adjacent buses: $\underline{\delta}_{ij} \le \delta_{ij} \le \overline{\delta}_{ij} \ \ \ \ \forall (i, j) \in \mathcal{E}$,
  • lower and upper bounds on output active power of generators: $\underline{p}_{i}^{G} \le p_{i}^{G} \le \overline{p}_{i}^{G} \ \ \ \ \forall i \in \mathcal{N}$,
  • lower and upper bounds on output reactive power of generators: $\underline{q}_{i}^{G} \le q_{i}^{G} \le \overline{q}_{i}^{G} \ \ \ \ \forall i \in \mathcal{N}$,
  • upper bounds on the absolute value of active power flow of transmission lines: $\lvert p_{ij}\rvert \le \overline{p}_{ij} \ \ \ \ \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}$,
  • upper bounds on the absolute value of reactive power flow of transmission lines: $\lvert q_{ij}\rvert \le \overline{q}_{ij} \ \ \ \ \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}$,
  • upper bound on the absolute (square) value of complex power flow of transmission lines: $\lvert S_{ij}\rvert^{2} = p_{ij}^{2} + q_{ij}^{2} \le \overline{S}_{ij}^{2} \ \ \ \ \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}$.
We note that, in actual implementation, some inequality constraints can appear in alternative forms. For instance, lower and upper bounds on power flows can be replaced by a corresponding set of constraints on the currents.

The economic dispatch problem using BFM

Putting everything together, the form of the economic dispatch OPF using BFM is the following: \begin{equation} \begin{aligned} & \underline{\mathrm{variables}}{:} \\ & \quad X^{\mathrm{BFM}} \begin{cases} V_{i} = (v_{i}, \delta_{i}) \hspace{11.0em} \forall i \in \mathcal{N}, \\ S_g^{G} = (p_{g}^{G}, q_{g}^{G}) \hspace{9.75em} \forall g \in \mathcal{G}, \\ S_{ij} = (p_{ij}, q_{ij}) \hspace{8.5em} \forall (i, j) \in \mathcal{E}, \\ S_{ji} = (p_{ji}, q_{ji}) \hspace{8.5em} \forall (j, i) \in \mathcal{E}^{R}, \\ I_{ij}^{s} = (i_{ij}^{s}, \gamma_{ij}^{s}) \hspace{8.75em} \forall (i, j) \in \mathcal{E}; \\ \end{cases} \\ & \underline{\mathrm{objective\ funcion}}{:} \\ & \quad \min \limits_{X^{\mathrm{BFM}}} \sum \limits_{g \in \mathcal{G}} C_{g}(p_{g}^{G}); \\ & \underline{\mathrm{subject\ to}}{:} \\ & \quad \mathcal{C}^{\mathrm{PB}} \hspace{1.0em} \begin{cases} \sum \limits_{g \in \mathcal{G}_{i}} S_{g}^{G} - \sum \limits_{l \in \mathcal{L}_{i}} S_{l}^{L} - \sum \limits_{s \in \mathcal{S}_{i}} \left( Y_{s}^{S} \right)^{*} \lvert V_{i} \rvert^{2} \\ \hspace{2.6em} \ = \sum \limits_{(i, j) \in \mathcal{E}} S_{ij} + \sum \limits_{(k, i) \in \mathcal{E}^{R}} S_{ki} \hspace{3.25em} \forall i \in \mathcal{N}; \\ \end{cases} \\ & \quad \mathcal{C}^{\mathrm{REF}} \hspace{0.75em} \begin{cases} \hspace{0.25em} \delta_{r} = 0 \hspace{13.4em} \forall r \in \mathcal{R}; \\ \end{cases} \\ & \quad \mathcal{C}^{\mathrm{BFM}} \hspace{0.25em} \begin{cases} I_{ij}^{s} = Y_{ij} \left( \frac{V_{i}}{T_{ij}} - V_{j} \right) \hspace{6.5em} \forall (i, j) \in \mathcal{E}, \\ S_{ij} = V_{i} \left( \frac{I_{ij}^{s}}{T_{ij}^{*}} + Y_{ij}^{C} \frac{V_{i}}{\lvert T_{ij} \rvert^{2}}\right)^{*} \hspace{3.75em} \forall (i, j) \in \mathcal{E}, \\ S_{ji} = V_{j} \left( -I_{ij}^{s} + Y_{ji}^{C}V_{j} \right)^{*} \hspace{4.25em} \forall (j, i) \in \mathcal{E}^{R}; \\ \end{cases} \\ & \quad \mathcal{C}^{\mathrm{INEQ}} \begin{cases} \underline{v}_{i} \le v_{i} \le \overline{v}_{i} \hspace{11.25em} \forall i \in \mathcal{N}, \\ \underline{p}_{i}^{G} \le p_{i}^{G} \le \overline{p}_{i}^{G} \hspace{10.25em} \forall i \in \mathcal{N}, \\ \underline{p}_{i}^{G} \le p_{i}^{G} \le \overline{p}_{i}^{G} \hspace{10.25em} \forall i \in \mathcal{N}, \\ \underline{\delta}_{ij} \le \delta_{ij} \le \overline{\delta}_{ij} \hspace{8.75em} \forall (i, j) \in \mathcal{E}, \\ \lvert p_{ij}\rvert \le \overline{p}_{ij} \hspace{10.5em} \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}, \\ \lvert q_{ij}\rvert \le \overline{q}_{ij} \hspace{10.6em} \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}, \\ \lvert S_{ij}\rvert^{2} = p_{ij}^{2} + q_{ij}^{2} \le \overline{S}_{ij}^{2} \hspace{5.2em} \forall (i, j) \in \mathcal{E} \cup \mathcal{E}^{R}; \\ \end{cases} \\ \end{aligned} \label{opf_bmf} \end{equation} where $\mathcal{C}^{\mathrm{PB}}$, $\mathcal{C}^{\mathrm{REF}}$, and $\mathcal{C}^{\mathrm{BFM}}$ denote the equality sets for the power balance, the voltage angle of reference buses, and the improved BFM, respectively, while $\mathcal{C}^{\mathrm{INEQ}}$ is the set of inequality constraints.

Solving the ED problem by the interior-point method

Economic dispatch is a non-linear and non-convex optimization problem. One of the most successful techniques to solve such large-scale programming is the interior-point method [Nocedal06][Wachter06]. In this section, we briefly discuss the key aspects of the method.

The ED problem can be expressed by the concise form of the general mathematical optimization, where the objective function $f(x)$ is minimized with respect to the optimization variables $x$ and subject to a set of equality and inequality constraints, represented by the $c^{\mathrm{E}}(x)$ and $c^{\mathrm{I}}(x)$ vectors, respectively: \begin{equation} \begin{aligned} & \min \limits_{x} \ f(x), \\ & \mathrm{s.t.} \; \; c^{\mathrm{E}}(x) = 0, \\ & \quad \; \; \; c^{\mathrm{I}}(x) \ge 0, \\ \end{aligned} \label{opt_problem} \end{equation} For example, using the BFM formalism, $x = X^{\mathrm{BFM}}$, $\ f(x) = C(x)$, with $\ c^{\mathrm{E}}_{i} \in \mathcal{C}^{E}$, where $\mathcal{C}^{\mathrm{E}} = \mathcal{C}^{\mathrm{PB}} \cup \mathcal{C}^{\mathrm{REF}} \cup \mathcal{C}^{\mathrm{BFM}}$, and $c^{\mathrm{I}}_{j} \in \mathcal{C}^{I}$, where $\mathcal{C}^{\mathrm{I}} = \mathcal{C}^{\mathrm{INEQ}}$.

Interior-point methods are also referred to as barrier methods as they introduce a surrogate model of the above optimization problem by converting the inequality constraints to equality ones and replacing the objective function with a logarithmic barrier function: \begin{equation} \begin{aligned} & \min \limits_{x, s} \ f(x) - \mu \sum \limits_{i=1}^{\lvert \mathcal{C}^{\mathrm{I}} \rvert} \log s_{i}, \\ & \mathrm{s.t.} \; \; c^{\mathrm{E}}(x) = 0, \\ & \quad \; \; \; c^{\mathrm{I}}(x) - s = 0, \\ \end{aligned} \label{barrier_problem} \end{equation} where $s$ is the vector of slack variables, and the second term in the objective function is called the barrier term, with $\mu > 0$ barrier parameter. We note that the barrier term implicitly imposes an inequality constraint on the slack variables, due to the logarithmic function: $s \ge 0$. It is easy to see that the barrier problem is not equivalent to the original non-linear program for $\mu > 0$. However, as $\mu \to 0$, the solution of the barrier problem approaches that of problem $\ref{opt_problem}$. Therefore, the idea of the barrier method is to find an approximate solution of the original problem by using an appropriate sequence of positive barrier parameters that converges to zero. The basic algorithm is based on the Karush–Kuhn–Tucker (KKT) conditions of the barrier problem. It consists of four sets of equations: the first two sets of equations are derived from the first-order condition of the system Lagrangian, i.e. its derivatives with respect to $x$ and $s$ are equal to zero, while the third and fourth sets of equations correspond to the equality and inequality constraints, respectively. \begin{equation} \begin{aligned} \nabla f(x) - J_{\mathrm{E}}^{T}(x)y - J_{\mathrm{I}}^{T}(x)z = 0, \\ Sz - \mu e = 0, \\ c^{\mathrm{E}}(x) = 0, \\ c^{\mathrm{I}}(x) - s = 0, \\ \end{aligned} \label{kkt} \end{equation} where the $J_{\mathrm{E}}(x)$ and $J_{\mathrm{I}}(x)$ are the Jacobian matrices of the equality and of the inequality functions, respectively, and $y$ and $z$ are their corresponding Lagrange multipliers. Also, we use the following additional notations: $S = \mathrm{diag}(s)$ and $Z = \mathrm{diag}(z)$ are diagonal matrices, $I$ is the identity matrix, and $e = (1, \dots, 1)^{T}$.

The KKT conditions establish a relationship between the primal ($x$, $s$) and dual ($y$, $z$) variables, resulting in a non-linear system (eq. \ref{kkt}) that can be expressed by a vector-valued function: $F_{\mathrm{KKT}}(x, s, y, z) = 0$. In order to solve this non-linear system, we can apply the Newton-Raphson method and look for a search direction based on the following linear problem: \begin{equation} J_{\mathrm{KKT}}(x, s, y, z)p = -F_{\mathrm{KKT}}(x, s, y, z), \label{primal_dual} \end{equation} where $J_{\mathrm{KKT}}(x, s, y, z)$ and $p$ are the Jacobian matrix of the KKT system and search direction vector, respectively: \begin{equation} \begin{aligned} F_{\mathrm{KKT}}(x, s, y, z) & = \begin{pmatrix} \nabla f(x) - J_{\mathrm{E}}^{T}(x)y - J_{\mathrm{I}}^{T}(x)z \\ Sz - \mu e \\ c^{\mathrm{E}}(x) \\ c^{\mathrm{I}}(x) - s \\ \end{pmatrix}, \\ J_{\mathrm{KKT}}(x, s, y, z) & = \begin{pmatrix} \nabla_{xx}L & 0 & -J_{\mathrm{E}}^{T}(x) & -J_{\mathrm{I}}^{T}(x) \\ 0 & Z & 0 & S \\ J_{\mathrm{E}}(x) & 0 & 0 & 0 \\ J_{\mathrm{I}}(x) & -I & 0 & 0 \\ \end{pmatrix}, \\ p & = \begin{pmatrix} p_{x} \\ p_{s} \\ p_{y} \\ p_{z} \\ \end{pmatrix}, \end{aligned} \end{equation} where the Lagrangian of the system is defined as \begin{equation} L(x, s, y, z) = f(x) - y^{T} c^{\mathrm{E}}(x) - z^{T} \left( c^{\mathrm{I}}(x) - s \right), \label{lagrangian} \end{equation} and $\nabla_{xx}L$ is the Hessian of the Lagrangian with respect to the optimization variables. Eq. $\ref{primal_dual}$ is called the primal-dual system, as it includes the dual variables besides the primal ones. The primal-dual system is solved iteratively, i.e. in the $(n+1)$th iteration, after obtaining $p^{(n+1)}$, we can update the variables and the barrier parameter: \begin{equation} \begin{aligned} & x_{n+1} = x_{n} + \alpha_{x}^{(n+1)} p_{x}^{(n+1)}, \\ & s_{n+1} = s_{n} + \alpha_{s}^{(n+1)} p_{s}^{(n+1)}, \\ & y_{n+1} = y_{n} + \alpha_{y}^{(n+1)} p_{y}^{(n+1)}, \\ & z_{n+1} = z_{n} + \alpha_{z}^{(n+1)} p_{z}^{(n+1),} \\ & \mu_{n+1} = g(\mu_{n}, x_{n+1}, s_{n+1}, y_{n+1}, z_{n+1}), \\ \end{aligned} \end{equation} where $\alpha_i$ denotes the line search parameter for variable $i$, and for the barrier parameter we used a fairly general function $g$. There are several modern interior-point methods that propose different heuristic functions for computing line search and barrier parameters, in order to cope with non-convexities and non-linearities and improve convergence rate.

OPF variants

So far we have been focusing on one specific OPF problem, a variant of the economic dispatch that is one of the most basic types in the realm of OPF. In this section, we briefly discuss other models applied in the field.

In a more general setting, the task of economic dispatch is to optimize some economic utility function for participants, including both loads and generators. For instance, incorporating both the power supply and demand offers into the objective function leads to the basic concept of the electricity market. At the optimum, the amount of supply and demand power reflects simple market principles subject to the physical constraints of the system: power is primarily bought from generators offering the lowest prices and primarily sold to loads offering the highest prices.

In the so-called volt/var control we optimize the combination of power loss, peak demand, and power consumption, with respect to the voltage level and reactive powers.

In the above problems, the set of generators is specified and all of these generators are supposed to contribute to the total power generation. However, in real industrial problems, not all generators available in the grid need to be active. Generators differ in many ways besides their generation curve: they have different start up and shut down time scales (i.e. how quickly a generator can reach its operating level and how quickly it can stop running) and costs, ramp rates (i.e. the rate at which a generator can increase or decrease its output), minimum and maximum power generation, etc., that all need to be taken into account in the planning procedure. Therefore, it is very much desirable to allow specific generators to not run in the optimization model. In order to select the active generators, electricity grid operators solve the unit commitment (UC) problem. The unit commitment economic dispatch (UCED) problem can be considered as an extension of the ED problem, where the set of variables is extended with binary variables, each responsible for including or excluding the actual generator in the optimization. UC is, therefore, a mixed-integer non-linear programming, which is an even harder problem to solve computationally.

Another aspect that needs to be considered is the security of the grid. Power grids are complex systems that include thousands of buses and transmission lines, where failure of any element can happen occasionally and impact the operation of the power grid significantly. Thus, it is very important to account for such contingencies. The corresponding set of problems is called security constrained optimal power flow (SCOPF) and in combination with ED it is called security constrained economic dispatch (SCED). SCOPF provides a solution not only for the ideal unfailing case, but at the same time for some contingency scenarios as well. These problems are, therefore, much larger scale problems, including multiple times the number of variables and constraints, resulting again in a mixed-integer, non-linear programming. For example, the widely used $N-1$ contingency problem addresses all scenarios where a single component fails.

Finally, we note that a more general (and very large-scale) problem can be formulated by using both UC and SCOPF, leading to the security constrained unit commitment problem (SCUC).

Solving OPF problems

Solving AC-OPF problems by interior-point methods is widely used as a reliable baseline approach. However, solving such problems by interior-point methods is not cheap, as it requires the computation of the Hessian of the Lagrangian $\left( \nabla_{xx} L(x_{n}, s_{n}, y_{n}, z_{n}) \right)$ at each iteration step $n$ (eq. $\ref{primal_dual}$). Since the required computational time has a disadvantageous superlinear scaling with system size, this can be prohibitive for solving large scale problems. Moreover, the scaling is even worse for OPF problems that are mixed-integer programming.

There are several approaches that address this issue by either simplifying the original AC-OPF problem or by using other techniques to compute exact or approximate solutions.

Convex relaxations

The first set of approximations we discuss is called convex relaxations. The main idea of these approaches is to approximate the original non-convex optimization problem with a convex one. Unlike in non-convex problems, any local minimum is a global minimum as well in convex problems (since convex problems, by definition, can only have a single minimum). Because of this and other useful features, there are very straightforward algorithms with excellent convergence rate to solve convex optimization problems [Boyd04][Nocedal06]. To demonstrate convex relaxations, let us consider now the economic dispatch problem using the basic BIM formulation. As it can be shown, this problem—with some slight modification—can be reformulated as a quadratically constrained quadratic program (QCQP) using the complex voltages of the buses [Low14]. In the reformulated optimization problem, the objective function and all (inequality) constraints have a quadratic form in the complex voltages: \begin{equation} \begin{aligned} & \min \limits_{V \in \mathbb{C}^{N}} V^{H}CV, \\ & \mathrm{s.t.} \quad V^{H}M_{k}V \le m_{k} \quad k \in \mathcal{K}, \\ \end{aligned} \label{qcqp} \end{equation} where the superscript $H$ denotes conjugate transpose, $M_{k}$s are Hermitian matrices (i.e. $M_{k}^{H} = M_{k}$), and $m_{k}$s are real-valued vectors specifying the corresponding inequality constraints. QCQP is still a non-convex problem that is also an NP-hard problem. Let $W = VV^{H}$ denote a rank–1 (i.e. rank$(W) = 1$) positive semidefinite (i.e. $W \succeq 0$) matrix. Using the trace function and the identity $V^{H}MV = \mathrm{tr}(MVV^{H}) = \mathrm{tr}(MW)$ for any Hermitian matrix, we can transform the quadratic forms of eq. $\ref{qcqp}$ to expressions that include $W$: \begin{equation} \begin{aligned} & \min \limits_{W \in \mathbb{S}^{N}} \mathrm{tr}(CW), \\ & \mathrm{s.t.} \quad \mathrm{tr}(M_{k}W) \le m_{k} \quad k \in \mathcal{K}, \\ & \quad \quad \ W \succeq 0, \\ & \quad \quad \ \mathrm{rank}(W) = 1, \\ \end{aligned} \end{equation} where $\mathbb{S}^{N}$ denotes the space of $N \times N$ symmetric matrices. The above problem is convex in $W$, except for the rank–1 equality constraint. Removing this constraint leads to a semidefinite programming (SDP), and therefore it is called the SDP relaxation of the OPF. SDP introduces a convex superset of the optimization variables. Other choices of convex supersets yield different relaxations, like chordal or second-order cone programming (SOCP), that can be solved even more efficiently than SDP problems. We note that all of these relaxations provide a lower bound to OPF and for some sufficient conditions these relaxations can be even exact [Low14a]. Also, an advantage of convex relaxation approaches over other approximations is that if a relaxed problem is infeasible, then the original problem is also infeasible.

DC-OPF

One of the most widely used approximations is the DC-OPF approach. Its name is pretty misleading, as this method does not assume a direct current (instead of alternating current) setup, although it removes the reactive power flow equations (and corresponding variables and constraints). The DC approximation is, in fact, a linearization of the AC problem.

In order to demonstrate the simplifications introduced by the DC-OPF approach, we recall that for a simple power flow model the active and reactive powers flowing from bus $i$ to $j$ can be writtes as: \begin{equation} \left\{ \begin{aligned} p_{ij} & = \frac{1}{r_{ij}^{2} + x_{ij}^{2}} \left[ r_{ij} \left( v_{i}^{2} - v_{i} v_{j} \cos(\delta_{ij}) \right) + x_{ij} \left( v_{i} v_{j} \sin(\delta_{ij}) \right) \right], \\ q_{ij} & = \frac{1}{r_{ij}^{2} + x_{ij}^{2}} \left[ x_{ij} \left( v_{i}^{2} - v_{i} v_{j} \cos(\delta_{ij}) \right) + r_{ij} \left( v_{i} v_{j} \sin(\delta_{ij}) \right) \right]. \\ \end{aligned} \right. \label{power_flow_z} \end{equation} DC-OPF has three major assumptions [Liu09]:

  1. The resistance $r_{ij}$ of each branch is negligible, i.e. $r_{ij} \approx 0$.
  2. The bus voltage magnitudes are everywhere approximated by one per unit, i.e. $v_{i} \approx 1$.
  3. The voltage angle difference of each branch is very small, i.e. $\cos(\delta_{ij}) \approx 1$ and $\sin(\delta_{ij}) \approx \delta_{ij}$.
Starting from eq. $\ref{power_flow_z}$ and applying the above approximations, we get the following simplified expressions: \begin{equation} \left\{ \begin{aligned} p_{ij} & = \frac{1}{r_{ij}^{2} + x_{ij}^{2}} \left[ r_{ij} \left( v_{i}^{2} - v_{i} v_{j} \cos(\delta_{ij}) \right) + x_{ij} \left( v_{i} v_{j} \sin(\delta_{ij}) \right) \right] \\ & \overset{\mathrm{1.}}{\approx} \frac{1}{x_{ij}} \left( v_{i}v_{j}\sin(\delta_{ij}) \right) \overset{\mathrm{2.}}{\approx} \frac{\sin(\delta_{ij})}{x_{ij}} \overset{\mathrm{3.}}{\approx} \frac{\delta_{ij}}{x_{ij}}, \\ q_{ij} & = \frac{1}{r_{ij}^{2} + x_{ij}^{2}} \left[ x_{ij} \left( v_{i}^{2} - v_{i} v_{j} \cos(\delta_{ij}) \right) + r_{ij} \left( v_{i} v_{j} \sin(\delta_{ij}) \right) \right] \\ & \overset{\mathrm{1.}}{\approx} \frac{1}{x_{ij}} \left( v_{i}^{2} - v_{i} v_{j} \cos(\delta_{ij}) \right) \overset{\mathrm{2.}}{\approx} \frac{1 - \cos(\delta_{ij})}{x_{ij}} \overset{\mathrm{3.}}{\approx} 0. \end{aligned} \right. \end{equation} As we can see, DC-OPF reduces the original power flow problem significantly by removing several variables and constraints. Also, the remaining set of constraints with a linear objective function results in a linear programming that can be solved very efficiently by interior-point or simplex methods. DC-OPF has been demonstrated pretty useful for a wide variety of applications. However, it also has several limitations among which the most severe one is that the solution of DC-OPF may not even be a feasible solution of AC-OPF. In such situations, some constraints need to be tightened, and the DC-OPF rerun.

Computational intelligence and machine learning

The increasing uncertainty in grid conditions due to the integration of renewable resources (such as wind and solar) requires OPF problems to be solved near real-time, in order to have the most accurate inputs reflecting the latest state of the system. This in turn requires the grid operators to have the computational capacity of running consecutive instances of OPF problems with fast convergence time. However, this is a rather challenging task for large-scale OPF problems, even using the DC-OPF approximation. Therefore, as an alternative to interior-point methods, there has been intense research in the computational intelligence [AlRashidi09] and machine learning [Hasan20] fields. Below, we give a brief list of the most important approaches.

Computational intelligence [AlRashidi09] techniques for solving OPF include:

Among machine learning algorithms [Hasan20], reinforcement, unsupervised, and supervised (both regression and classification) learning approaches are used to solve OPF problems. In general, these techniques shift the computational effort away from real-time to offline training. The advantage of taking this approach is that, while training the machine learning algorithms may be a time-consuming process, once trained the evaluation is very fast. In this way, once the algorithm is trained once it can be used multiple times for extremely fast evaluation. Their most widely used applications include:
  • direct prediction of OPF variables;
  • prediction of set-points of OPF variables for subsequent warm-start of interior-point methods;
  • prediction of the active constraints of OPF in order to build and solve significantly smaller reduced OPF problems.
In a subsequent post we present a detailed description of the above techniques using neural networks.

References

[Coffrin18]: C. Coffrin, R. Bent, K. Sundar, Y. Ng and M. Lubin, "PowerModels. JL: An Open-Source Framework for Exploring Power Flow Formulations," Power Systems Computation Conference (PSCC), Dublin, pp. 1, (2018).

[AlRashidi09]: M. R. AlRashidi and M. E. El-Hawary, "Applications of computational intelligence techniques for solving the revived optimal power flow problem," in Electric Power Systems Research, 79, (2009).

[Gholami14]: A. Gholami, J. Ansari, M. Jamei and A. Kazemi, "Environmental/economic dispatch incorporating renewable energy sources and plug-in vehicles," in IET Generation, Transmission & Distribution, 8, pp. 2183, (2014).

[Nocedal06]: J. Nocedal and S. J. Wrigh, "Numerical Optimization," New York: Springer, (2006).

[Wachter06]: A. Wächter and L. Biegler, "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming," Math. Program. 106, pp. 25, (2006).

[Boyd04]: S. Boyd and L. Vandenberghe, "Convex Optimization, " New York: Cambridge University Press, (2004).

[Low14]: S. H. Low, "Convex Relaxation of Optimal Power Flow—Part I: Formulations and Equivalence," in IEEE Transactions on Control of Network Systems, 1, pp. 15, (2014).

[Low14a]: S. H. Low, "Convex Relaxation of Optimal Power Flow—Part II: Exactness," in IEEE Transactions on Control of Network Systems, 1, pp. 177, (2014).

[Liu09]: H. Liu, L. Tesfatsion and A. A. Chowdhury, "Locational marginal pricing basics for restructured wholesale power markets," IEEE Power & Energy Society General Meeting, Calgary, AB, pp. 1, (2009).

[Hasan20]: F. Hasan, A. Kargarian and A. Mohammadi, "A Survey on Applications of Machine Learning for Optimal Power Flow," IEEE Texas Power and Energy Conference (TPEC), College Station, TX, USA, pp. 1, (2020).