1 Introduction

Power flow is the underlining problem for power system analysis. Integration of intermittent renewable energy resources and possible network contingencies further highlight the merit of providing an efficient framework to solve this nonlinear problem. The power flow problem is formulated as a set of nonlinear equations. Several iterative approaches including Gauss–Seidel (GS) and Newton–Raphson (NR) were adopted to solve this problem. However, the convergence and stability of these approaches could not be guaranteed. NR is the most popular approach to solve this system of equations as it provides a better convergence rate compared to other techniques. Employing NR to solve the power flow problem leads to the following scenarios: ① a unique solution exists that can be found regardless of the initial point; ② multiple solutions exist and one of the solutions is returned based on the initial point; ③ no solution exists; ④ unique or multiple solutions exist while no solution can be procured because of improper initial point that renders a singular Jacobean matrix in the iterative process. The latter scenario presents an ill-conditioned power flow problem.

Several research works were focused to solve the power flow problem. Lower-upper (LU) factorization is a direct approach to solve NR that is computationally expensive and impractical for large-scale problems. Krylov subspace method [1] is a nonstationary iterative method that converges in at most NB iterations, where NB is the number of buses in the network. The Newton-generalized minimal residual (GMRES) method is one of the Krylov methods utilized to solve the power flow problem. The adaptive preconditioning schemes to update the preconditioners for the linearized equations of the next iteration and a flexible inner-outer preconditioned GMRES were proposed in [2] and [3], respectively, to improve the convergence of the Krylov methods. Another approach is to utilize incomplete LU factorization as a preconditioner for the Krylov–Newton methods as presented in [4]. The continuous Newton’s method proposed in [5] formulates the power flow problem as a set of ordinary differential equations.

Finding a proper preconditioning for the Krylov methods to procure a feasible solution for the power flow problem is challenging. Although many research works addressed the feasibility of the power flow problem [6,7,8], the proposed approaches for finding a feasible solution or providing a certificate for the infeasibility of the problem are not yet effective for ill-conditioned power flow problems. Iterative approaches are incapable of handling the ill-conditioned power flow problems [5, 9].

The analytical approaches (e.g. convex relaxation [10]) can address the singularity issue with the ill-conditioned power flow problems. However, improving the scalability of these approaches is a challenging task [11]. The solution rendered by leveraging the semi-definite programming (SDP) relaxation may not be a tight one that is feasible for the non-convex power flow problem. Therefore, more cutting planes may be needed to obtain a feasible solution with reasonable computational time. A high order of the Lasserre hierarchy could be used to find a feasible solution for the relaxed power flow problem [12]. However, the large computation burden of leveraging higher orders of the Lasserre hierarchy restricts the application of this approach for small-scale power networks. Nevertheless, exploiting the sparsity in power flow equations facilitates large-scale application of the SDP relaxation i.e. first order of Lasserre hierarchy [13].

In this paper, an efficient convex relaxation approach is presented to find a feasible solution for the power flow problem. The convex relaxation is tightened using appropriate perturbation functions along with a network reconfiguration scheme. The solution rendered by the tight convex relaxation problem is feasible for the original non-convex problem. The contributions of this paper are as follows:

  1. 1)

    A convex relaxation formulation for the power flow problem is presented, where a set of lifting variables is defined for the nonlinear terms in the power flow formulation. The introduced lifting variables are utilized in the lowest order of moment relaxation to find a solution to a relaxed power flow problem.

  2. 2)

    By exploiting the sparsity in the power network, the size of moment relaxation matrices employed to formulate the relaxed problem is substantially reduced. This reduction in the size of the problem facilitates the scalability of the presented convex relaxation approach.

  3. 3)

    A perturbation scheme is presented to improve the tightness of the presented relaxed problem. Once the original power flow problem is formulated as an optimization problem, the objective is zero. The presented relaxation for this optimization problem might not be tight with the lowest order of the moment relaxation. Presenting a perturbation function to tighten to relaxation does not negatively impact the procured solution to the relaxed problem as the objective of the relaxed optimization problem is zero.

  4. 4)

    A network reconfiguration scheme is implemented to improve the tightness of the presented relaxation. The network reconfiguration eliminates the zero-injection buses within the network. The procured voltages with the original topology are the same as those with the reconfigured topology. It is shown that removing the zero-injection buses along with perturbation improves the tightness of the presented sparse convex relaxation for the power flow problem.

  5. 5)

    A specific tightness measure is introduced to evaluate the effectiveness of the presented solution method to obtain a feasible solution to the original power flow problem. Moreover, a recovery process is presented to procure the solution to the original power flow problem from the solution provided by the relaxed problem.

A set of work addressed the optimal power flow problem in the literature [14,15,16,17,18,19], using SDP and SOCP relaxations, however, the objective of optimal power flow problem is different from power flow feasibility problem. The focus of the optimal power flow problem is to find the dispatch of generation units with minimum cost. Employing perturbation for the optimal power flow problem will undermine the optimality of the procured solution, while it is not an issue for the power flow problem.

This paper is organized as follows: Section 2 presents the power flow problem formulation. Section 3 presents a solution methodology to evaluate the existence of a feasible solution and find a solution to the power flow problem. Section 4 presents the numerical results to show the effectiveness of the proposed approach, and the conclusions are presented in Section 5.

