An Iterative Approach to Finding Global Solutions of AC Optimal Power Flow Problems

Optimal power flow (OPF) is a fundamental resource allocation problem in power system operations that minimizes the cost of power generation while satisfying demand. The ACOPF formulation of the problem uses nonlinear power flow equations, resulting in nonlinear and nonconvex optimization problems [1, 2, 3]. The consequence of the nonconvexity of ACOPF we study in this paper is the presence of multiple solutions. The existence of multiple local solutions of the OPF problem has been well-known for several decades [4]. A common assumption is that OPF problems tend to have a single “practical” solution, and therefore the fact that multiple solutions can exist do not impact day-to-day operations [5, 6]. However, an increasingly large body of work have pointed to that multiple solutions do occur under reasonable conditions and cannot easily ruled out [7, 8, 2]. For example, [7] shows how modifications of the standard IEEE benchmarks can lead to each having more than one local solutions.

Most ACOPF problems are solved via variations of nonlinear optimization algorithms, including Newton-Raphson, sequential programming, interior points and others (see [2, 9, 10] and the references within). These algorithms are in general only able to certify a solution is locally optimal. These solutions satisfy some necessary optimality conditions (e.g., the KKT conditions) but may not globally optimal, that is, the least cost solution. An open question in the field is to develop algorithms that can find global optimal solutions. In addition, understanding and distinguishing the local and the global optimal solutions can lead to important theoretical discoveries about the ACOPF problem.

In this paper, we propose a simple algorithm that can effectively escape from local solutions to seek global optimal solutions. The process is outlined in Fig. 1. First, we solve the ACOPF problem using an existing solver (e.g., IPOPT [11]). From the solution and its associated dual variables, we form a partial Lagrangian. This partial Lagrangian serves to “flatten” the geometry of the optimization problem. We then optimize this partial Lagrangian, which can lead to a different solution. Using this second solution as a warm start, we again call the solver for the ACOPF problem. Repeating this iterative process, we can successively improve the solution quality, moving from local solutions to globally optimal ones.

We show the performance of our algorithm on standard 9-bus, 22-bus and 39-bus networks. We show that our algorithm can escape from local solutions [7] and find the global optimal solution a single iteration for most instances. We also sketch the intuition behind the algorithm using a two bus system.

Our approach can be seen as a way to provide good warm starts to nonlinear optimization solvers. Existing approaches along this line either randomizes [12] or uses a previous solution as the starting point [13]. The former tends to be time consuming, while the latter does not guarantee the resulting solution would be globally optimal [7]. The work in [14] suggests that solutions can escape local minima if load undergoes random fluctuations. Our approach can be seen as providing an explicit and deterministic algorithm to search for global solutions by using dual variables at local solutions.

Consider a power system network where $n$ buses are connected by $m$ edges. For bus $i$, let $P^{G}_{i}$ and $Q^{G}_{i}$ denote the active and reactive output of the generated and $P^{D}_{i}$ and $Q^{D}_{i}$ denote the active and reactive load. We use $P^{f}_{ij}$ and $Q^{f}_{ij}$ to denote the active and reactive power flowing from bus $i$ to bus $j$. Note that if buses $i$ and $j$ are not connected, then the flows are zero. We use $\theta_{ij}$ as a shorthand for $\theta_{i}-\theta_{j}$.

The OPF problem is to minimize the total of active power generation costs while satisfying a number of constraints:

	$\displaystyle\min_{\mathbf{V},\boldsymbol{\theta}}$	$\displaystyle\textstyle\sum_{i}c_{i}(P^{G}_{i})$		(1a)
	s.t.	$\displaystyle P^{G}_{i}=P^{D}_{i}+\textstyle\sum_{j=1}^{N}P^{f}_{ij}$		(1b)
		$\displaystyle Q^{G}_{i}=Q^{D}_{i}+\textstyle\sum_{j=1}^{N}Q^{f}_{ij}$		(1c)
		$\displaystyle P^{f}_{ij}=V_{i}^{2}g_{ij}-V_{i}V_{j}(g_{ij}\cos(\theta_{ij})-b_{ij}\sin(\theta_{ij}))$		(1d)
		$\displaystyle Q^{f}_{ij}=V_{i}^{2}\hat{b}_{ij}-V_{i}V_{j}(b_{ij}\cos(\theta_{ij})+g_{ij}\sin(\theta_{ij}))$		(1e)
		$\displaystyle\underline{V}_{i}\leq V_{i}\leq\bar{V}_{i}$		(1f)
		$\displaystyle\underline{P}^{G}_{i}\leq P^{G}_{i}\leq\bar{P}^{G}_{i}$		(1g)
		$\displaystyle\underline{Q}^{G}_{i}\leq Q^{G}_{i}\leq\bar{Q}^{G}_{i}$		(1h)
		$\displaystyle(P^{f}_{ij})^{2}+(Q^{f}_{ij})^{2}\leq(S_{ij}^{\max})^{2}$		(1i)