2 Problem formulation

The problem formulation for the power flow problem is presented in (1). If the power network has NB buses, there are 2NB known parameters and 2NB unknown variables in the power flow equations. The magnitude and angle of the voltage for the slack/reference bus, the voltage magnitude and real power injection at voltage-controlled buses, and real and reactive power injections at load buses are the known parameters for the power flow problem. The unknown variables are the real and reactive power injections for the slack buses, voltage angle and reactive power injection for the voltage-controlled buses, and the voltage magnitude and voltage angle for the load buses. Real and reactive power injections are not usually enforced for the slack bus as they supposed to compensate for the real and reactive power mismatches in the network. The slack bus is considered as a reference bus to calculate the voltage angle for the load and voltage-controlled buses. The voltage for the slack/reference bus is enforced as shown in (1). The real power balance for the voltage-controlled and load buses is given in (2). The reactive power balance for the load buses is shown in (3). The voltage magnitudes for the voltage-controlled buses are enforced by (4). The reactive power generation limits for the generation units connected to voltage-controlled buses are shown in (5). If the reactive power generation of a unit reaches its limits, the reactive power generation is fixed to the limit and the corresponding bus becomes a load bus.

$$\left\{ {\begin{array}{*{20}l} {V_{i}^{d} = V_{ref}^{d} } \hfill & {i \in S} \hfill \\ {V_{i}^{q} = V_{ref}^{q} } \hfill & {i \in S} \hfill \\ \end{array} } \right.$$
(1)
$$\begin{aligned} P_{i}^{G\prime } - P_{i}^{D} = & \,\sum\limits_{j = 1}^{{N_{B} }} {\left[ {G_{ij} \left( {V_{i}^{d} V_{j}^{d} + V_{i}^{q} V_{j}^{q} } \right)} \right.} \\ & \left. { - B_{ij} \left( {V_{i}^{d} V_{j}^{q} - V_{i}^{q} V_{j}^{d} } \right)} \right]\quad \forall i \in PV\cup PQ \\ \end{aligned}$$
(2)
$$\begin{aligned} Q_{i}^{D} = & \,\sum\limits_{j = 1}^{{N_{B} }} {\left[ {B_{ij} \left( {V_{i}^{d} V_{j}^{d} + V_{i}^{q} V_{j}^{q} } \right)} \right.} \\ & \left. { +\, G_{ij} \left( {V_{i}^{d} V_{j}^{q} - V_{i}^{q} V_{j}^{d} } \right)} \right]\quad \forall i \in PQ \\ \end{aligned}$$
(3)
$$\sqrt {\left( {V_{i}^{d} } \right)^{2} + \left( {V_{i}^{q} } \right)^{2} } = \left| {V_{i}^{\prime } } \right|\quad \forall i \in PV$$
(4)
$$\begin{aligned} Q_{i,\hbox{min} }^{G} \le & \,Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} \left( {V_{i}^{d} V_{j}^{d} + V_{i}^{q} V_{j}^{q} } \right)} \right.} \\ & \left. { - G_{ij} \left( {V_{i}^{d} V_{j}^{q} - V_{i}^{q} V_{j}^{d} } \right)} \right] \le Q_{i,\hbox{max} }^{G} \quad \forall i \in PV \\ \end{aligned}$$
(5)

where \(V_{i}^{d}\) and \(V_{i}^{q}\) are the real and imaginary parts of the voltage phasor at bus i, respectively; \(V_{ref}^{d}\) and \(V_{ref}^{q}\) are the real and imaginary parts of the voltage phasor at reference bus, respectively; \(\left|V_{i}^{\prime}\right|\) is the voltage magnitude at voltage-controlled bus i; \(P_{i}^{G\prime }\) is the real power generation at bus i; \(P_{i}^{D}\) is the real power demand at bus i; Bij is the element of the susceptance matrix; Gij is the element of the conductance matrix; PQ is the set of load buses; PV is the set of voltage-controlled buses; \(Q_{i}^{D}\) is the reactive power demand at bus i; \(Q_{i,\hbox{min} }^{G}\) and \(Q_{i,\hbox{max} }^{G}\) are the minimum and maximum reactive power generations at bus i.

Solving (1)–(5) could be challenging under certain circumstances as discussed earlier and iterative approaches may fail to find a feasible solution. Thus, a non-iterative solution methodology is presented in the next section which tried to address this challenge.

3 Solution methodology

To find the solution for (1)–(5), it is reformulated as a convex optimization problem that could render a feasible solution for the power flow problem in polynomial time. Particularly, a sparse SDP relaxation for the problem presented in (1)–(5) is formulated with a suggested perturbation function. A topology reconfiguration scheme is proposed to improve the tightness of the presented convex relaxation. The details of the proposed solution methodology are presented in the following subsections.

3.1 Convex relaxation

A promising approach for solving polynomial optimization problems is the Lasserre hierarchy of moment relaxation, where SDP relaxation represents the first order of relaxation in this hierarchy [12]. Theoretically, with the increase in the order of the relaxation to infinity, the relaxation becomes tighter and the procured solution is feasible for the original polynomial optimization problem. However, the large computation burden of introducing higher orders of moment relaxation makes this approach impractical for large-scale applications. The computation complexity of first order of moment relaxation (i.e. SDP relaxation) is O(n3), where n is the number of monomials which is twice as the number of buses for the relaxed power flow problem. Exploiting the network sparsity mitigates the computation burden. Several sparse moment matrices associated with each maximal clique are defined. A clique, by definition, is a set of nodes within a graph that are all adjacent to each other. The maximal clique is a clique that its set of nodes is not a subset of any other clique. Here, the problem is reformulated as a first-order moment relaxation problem (SDP problem). The nonlinear terms in (1)–(5) are represented by respective lifting variables in the SDP relaxation matrix as formulated for the relaxed problem in (6)–(19). If all sparse SDP relaxation matrices are near-rank-1, the presented relaxation in (6)–(19) is tight and a feasible solution for the power flow problem in (1)–(5) is procured.

$$\mathop {\hbox{min} }\limits_{{X_{c} }} \sum\limits_{i} {\left( {\gamma_{{\left( {V_{i}^{d} } \right)^{2} }} + \gamma_{{\left( {V_{i}^{q} } \right)^{2} }} - 2V_{i}^{d} } \right)}$$
(6)

s.t.

$$\left\{ {\begin{array}{*{20}l} {\gamma_{{\left( {V_{i}^{d} } \right)^{2} }} = \left( {V_{ref}^{d} } \right)^{2} } \hfill & {\forall i \in S} \hfill \\ {\gamma_{{\left( {V_{i}^{q} } \right)^{2} }} = \left( {V_{ref}^{q} } \right)^{2} } \hfill & {\forall i \in S} \hfill \\ \end{array} } \right.$$
(7)
$$\begin{aligned} P_{i,\hbox{min} }^{G} \le & \,P_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ {G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.} \\ & \left. { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right] \le P_{i,\hbox{max} }^{G} \quad \forall i \in S \\ \end{aligned}$$
(8)
$$\begin{aligned} Q_{i,\hbox{min} }^{G} \le & Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.} \\ & \left. { - G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right] \le Q_{i,\hbox{max} }^{G} \quad \forall i \in S \\ \end{aligned}$$
(9)
$$\begin{aligned} P_{i}^{G\prime } - P_{i}^{D} = & \sum\limits_{j = 1}^{{N_{B} }} {\left[ {G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.} \\ & \left. { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right]\quad \forall i \in PV\cup PQ \\ \end{aligned}$$
(10)
$$\begin{aligned}Q_{i}^{D} =& \sum\limits_{j = 1}^{{N_{B} }} {\left[ {B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.}\\& \left. { + G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right]\quad \forall i \in PQ\end{aligned}$$
(11)
$$V_{i,\hbox{min} }^{2} \le \gamma_{{\left( {V_{i}^{d} } \right)^{2} }} + \gamma_{{\left( {V_{i}^{q} } \right)^{2} }} \le V_{i,\hbox{max} }^{2} \quad \forall i \in PQ$$
(12)
$$\begin{aligned} & \left| {V^{\prime}_{i} } \right|^{2} - (u_{i,\hbox{max} }^{Q} + u_{i,\hbox{min} }^{Q} )(V_{i,\hbox{max} }^{2} - V_{i,\hbox{min} }^{2} ) \\ & \quad \le \gamma_{{\left( {V_{i}^{d} } \right)^{2} }} + \gamma_{{\left( {V_{i}^{q} } \right)^{2} }} \quad \forall i \in PV \\ \end{aligned}$$
(13)
$$\begin{aligned} & \gamma_{{\left( {V_{i}^{d} } \right)^{2} }} + \gamma_{{\left( {V_{i}^{q} } \right)^{2} }} \le \left| {V^{\prime}_{i} } \right|^{2} \\ & \quad + (u_{i,\hbox{min} }^{Q} + u_{i,\hbox{max} }^{Q} )(V_{i,\hbox{max} }^{2} - V_{i,\hbox{min} }^{2} )\quad \forall i \in PV \\ \end{aligned}$$
(14)
$$\begin{aligned} & Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.} \left. { - G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right] \\ & \quad - \varepsilon u_{i,\hbox{max} }^{Q} \le Q_{i,\hbox{max} }^{G} - \varepsilon \quad \forall i \in PV \\ \end{aligned}$$
(15)
$$\begin{aligned} & Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} (\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} ) - G_{ij} (\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} )} \right]} \\ & \quad + (2 + \varepsilon )Q_{i,\hbox{max} }^{G} u_{i,\hbox{max} }^{Q} \ge ( - 1 - \varepsilon )Q_{i,\hbox{max} }^{G} \quad \forall i \in PV \\ \end{aligned}$$
(16)
$$\begin{aligned} & Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right)} \right.} \left. { - G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right] \\ & \quad + \varepsilon u_{i,\hbox{min} }^{Q} \ge Q_{i,\hbox{min} }^{G} + \varepsilon \quad \forall i \in PV \\ \end{aligned}$$
(17)
$$\begin{aligned} & Q_{i}^{D} + \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} (\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} ) - G_{ij} (\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} )} \right]} \\ & \quad - (2 + \varepsilon )Q_{i,\hbox{min} }^{G} u_{i,\hbox{min} }^{Q} \le ( - 1 - \varepsilon )Q_{i,\hbox{min} }^{G} \quad \forall i \in PV \\ \end{aligned}$$
(18)
$$\varvec{X}_{c} = \left[ {\begin{array}{*{20}c} 1 & {V_{i}^{d} } & \cdots & {V_{\left| c \right|}^{d} } & {V_{i}^{q} } & \cdots & {V_{\left| c \right|}^{q} } \\ {V_{i}^{d} } & {\gamma_{{( {V_{i}^{d} } )^{2} }} } & \cdots & {\gamma_{{V_{i}^{d} V_{\left| c \right|}^{d} }} } & {\gamma_{{V_{i}^{d} V_{i}^{q} }} } & \cdots & {\gamma_{{V_{i}^{d} V_{\left| c \right|}^{q} }} } \\ \vdots & \vdots & {} & \vdots & \vdots & {} & \vdots \\ {V_{\left| c \right|}^{d} } & {\gamma_{{V_{\left| c \right|}^{d} V_{i}^{d} }} } & \cdots & {\gamma_{{( {V_{\left| c \right|}^{d} } )^{2} }} } & {\gamma_{{V_{i}^{d} V_{i}^{q} }} } & \cdots & {\gamma_{{V_{\left| c \right|}^{d} V_{\left| c \right|}^{q} }} } \\ {V_{i}^{q} } & {\gamma_{{V_{i}^{q} V_{i}^{d} }} } & \cdots & {\gamma_{{V_{i}^{q} V_{\left| c \right|}^{d} }} } & {\gamma_{{( {V_{i}^{q} } )^{2} }} \gamma_{{( {V_{i}^{q} } )^{2} }} } & \cdots & {\gamma_{{V_{i}^{q} V_{\left| c \right|}^{q} }} } \\ \vdots & \vdots & {} & \vdots & \vdots & {} & \vdots \\ {V_{\left| c \right|}^{q} } & {\gamma_{{V_{\left| c \right|}^{q} V_{i}^{d} }} } & \cdots & {\gamma_{{V_{\left| c \right|}^{q} V_{\left| c \right|}^{d} }} } & {\gamma_{{V_{\left| c \right|}^{q} V_{i}^{q} }} } & \cdots & {\gamma_{{( {V_{\left| c \right|}^{q} } )^{2} }} } \\ \end{array} } \right]\underline{ \succ }\, {\mathbf{0}}\quad \forall c$$
(19)