where $\hat{b}_{ij}=b_{ij}-0.5b_{ij}^{C}$ and $b_{ij}^{C}$ is the line charging susceptance. The constraints (1b) and (1c) enforce power balance, (1d) and (1e) are the AC power flow equations, (1f) limits the bus voltage magnitudes, (1g) and (1h) represent the active and reactive limits and (1i) are the line flow limits. The cost at bus $i$ is $c_{i}(\cdot)$ and can be linear, quadratic or other functions.

A pair of $(\mathbf{V},\boldsymbol{\theta})$ that minimizes the objective cost among all feasible solutions to (1) is called a global solution. In practice, a nonlinear programming (NLP) solver may only return a local solution, which is a feasible point that satisfy some local optimality conditions (e.g., KKT). Since a local solution is not necessarily global, in the following subsection, we propose an iterative approach to find global solutions by alternatively solving (1) and its partial Lagrangian. Any existing NLP can be used, and and we use IPOPT in this paper.

Our algorithm starts with a call to a NLP solver with an initial guess, denoted by $\boldsymbol{\theta}_{\text{init}}$, $\mathbf{V}_{\text{init}}$. For example, this can be the standard flat start. Then we assume the solver returns a feasible solution. At this solution, we record the dual variables associated with the power balance equations (1b) and (1c), denoted as $\bar{\boldsymbol{\mu}}^{P}$ and $\bar{\boldsymbol{\mu}}^{Q}$. Using these dual variables, we form the following partial Lagrangian by dualizing the power balance equations:

	$\displaystyle\mathcal{L}(\mathbf{V},\boldsymbol{\theta},\boldsymbol{\mu}^{P},\boldsymbol{\mu}^{Q})=$	$\displaystyle\textstyle\sum_{i}c_{i}P^{G}_{i}+\textstyle\sum_{i}\mu_{i}^{P}(P^{D}_{i}+\textstyle\sum_{j=1}^{N}P^{f}_{ij}-P^{G}_{i})$
		$\displaystyle+\textstyle\sum_{i}\mu_{i}^{Q}(Q^{D}_{i}+\sum_{j=1}^{N}Q^{f}_{ij}-Q^{G}_{i}).$

We then minimize the partial Lagrangian by solving

	$\displaystyle\min_{\mathbf{V},\boldsymbol{\theta}}\;$	$\displaystyle\mathcal{L}(\mathbf{V},\boldsymbol{\theta},\boldsymbol{\mu}^{P},\boldsymbol{\mu}^{Q})$		(3)
	s.t.	$\displaystyle~{}(\ref{PEq})-(\ref{flimits}).$

We solve the problem in (3) at the same initial point $(\mathbf{V}_{\text{init}},\boldsymbol{\theta}_{\text{init}})$ that was used to solve the original primal problem in (1). Denote this solution to (3) by $(\bar{\mathbf{V}},\bar{\boldsymbol{\theta}})$. Then we start the solver again to solve (1) but with the initial point $(\bar{\mathbf{V}},\bar{\boldsymbol{\theta}})$. This process can be repeated until the solutions stop changing or up to a predefined number of iterations.

It tures out that the initial point $(\bar{\mathbf{V}},\bar{\boldsymbol{\theta}})$ that is given by solving the partial Lagrangian is often a much better starting point than the original choice of $(\mathbf{V}_{\text{init}},\boldsymbol{\theta}_{\text{init}})$. Therefore, by repeating these steps, we can iteratively improve the quality of the starting point and get the global solution. The algorithm is summarized below as Algorithm 1. We illustrate the intuition behind this algorithm in the next section using a two bus network and present the numerical results in Section V.