where \(P_{i,\hbox{min} }^{G}\) and \(P_{i,\hbox{max} }^{G}\) are the minimum and maximum real power generations at bus i; \(Q_{i}^{G\prime }\) is the reactive power generation at bus i; S is the set of slack/reference buses; \(V_{i,\hbox{min} }^{{}}\) and \(V_{i,\hbox{max} }^{{}}\) are the minimum and maximum voltage magnitudes at load bus i; c is the index for each maximal clique within the network graph; |c| is the number of buses within clique c; \(u_{i,\hbox{max} }^{Q}\) is the auxiliary binary variable, 1 if the reactive power generation at bus i reaches its maximum value, otherwise 0; \(u_{i,\hbox{min} }^{Q}\) is the auxiliary binary variable, 1 if the reactive power generation at bus i reaches its minimum value, otherwise 0; ε is an arbitrary small constant.

The objective of the perturbed convex relaxation is given in (6), where the lifting terms associated with the square of real and imaginary parts of voltage on each bus is employed. The choice of perturbation plays an important role in procuring a near-rank-1 solution. An SDP relaxation for the rank minimization problem is presented in [20], which may not render a feasible solution for the power flow problem. The choice of the perturbation is not unique; however, various functions may lead to various near-rank-1 solutions [21]. This choice depends on the system operator preferences to obtain a particular solution among multiple solutions that may exist for the power flow problem. The system operator does not need to know about the superiority of one of the perturbations over another. However, they might have various preferences and technical considerations to choose a perturbation function. Here, the perturbation matrix is employed in the objective function to determine a feasible solution in which the voltage magnitudes on the buses are close to 1 p.u.. This choice can be the last known voltage of the system, a desired voltage profile, 1 p.u. voltage for all buses, etc. The nonlinear terms in (1)–(5) are presented by their associated lifting variables in (2)–(19), and the SDP matrices associated with each maximal clique given in (19) contain the lifting variables. For example, \(V_{i}^{q} V_{j}^{d}\) is a nonlinear term in (1)–(5) which is replaced by a lifting variable \(\gamma_{{V_{i}^{q} V_{j}^{d} }}\) as defined in the semi-definite matrix constraint (19).

The voltage for the slack/reference bus is enforced by (7). Although, the real power balance for the slack bus is ignored in (1)–(5), enforcing the generation capacity limits for the slack bus will ensure the feasibility of the solution procured by solving the relaxed problem. Enforcing these limits will avoid procuring a solution that is impractical. The power flow problem demonstrates the state of the system, where slack bus cannot provide real and reactive power beyond its generation capacity limits. Thus, the real and reactive power generation capacities of the unit connected to the slack bus is enforced by (8) and (9), respectively. The real power balance for the voltage-controlled and load buses is presented in (10). The reactive power balance for the load buses is shown in (11). The voltage limits for load buses are not usually considered for the power flow problem. However, the power flow problem might have multiple solutions, where some of them are low voltage solution vulnerable to voltage collapse [22]. Thus, to ensure system security and the technical feasibility of the procured solution, the voltage limits for the load buses is presented in (12). For the voltage-controlled buses, the voltage magnitude is enforced by (13) and (14). Once the reactive power of the generation units connected to a voltage-controlled bus reaches its limits, the voltage-controlled bus will transform into a load bus with a fixed reactive power and unknown voltage magnitude. This condition is captured by two auxiliary binary variables for each voltage-controlled bus to check if any of the upper and lower limits for the reactive power generation of the generation units are reached. The reactive power generation of the units connected to voltage-controlled buses is enforced by (15)–(18). Here once the reactive power generation reaches the upper or lower limits, the auxiliary binary variable becomes 1.

Considering the auxiliary binary variable for voltage-controlled buses, the convex relaxation problem is formulated as a mixed-integer semi-definite programming (MISDP) problem. The solver developed for solving the MISDP problem is not very efficient for solving large-scale problems. There are two possible separate sets of actions to address this challenge. The first way is to utilize a branch and bound approach similar to those adopted in mixed-integer linear programming and mixed-integer quadratic programming problems [23]. To solve the MISDP problem, a branch and bound algorithm can be employed. The sub-problem in each node of the branch and bound algorithm is an SDP problem with sparse SDP relaxation matrices. Although employing the branch and bound algorithm is promising for small to medium size networks, its performance will be deteriorated with the increase in the number of voltage-controlled buses in the network. The alternative technique to address this challenge is to relax the binary variable, u, as a continuous variable, ur, as given in (20). Then enforce a non-convex constraint to ensure the continuous variable takes values 0 and 1 as given in (21). Then, employing a convex relaxation approach to convexify the non-convex feasibility region. To tighten such relaxation regularization linearization technique (RLT) and valid constraints [24] are leveraged. As the convex relaxation problem formulated in the sparse form, relaxing the binary variables would be a better choice for large-scale applications.

$$0 \le u^{r} \le 1$$
(20)
$$\left( {u^{r} } \right)^{2} = u^{r}$$
(21)

3.2 Tightness measure

The solution procured from the problem presented in (6)–(19) is feasible for the original power flow problem when the convex relaxation is tight, and the rank of all SDP relaxation matrices is near one. The rank of a matrix is equal to the number of its non-zero eigenvalues. The rank of a matrix is exactly one when it has only one non-zero eigenvalue. The measure for the tightness of the solution procured from (6)–(19) is the ratio of two largest eigenvalues of the procured SPD relaxation matrices as shown in (22), where TRc is the tightness measure for each clique c and \(\lambda_{\left| c \right|}^{c}\) is the eigenvalue of the SDP matrix associated with maximal clique c.

$$TR_{c} = \lg \left( {\frac{{\lambda_{\left| c \right|}^{c} }}{{\lambda_{\left| c \right| - 1}^{c} }}} \right)$$
(22)

A large ratio in (22) indicates that the second eigenvalue is very small compared to the first eigenvalue and can be neglected. Thus, if the ratio is a large number, the rank of the SDP matrix is near-one and the relaxation is tight. Alternatively, the gap between the derived solution and the original solution can be procured using the reciprocal of the presented tightness measure.

3.3 Recovering solution to original power flow problem

To recover the solution to the original power flow problem, the matrix of lifting variables procured as a solution to the relaxed problem is utilized. The voltages for the original problem can be procured using the Cholesky decomposition of the sparse matrices of lifting variables [25]. The vector of real and imaginary parts of voltages within a clique is procured by (23). If the rank of the sparse SDP matrix is near one, the largest eigenvalue of each sparse matrix will dominate the other eigenvalues. Thus, to recover the solution to the original problem within a maximal clique, the dominant eigenvalue is multiplied by the associated eigenvectors. The procured vector is the vector of voltages for the maximal clique which is the solution to the original power flow problem in (1)–(5).

$$\begin{aligned} & \left[ {1,V_{i}^{d} , \ldots ,V_{\left| c \right|}^{d} ,V_{i}^{q} , \ldots ,V_{\left| c \right|}^{q} } \right]^{\text{T}} \\ & \quad \cdot \left[ {1,V_{i}^{d} , \ldots ,V_{\left| c \right|}^{d} ,V_{i}^{q} , \ldots ,V_{\left| c \right|}^{q} } \right] = \varvec{q}_{\left| c \right|}^{c} \lambda_{\left| c \right|}^{c} \mathop {\left( {\varvec{q}_{\left| c \right|}^{c} } \right)}\nolimits^{\text{T}} \\ \end{aligned}$$
(23)

where \(\varvec{q}_{\left| c \right|}^{c}\) is the eigenvector of the SDP matrix associated with maximal clique c.

3.4 Network reconfiguration