Algorithm 1 is successful because of the geometry of the partial Lagrangian is “better” than the geometry of the original problem. We use a two-bus network as an illustrative example. For simplicity, we set both voltage magnitudes to 1 p.u. and optimize over the angle. Suppose bus 1 is a generator and the reference (slack) bus with linear cost at $1/MW, and bus 2 is the load bus with angle $-\theta$. The line admittance is $g-jb$. Given a load of $l$ at bus 2, the ACOPF is:

	$\displaystyle\min_{\theta}\;$	$\displaystyle g-g\cos(\theta)+b\sin(\theta)$		(4a)
	s.t.	$\displaystyle l+g-g\cos(\theta)-b\sin(\theta)=0.$		(4b)

For a feasible load $l$, there are two solutions to (4b). To make visualize how an nonlinear solver would approach (4), we convert it to a penalized unconstrained problem [15]:

	$\displaystyle\mathcal{L}_{\rho}=$	$\displaystyle g-g\cos(\theta)+b\sin(\theta)$		(5)
		$\displaystyle+\rho/2(l+g-g\cos(\theta)-b\sin(\theta))^{2},$

where $\rho$ is a (large) penalty parameter. The function $\mathcal{L}_{\rho}$ is plotted in Fig. 2 and there are two local minima, with the left one being global. If a NLP is initialized in the right valley, the suboptimal solution would be found.

Now consider the partial Lagrangian for (4), given by

	$\displaystyle\mathcal{L}_{\mu}=$	$\displaystyle g-g\cos(\theta)+b\sin(\theta)$		(6)
		$\displaystyle+\mu(l+g-g\cos(\theta)-b\sin(\theta)),$

where $\mu$ is the multiplier corresponding to the equality constraint (4b) at a local solution. The blue line in Fig. 2 plots $\mathcal{L}_{\mu}$ when dual variable is computed at the right (suboptimal) local solution. In contrast to $\mathcal{L}_{\rho}$, $\mathcal{L}_{\mu}$ has a unique minimum that falls into the left valley, and can be obtained by minimizing $\mathcal{L}_{\mu}$ from any initialization. Therefore, if we use the minimum of $\mathcal{L}_{\mu}$ as an initialization point for the nonlinear solver, we would reach the global optimal solution.

The above example illustrates that the Lagrangian has a more advantageous geometry for optimization than the original primal problems. This is not surprising, since equality constraints tend to complicate the optimization landscape [2]. The key insight is to observe that it suffices to use the dual variables at suboptimal local solutions to form the Lagrangian.

In this section we report the simulation results to validate the effectiveness of our algorithm. The NLP solver used is IPOPT [11] and the convergence tolerance is set to $0.0001$. We assume IPOPT returns a feasible solution, which may not be a global optimum. We test our algorithm on different sizes of IEEE networks. For the 9-bus and 22-bus networks, we utilize the local solutions found in [7] as starting points to launch IPOPT. Specifically, there are 4 local solutions in the 9-bus network, and 2 local solutions in the 22-bus network. For each network, starting from different local solutions, our algorithm can achieve global solutions within one iteration.

For the 39-bus network, we take a set of $600$ initial points for $\boldsymbol{\theta}$ and $\mathbf{V}$. To find the global solution, we do an exhaustive search within the bounds of each variable. In Fig. 3, we illustrate the improvement of solution quality for the 39-bus network. The x-axis represents the number of iterations that Algorithm 1 is ran, and y-axis represents the percentage of globally optimal solutions after each iteration. Without Algorithm 1, less than half of the solutions are globally optimal. One application of Algorithm 1 increases the percentage of globally optimal solutions to $98\%$. After two iterations, only four cases do not reach global optimas. All solutions are globally optimal after three iterations.

In this paper, we propose a simple algorithm to find globally optimal solutions to ACOPF problems iteratively. First, we solve the ACOPF problem using an existing nonlinear solver. From the solution and its associated dual variables, we construct a partial Lagrangian. Optimizing this partial Lagrangian leads to a new solution. With this solution as an initial point, we again call the solver for the ACOPF problem. By repeating these steps, we can iteratively improve the solution quality, escaping from local solutions to achieve global optimums. We illustrate the intuition behind our algorithm using a two-bus network, which shows that the partial Lagrangian has a flatter optimization landscape compared to the original primal problem. We validate the effectiveness of our algorithm on standard 9-bus, 22-bus and 39-bus networks. Regardless of the initial points, our algorithm always finds the global optimum within at most three iterations.