Formulating the convex relaxation for the power flow problem in sparse form along with the proposed perturbation would improve the tightness of the proposed relaxation to find a feasible solution. However, to further improve the tightness, a network reconfiguration is proposed to avoid employing the higher order moment relaxation. A number of transshipment buses exist in the power network to facilitate the connection between the generation and load buses while they have zero power injection. In many occasions, the outcome of the relaxed problem is not near-rank-1 if the maximal cliques have load buses with zero real and reactive power injections. The reason is that the zero injection buses increase the degree of freedom for the set of power flow equations by introducing (24) and (25). This leads to a solution with a rank higher than one for the SDP matrix corresponding to the maximal clique containing these buses.

$$0 = \sum\limits_{j = 1}^{{N_{B} }} {\left[ {G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right) - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right]}$$
(24)
$$0 = \sum\limits_{j = 1}^{{N_{B} }} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{d} }} + \gamma_{{V_{i}^{q} V_{j}^{q} }} } \right) - G_{ij} \left( {\gamma_{{V_{i}^{d} V_{j}^{q} }} - \gamma_{{V_{i}^{q} V_{j}^{d} }} } \right)} \right]}$$
(25)

To tackle this challenge, a network reconfiguration is proposed to eliminate the load buses with zero injection from the network. The procured topology is equivalent to the original network topology.

The following cases in which the degree of the load bus with zero power injection in the graph associated with the power network is equal to 1, 2, and 3 are considered. If the degree of the zero-injection load bus j is 1, as shown in Fig. 1, the bus can be removed from the network. The flow of the line connected to this bus is zero, and the voltage of this bus is equal to the bus connected to it. If the set of zero-injection buses with connectivity degree of 1 is Ω1 and \(\left\| \cdot\right\|_{0}\) indicates the number of nonzero elements, the voltage of the zero-injection bus can be further recovered as a function of the adjacent bus voltage using (26) and (27), where \(\gamma_{{V_{i}^{d} }}\) and \(\gamma_{{V_{i}^{q} }}\) are procured from the solution of the reconfigured network.

$$G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right) = 0$$
(26)
$$- B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right) = 0\,\,\,$$
(27)

where \(j \in \varOmega_{1} ; \, G_{ij} \ne 0; \, B_{ij} \ne 0; \, \left\| {G_{ij} } \right\|_{0} = \left\| {B_{ij} } \right\|_{0} = 1\).

Fig. 1
figure 1

Zero-injection load bus with degree 1

Full size image

If the degree of the zero-injection load bus j is 2, as shown in Fig. 2a, the load bus with zero injection can be removed from the network while the two lines connected to this bus will merge into a single line in the reconfigured topology, as shown in Fig. 2b. Here, impedance zik= zij+ zjk and the Y bus of the configured network is further adjusted. The voltage of this bus can be further recovered as a function of the adjacent bus voltages using (28) and (29), where Ω2 is the set of zero-injection load buses with connectivity degree of 2. Here, two unknowns i.e. \(V_{j}^{d}\) and \(V_{j}^{q}\), could be found once \(\gamma_{{V_{i}^{d} }}\) and \(\gamma_{{V_{i}^{q} }}\) are procured from the power flow solution of the reconfigured network.

$$\sum\limits_{i} {\left[ {G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right)} \right]} = 0\,\,\,$$
(28)
$$\sum\limits_{i} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right)} \right]} = 0\,\,$$
(29)

where \(j \in \varOmega_{2} ; \, G_{ij} \ne 0; \, B_{ij} \ne 0; \, \left\| {G_{ij} } \right\|_{0} = \left\| {B_{ij} } \right\|_{0} = 2\).

Fig. 2
figure 2

Zero-injection load bus with degree 2 and its reconfigured network

Full size image

If the degree of the zero-injection load bus j is 3, as shown in Fig. 3a, the load bus with zero injection can be removed from the network by changing the network topology to its equivalent shown in Fig. 3b. The branch impedances are procured using Y–△ conversion shown in (30) [26] and Y bus of the reconfigured network is constructed accordingly.

$$\left\{ {\begin{array}{*{20}l} {Z_{im} = Z_{ij} Z_{jm} (Z_{ij} + Z_{jk} + Z_{jm} )^{ - 1} } \\ {Z_{ik} = Z_{ij} Z_{jk} (Z_{ij} + Z_{jk} + Z_{jm} )^{ - 1} } \\{Z_{km} = Z_{kj} Z_{jm} (Z_{ij} + Z_{jk} + Z_{jm} )^{ - 1} } \end{array}} \right.$$
(30)
Fig. 3
figure 3

Zero-injection load bus with degree 3 and its reconfigured network

Full size image

The voltage of the removed buses can be recovered as a function of the adjacent bus voltages using the (31) and (32), where Ω3 is the set of zero-injection load buses with connectivity degree of 3. Here, two unknowns, i.e. \(V_{j}^{d}\) and \(V_{j}^{q}\), could be found once \(\gamma_{{V_{i}^{d} }}\) and \(\gamma_{{V_{i}^{q} }}\) are procured from the solution of the reconfigured network.

$$\sum\limits_{i} {\left[ {G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right)} \right]} = 0$$
(31)
$$\sum\limits_{i} {\left[ { - B_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{d} + \gamma_{{V_{i}^{q} }} V_{j}^{q} } \right) + G_{ij} \left( {\gamma_{{V_{i}^{d} }} V_{j}^{q} - \gamma_{{V_{i}^{q} }} V_{j}^{d} } \right)} \right]} = 0\,\,$$
(32)

where \(j \in \varOmega_{3} ; \, G_{ij} \ne 0; \, B_{ij} \ne 0; \, \left\| {G_{ij} } \right\|_{0} = \left\| {B_{ij} } \right\|_{0} = 3\).

4 Numerical results

To illustrate the effectiveness of the proposed methodology, several case studies are presented. The presented problem formulation is solved using MOSEK [27]. The tightness of the solutions to the power flow problem is compared in the following cases: ① Case 1, relaxation without perturbation and network reconfiguration; ② Case 2, relaxation with perturbation but without network reconfiguration; ③ Case 3, relaxation without perturbation but with network reconfiguration; ④ Case 4, relaxation with perturbation and network reconfiguration.

The comparison among the results for Case 2 and Case 3 with those for Case 4 presents the impact of perturbation and network reconfiguration individually.

4.1 Ill-conditioned 13-bus system

The power flow for this network is addressed in [9] as an ill-conditioned problem. An inappropriate choice of initial point leads to the divergence of power flow problem. The numbers of maximal cliques are 13, 13, 11, and 11 for Case 1, Case 2, Case 3, and Case 4, respectively. The maximal cliques within the network before and after network reconfiguration are respectively shown in Figs. 4 and 5 with dotted lines. Here, buses 2, 3, 4, 7, and 13 are the zero-injection buses that are removed in the reconfigured topology.

Fig. 4
figure 4

13-bus system before reconfiguration

Full size image
Fig. 5
figure 5

13-bus system after reconfiguration

Full size image

The tightness measure results of the sparse relaxation matrices given in Fig. 6 show the effectiveness of the proposed methodology to tackle the ill-conditioned power flow problem. The impact of perturbation is notable by comparing the tightness of various cases in Fig. 6. Here, the relaxation of Case 2 and Case 4 is much tighter than that for Case 1 and Case 3. It should be noted that cliques 1, 2, 3, 4, and 11 in Case 3 and Case 4 contain the load buses with zero injection that are removed from the network. Here, the essence of network reconfiguration and perturbation techniques to determine a feasible power flow solution is shown. Applying the perturbation in Case 2 leads to a tighter relaxation in clique 6, compared to that in Case 1. It is worth noting that the cliques in Case 2 and Case 4 are different as Case 4 shows the clique after network reconfiguration. Thus, pair-wise comparison among the cliques is not possible for these cases. However, the average tightness measure results for Case 4 are improved compared to those for Case 2. The voltage magnitude and angle for Case 2 and Case 4 are given in Fig. 7a and b, respectively. Comparison of the voltage profile procured by the presented approach with traditional approach is not available for this ill-conditioned test case [9]. The average solution time of for this case is 0.75 s.

Fig. 6
figure 6

Tightness measure results of Cases 1–4 in 13-bus system

Full size image
Fig. 7
figure 7

Profiles of voltage magnitude and angle for 13-bus system

Full size image

To illustrate the voltage recovery process from the solution of the relaxed problem, the maximal clique which includes bus 11 and bus 12 is selected. The procured solution for the sparse SDP matrix associated with this clique is presented in (33).

$$\left[ {\begin{array}{*{20}c} {1.0000} & {1.0896} & {1.0257} & {0.0387} & { - 0.0288} \\ {1.0896} & {1.1873} & {1.1177} & {0.0421} & { - 0.0314} \\ {1.0257} & {1.1177} & {1.0521} & {0.0397} & { - 0.0295} \\ {0.0387} & {0.0421} & {0.0397} & {0.0015} & { - 0.0011} \\ { - 0.0288} & { - 0.0314} & { - 0.0295} & { - 0.0011} & {0.0008} \\ \end{array} } \right]$$
(33)

The eigenvalues of the procured solution of the clique are 2.2617 × 10−9, 1.2243 × 10−8, 9.1753 × 10−8, 9.7194 × 10−8, and 3.2417. Here, there is only one eigenvalue which is much greater than zero. The tightness ratio for this clique defined in (22) is 7.52. Therefore, the presented solution to the relaxed problem is near-rank-1. The eigenvector associated with the largest eigenvalue is [− 0.5554, − 0.6051, − 0.5696, − 0.0214, 0.0159]. According to (23), the recovered voltages from the dominant eigenvalue and its associated eigenvector are the same as (33) up to seven-digit precision. The real and imaginary parts of the voltage of bus 11 are 1.0896 and 0.0387, respectively. The real and imaginary parts of the voltage of bus 12 are 1.0257 and − 0.0288, respectively. Converting the rectangular form of voltages to polar form returns \(1.0 90 3\angle 2.0 3 2 2^\circ\) and \(1.0 2 6 1\angle 1. 60 7 1{^\circ }\) for bus 11 and bus 12, respectively.

4.2 IEEE 30-bus system

The data provided in [28] is utilized in this case. The numbers of maximal cliques are 29, 29, 27, and 27 for Case 1, Case 2, Case 3, and Case 4, respectively. The tightness measure results for the sparse SDP matrices are given in Fig. 8. As shown in this figure, the perturbation in Case 4 leads to a tighter relaxation for the power flow problem compared to Case 1, Case 2, and Case 3. By comparing Case 1 with Case 3, and Case 2 with Case 4, it is shown that network reconfiguration provides a tighter measure for the SDP relaxation of the power flow problem. Here, the sparse SDP matrices associated with cliques 9, 25, and 26 are tighter compared to those cliques for Case 1. As shown in Fig. 9, comparing the voltage profiles of Case 4 and NR algorithm reveals that the network reconfiguration along with perturbation enable finding of a feasible solution for the power flow problem. The voltage profile of Case 4 is a solution to the original power flow problem. The average solution time for Case 4 is 0.81 s.

Fig. 8
figure 8

Tightness measure results of Cases 1–4 in IEEE 30-bus system

Full size image
Fig. 9
figure 9

Comparison of voltage profiles of Case 4 and NR algorithm in IEEE 30-bus system

Full size image

4.3 IEEE 57-bus system

The data provided in [28] is utilized in this case. The numbers of maximal cliques are 62, 62, 59, and 59 for Case 1, Case 2, Case 3, and Case 4, respectively. The tightness measure results of the sparse SDP matrices of various cases are given in Fig. 10. The effect of the perturbation on the tightness of the presented relaxation for the power flow problem is shown by comparing the tightness measure results for Case 4 with those for Case 3. By comparing Case 2 with Case 4, and Case 1 with Case 3, it is shown that network reconfiguration improves the tightness of SDP relaxation. Here, the voltage deviations are smaller in Case 4 compared to Case 3. Moreover, the voltage angles are closer to each other in Case 4 compared to Case 3. As shown in Fig. 11, comparing the voltage profiles of Case 4 and NR algorithm shows the merit of the proposed approach, where the voltage for the slack bus and voltage-controlled buses are very close in the two approaches. The average solution time for Case 4 is 0.85 s.

Fig. 10
figure 10

Tightness measure results of Cases 1–4 in IEEE 57-bus system

Full size image
Fig. 11
figure 11

Comparison of voltage profiles of Case 4 and NR algorithm in IEEE 57-bus system

Full size image

4.4 200-bus system

This case addresses a synthetic 200-bus network provided in [29]. The numbers of maximal cliques are 223, 223, 213, and 213 for Case 1, Case 2, Case 3, and Case 4, respectively. The tightness measure results of the sparse SDP matrices for Case 3 and Case 4 are given in Fig. 12. The average and standard deviations of voltage magnitude of the buses in Case 4 are 1.0109 and 0.0107 p.u., respectively. The average and standard deviations of voltage magnitude for the buses in Case 3 are 0.9794 and 0.0503 p.u., respectively. The average and standard deviations of voltage angle for Case 4 are 1.0974° and 1.1740°, respectively. As shown in Fig. 13, comparing the voltage profiles of Case 4 and NR algorithm demonstrates that the voltage magnitudes for the solution procured by the approach presented in this paper are closer to 1 p.u. with less deviation compared to the solution procured by the NR algorithm. The average solution time for Case 4 is 1.11 s.

Fig. 12
figure 12

Tightness measure results of Cases 1–4 in 200-bus system

Full size image
Fig. 13
figure 13

Comparison of voltage profiles of Case 4 and NR algorithm in 200-bus system

Full size image

4.5 2383-bus system

This case addresses the Polish power network with 2383 buses [28]. The number of maximal cliques is 2640 for Case 3 and Case 4. The tightness measure results of the sparse SDP matrices of Case 3 and Case 4 are given in Fig. 14. The effect of perturbation on the tightness of the relaxation is shown by comparing the tightness measure results in Case 4 with those in Case 3. This verifies the merit of proposed approach shown in Case 4. As shown in Fig. 15, comparing the voltage profiles of Case 4 and NR algorithm reveals that the solution procured by the NR algorithm violates the limits presented in the datasheet which should be between 0.94 and 1.12. The voltage profile of the solution procured by the approach presented in this paper remains within those limits. The voltage procured by the NR algorithm is a low-voltage solution that may threaten the security of the system as discussed in [22]. The average solution time for Case 4 is 11.63 s.

Fig. 14
figure 14

Tightness measure results of Cases 3 and 4 in 2383-bus system

Full size image
Fig. 15
figure 15

Comparison of voltage profiles of Case 4 and NR algorithm in 2383-bus system

Full size image

5 Conclusion

This paper presents a convex relaxation approach to determine the feasibility of the power flow and yields a solution for this problem. The presented approach is scalable for large-scale applications. The near-rank-1 relaxation is procured by using the presented perturbation and network reconfiguration techniques. Leveraging these techniques helps avoiding employment of higher orders of moment relaxation and improves the computation efficiency of the presented relaxation for large networks. The effectiveness of the proposed approach is illustrated for well-conditioned IEEE test cases as well as for an ill-conditioned case.