A Mixed-Integer Quadratic Formulation of the Phase-Balancing Problem in Residential Microgrids

Abstract

Phase balancing is a classical optimization problem in power distribution grids that involve phase swapping of the loads and generators to reduce power loss. The problem is a non-linear integer and, hence, it is usually solved using heuristic algorithms. This paper proposes a mathematical reformulation that transforms the phase-balancing problem in low-voltage distribution networks into a mixed-integer convex quadratic optimization model. To consider both conventional secondary feeders and microgrids, renewable energies and their subsequent stochastic nature are included in the model. The power flow equations are linearized, and the combinatorial part is represented using a Birkhoff polytope ${\mathcal{B}}_3$ that allows the selection of phase swapping in each node. The numerical experiments on the CIGRE low-voltage test system demonstrate the use of the proposed formulation.

1. Introduction

In recent years, electric distribution networks have become more intricate by connecting new technologies, especially in medium- and low-voltage levels [1,2]. These technologies include renewable energy sources (i.e., photovoltaic and wind farms), energy storage devices, and controllable and non-controllable loads [3,4]. A network that integrates several of these new technologies, with the possibility of operating as an individual unit under grid-connected or isolated modes, is known as a microgrid [5,6]. Microgrids may operate under unbalanced conditions due to single-phase residential loads and the connection of small single-phase solar-photovoltaic sources and batteries [7,8]. The main problem associated with unbalanced loads is the increment of the conductors’ energy losses due to the imbalance produced in the current magnitudes and the increase in current through the neutral. These characteristics in the network configuration lead to an interesting optimization problem known as optimal phase-balancing [9].

Models for phase-balancing have been developed in many applications, such as in aircraft electric distribution systems [10] and in power distribution grids with a high penetration of electric vehicles [11]. Due to its combinatorial nature, the phase-balancing problem is usually solved using heuristics [12] and meta-heuristics [13] as well as expert systems [14]. Modern approaches include the uncertainty associated with the load and generator [15].

While the phase–phase balance is a classical problem in power distribution networks, no previous scientific literature report exists that addresses this problem using integer-convex optimization. Mainly, the problem has been solved with metaheuristics due to its binary and non-linear nature [16]. Notwithstanding, we propose a mixed-integer quadratic programming (MI-QP) approximation, which has not been previously reported in the literature.

Unbalance in a network, it is important to note, can also be reduced using power electronics technologies. For example, an active filter with a four-leg converter can control the neutral current, and unbalance is reduced [17]. However, these devices imply an investment cost, while the solution of the optimization model proposed here only requires a suitable rearrangement of the loads. Nevertheless, both solutions are complementary and may be used together.

The problem of the optimal-phase balancing in electrical distribution networks and microgrids is complex for three main reasons, namely: (i) Power balance equations constitute a high-dimensional non-linear and non-convex set of constraints; (ii) the problem is combinatorial due to the binary variables, associated with the possible options for balancing loads and generators at each node [16,18]; and (iii) the problem is stochastic due to the high variation of the load and generation of renewable sources [15]. These complications in the solution of the phase-balancing problem motivate the development of the current research. It is possible to obtain a convex relaxation of it by using a linear representation of the power flow problem in three-phase networks [19]. For this representation, a linearization is applied to the relation between voltages and powers, allowing the rewriting of them as a set of complex affine functions [20]. We recur to the integer convex optimization theory regarding the binary variables that model the phase-balancing problem. The optimal solution of an integer-convex optimization problem can be guaranteed via the branch and bound method [21]. For each binary combination of variables, the internal problem is convex, implying that each internal problem has a global solution [22].

The remainder of the paper is organized as follows: Section 2 describes the main features of the problem; then, Section 3 presents the stationary state model of a microgrid and how the stochastic nature of the problem is managed. Section 4 shows the power flow linearization that convexifies the problem; after, Section 5 presents the Birkhoff polytope concepts, which allow an efficient geometric representation of the combinatorial part of the problem, with a reduced number of integer variables. Section 6 presents the numerical experiments performed and is followed by conclusions in Section 7 and the relevant references.

2. Problem Formulation

Microgrids are increasingly used in power distribution grids to integrate renewable energies [5]. They include single-phase components such as loads and solar-photovoltaic units. Therefore, microgrids are usually unbalanced, especially in low-voltage residential cases [7]. To reduce power loss, microgrid designers need to define the right placement of each single-phase unit (load or generation) using a process known as phase balancing [16]. Each three-phase node has six possible configurations as depicted in Figure 1. The problem involves defining each load or generator’s configuration to reduce the microgrid’s total losses. Therefore, since there are $6^n$ possible configurations, where n is the number of three-phase nodes, the problem is combinatorial. The problem is not trivial, even in small networks. For example, a microgrid with 10 nodes will have more than 60 million possible configurations.

A general representation of the problem is as follows:

$$\begin{array}{cc}\hfill min& \mathbb{E}\left(p_L\right)\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & f(v,M)=0\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & M_i\in \left\{M_1,M_2,M_3,M_4,M_5,M_6\right\}\hfill \end{array}$$

(3)

where $\mathbb{E}$ represents the expected value; $p_L$ is the active power loss; f is a vector function that represents the power flow equations; and M is the phase-swapping included in the set of possible configurations already depicted in Figure 1 (e.g., the configuration $M_i$ with $i=3$, i.e., $M_3$, implies that the original configuration ABC changes to BCA). Each of this model’s equations shows a different source of complexity: (i) the objective function (1) is stochastic; (ii) the equations of the power flow (2) are non-convex, and (iii) the configuration (3) leads to a combinatorial problem. This paper presents a mixed-integer quadratic formulation to this problem, which allows to find an optimal solution. In this way, by using the paradigm of disciplined convex programming, the non-linear integer stochastic problem is transformed into a tractable model [21,23].

3. Model of the Microgrid

A microgrid is represented by a connected graph with $3n$ nodes and a three-phase nodal admittance matrix $Y\in {\mathbb{C}}^{3n\times 3n}$. All the values are given per-unit under a well-defined basis. The set of nodes is divided into two subsets, namely: S for the slack node, and N for the rest of the nodes. Note that the slack node is defined in the three-phase domain, implying that it fixes the voltage output for each phase using positive sequence. Thus, the voltages in the slack node are given by (4),

$$V_S={(1,e^{-2\pi /3j},e^{2\pi /3j})}^{\top }.$$

(4)

The loads are represented by exponential models [24]. Therefore, the power flow equations take the following representation:

$${\left(\frac{s_k}{v_k}\right)}^{*}{∥v_k∥}^{\alpha }=\sum _{s=1}^3{y}_{ks}{v}_S+\sum _{m\in N}{y}_{km}{v}_m,\forall k\in N.$$

(5)

Note that (5) is non-linear and defined in the complex domain. This representation is general for both generators and loads since photovoltaic units can be introduced into the model as constant power injections (i.e., $\alpha =0$ and $real\left(s_k\right)\ge 0$). The total power losses are given by (6), a convex quadratic function since the real part of the nodal admittance matrix is positive semidefinite. In fact, the matrix is positive definite if the graph is connected. In that case, the function is not only convex but strictly convex.

$$p_L=real\left(\sum _{k=1}^{3n}\sum _{m=1}^{3n}{y}_{km}{v}_k{v}_m^{*}\right).$$

(6)

The main objective is to reduce the expected value of the active power losses. Therefore, the system is divided into scenarios of generation and load, and the following deterministic representation is obtained:

$$\mathbb{E}\left(p_L\right)=\sum _{i\in \mathbb{W}}{\rho }_i{p}_{L\left(i\right)},$$

(7)

where $\mathbb{W}$ is the set of scenarios, ${\rho }_i$ is the probability associated with the scenario i, and $p_{L\left(i\right)}$ is the corresponding power loss. Observe that this equation is convex.

4. Power Flow Linearization

A linear approximation is proposed for the power flow Equation (5), taking into account the model of the loads. Equation (5) can be rewritten in matrix form as,

$$S_N^{*}=diag\left(V_N^{-\frac{\alpha }{2}}\circ {\left(V_N^{1-\frac{\alpha }{2}}\right)}^{*}\right)\left(Y_S{V}_S+Y_N{V}_N\right).$$

(8)

Now, using the Laurent’s series expansion is linearized Equation (8) around the point $U_N=V_S\otimes {\mathbf{1}}_{n-1}$ where ${\mathbf{1}}_{n-1}$ represents a vector of size $n-1$ whose entries are 1, and ⊗ is the Kronecker product (see [19] for more details of this linearization approach). Therefore, (8) can be approximated by the following affine compact form [25]:

$$S_N^{*}=A+B{V}_N+C{V}_N^{*},$$

(9)

where

$$\begin{array}{cc}\hfill A& =T_N^{*}-{\left(U_N^{1-\alpha /2}\right)}^{*}\circ U_N^{-\alpha /2}\circ \left(Y_N{U}_N\right)-\left(1-\frac{\alpha }{2}\right)\circ {\left(U_N^{-\alpha /2}\right)}^{*}\circ U_N^{-\alpha /2}\circ J_N\circ {\left(U_N\right)}^{*}\hfill \\ & +\frac{\alpha }{2}\circ {\left(U_N^{1-\alpha /2}\right)}^{*}\circ U_N^{-\alpha /2-1}\circ J_N\circ U_N,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill B& =diag\left({\left(U_N^{1-\alpha /2}\right)}^{*}\circ U_N^{-\alpha /2}\right)Y_N-diag\left(\frac{\alpha }{2}\circ {\left(U_N^{1-\alpha /2}\right)}^{*}\circ U_N^{-1-\alpha /2}\circ J_N\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill C& =diag\left(\left(1-\frac{\alpha }{2}\right)\circ {\left(U_N^{-\alpha /2}\right)}^{*}\circ U_N^{-\alpha /2}\circ J_N\right),\hfill \end{array}$$

(12)

and $J_N,T_N$ are auxiliary vectors given by

$$\begin{array}{cc}\hfill J_N& =Y_S{V}_S+Y_N{U}_N,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill T_N& ={\left(U_N^{1-\alpha /2}\circ U_N^{-\alpha /2}\circ J_N\right)}^{*},\hfill \end{array}$$

(14)

in this case, ∘ represents the Hadamard product and ${(·)}^{*}$ is the complex conjugate. The proposed linearization must be implemented for each of the scenarios to be considered.

5. Birkhoff Polytope

A permutation group can represent all the possible configurations of the nodes. This finite group is wholly determined by matrices $M_i$, as given in

Table 1. Feasible permutations for the phase-balancing problem.

Matrix	Value	Permutation	Determinant
$M_1$	$\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$	ABC	+1
$M_2$	$\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right)$	BCA	+1
$M_3$	$\left(\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right)$	CAB	+1
$M_4$	$\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right)$	ACB	−1
$M_5$	$\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right)$	BAC	−1
$M_6$	$\left(\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right)$	CBA	−1

.

Note that this group is discrete; however, it can be represented as the set ${\mathcal{B}}_3\in {\mathbb{R}}^{3\times 3}$ that is schematically depicted in Figure 2. This set constitutes the convex hull of the six permutation matrices given in Table 1 and is called the Birkhoff polytope [26].

Each vertex corresponds to a feasible permutation, and the shaded part corresponds to the convex hull (i.e., the Birkhoff polytope). The vertices are a discrete set, but the entire polytope is a continuous and convex set given by (15).

$$\begin{array}{c}\hfill {\mathcal{B}}_3=\left\{M\in {\mathbb{R}}^{3\times 3}:M{\mathbf{1}}_3={\mathbf{1}}_3,M^{\top }{\mathbf{1}}_3={\mathbf{1}}_3,M\ge 0\right\}.\end{array}$$

(15)

There are many interesting properties of ${\mathcal{B}}_3$ that hold a direct physical interpretation. For example, $∥M∥=1$, meaning that the total power in the original system is equal to the permutation’s total power. In addition, the determinant of the matrix is $det\left(M\right)=\pm 1$. A positive determinant indicates a permutation that does change the sequence of the grid.

The feasible solutions of the problem are in the vertex of the Birkhoff polytope; therefore, the model must guarantee that each entry of $M_i$ is binary, resulting in the following quadratic mixed-integer model:

$$\begin{array}{cc}\hfill min& \sum _{i\in \mathbb{W}}{\rho }_i{P}_{L\left(i\right)}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill P_{L\left(i\right)}& \ge real\left(\sum _{k=1}^{3n}\sum _{m=1}^{3n}{y}_{km}{v}_{ki}{v}_{mi}^{*}\right)\forall i\in \mathbb{W}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill S_{Ni}^{*}& =A+B{V}_{Ni}+C{V}_{Ni}^{*}\forall i\in \mathbb{W}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill S_{ki}& =M_k{\overline{S}}_{ki}\forall k\in N,i\in \mathbb{W}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill M_k{\mathbf{1}}_3& ={\mathbf{1}}_3\forall k\in N\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill M_k^{\top }{\mathbf{1}}_3& ={\mathbf{1}}_3\forall k\in N\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill M_k& \ge 0\forall k\in N\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill M_k& \in {\mathbb{B}}^{3n\times 3n\times s}\hfill \end{array}$$

(23)

Stochastic non-linear mixed-integer problems, such as (1) to (3) are highly complicated to solve, so heuristic algorithms are commonly used to solve them [16]. However, the proposed approximations transform the problem from a general non-linear mixed-integer problem to a convex-quadratic mixed-integer model, given by (16) to (23), since the objective function is quadratic and all the constraints are linear expressions with binary components; this type of problem may be solved in practice, using the Branch and Bound (B&B) algorithm [27]. The B&B algorithm departs from a relaxed problem; in this case, the convex quadratic programming problem QP. After that, a rooted tree is generated by branching non-binary variables, as shown in Figure 3, where m refers to the entries of each matrix M that represents the Birkhoff polytope.

The B&B algorithm works efficiently in this case for two reasons: First, each node corresponds to a QP problem for whose solution there are high-speed algorithms; second, the Birkhoff polytope structure allows a geometric representation that helps to discard nodes in the tree. As a result, a problem that was commonly solved by heuristics is now solvable by using an exact technique, which is the main contribution of this research to the scientific literature.

6. Numerical Validation

This section presents the computational validation of the proposed mixed-integer quadratic problem evaluated in a low-voltage distribution network (i.e., CIGRE network). First, we give its main characteristics, and second, the numerical implementations and results found using the MATLAB software are presented and discussed (Supplementary Materials).

6.1. Test System

Figure 4 illustrates a modified version of the CIGRE low-voltage test system for microgrid applications. This system is a typical residential network composed of 19-nodes, with a peak power demand of 186.9 kW and a nominal voltage of 400 V. This modified CIGRE version was proposed in [28], which includes solar generations at nodes 8, 10, 12, and 16.

6.2. Scenario Generation

The number of scenarios is key to obtain an accurate but tractable model. In this case, three conditions for generation and three load scenarios were considered. These correspond to low, medium, and high generation/demand. These conditions were combined, generating nine generation and demand scenarios with their respective probability. A higher number of scenarios is possible, although the degree of improvement of the solutions to justify the additional computational effort was insufficient. The likelihood of each scenario is given in

Table 2. Scenarios of generation and load.

Scenario	Load/Generation	Probability
1	low/low	0.2210
2	low/medium	0.0443
3	low/high	0.0676
4	medium/low	0.2767
5	medium/medium	0.0554
6	medium/high	0.0845
7	high/low	0.0845
8	high/medium	0.0332
9	high/high	0.0507

.

6.3. Numerical Results

The numerical experiments were performed on a modified version of the CIGRE low-voltage test system [29]. The convex model was solved using CVX, a package for specifying and solving convex programs [30].

Table 3. Results of the configuration nodes.

Node	8	10	11	12	13	14	16	18	19
Permutation	$M_4$	$M_2$	$M_4$	$M_6$	$M_5$	$M_5$	$M_4$	$M_3$	$M_6$

shows the configuration of phase-balancing achieved by the proposed model in generation and demand nodes.

Figure 5 shows the initial power losses for each scenario and their corresponding power loss reduction with the configuration achieved by the MI-QP model.

From Figure 5, it can be observed that the MI-QP model reduced system power losses by, on average, 41.74 %; this indicates that, in any scenario, the system power losses always will reduce concerning the base case. The MI-QP model’s solution guarantees the global optimum for all scenarios with a robust configuration.

To demonstrate that the proposed approach reaches the global optimum of the problem, we evaluate all the possible phase-swapping configurations using an exhaustive search method based on nested loops, which takes about 46 h to evaluate 34 million combinations. However, with the proposed MI-QP, the global optimum finding takes about 104.17 s, demonstrating the efficiency of the proposed approach regarding the required processing times to solve this complex optimization problem.

Figure 6 depicts the box-plot of thee-phase voltages, before and after the phase-balancing. This type of plot, provides a visual representation of the results under different scenarios; it also shows the minimum, maximum, median and first quartile. After the phase-balancing, the voltage profiles are enhanced by reducing the average unbalanced index from 0.3354 to 0.2921. In addition, the voltages moved to acceptable ranges between 0.9 and 1.1 pu. For example, the average of the voltages in Nodes 8, 9 and 10 was below 0.9, in the original system, and changed to values between 0.9 to 1.0 pu in the balanced system. There are, some scenarios with voltages below 0.9, even after balancing; however, this behavior occurs only in few critical scenarios.

It is worth mentioning that the proposed MI-QP model can deal with planning and operating problems associated with phase balancing. However, the proposed model is most suitable for the planning stage, where a working group will configure all the phases before setting up the grid under operation. In the grid operation, all the grid nodes must have a three-phase device composed of switches to change the phase configuration. This option is possible in practice, but it is expensive for low-voltage applications and may introduce unwanted transient behavior.

7. Conclusions

A reformulation of the MINLP model that represents the optimal phase-balancing in low-voltage distribution networks has been proposed in this paper. A mixed-integer quadratic programming (MI-QP) model was developed, whose main advantages compared to the exact MINLP model are as follows: (i) the global optimum can be guaranteed by combining a B&B algorithm with the quadratic approximation of the three-phase power flow model. (ii) classical heuristic and metaheuristic approaches can be avoided to solve this large-scale combinatorial optimization problem; (iii) by preserving the convex structure of the optimal power flow problem, this model allows one to deal with the uncertainties caused by renewable generation and variable loads while working with probabilities; (iv) under the same simulating conditions, the solution reached the same global solution. The latter implies that it does not require statistical tests to prove its efficiency, which is not possible with metaheuristics.

As the numerical results have demonstrated, when the probability scenarios of load consumption and power generation are considered, the proposed MI-QP model’s solution is robust as it reaches the best configuration possible, thus minimizing the total power losses. This is not possible if the grid is optimized only for a particular load/generation condition.

For future works, we propose using local measures in load and generation points and extending our approach to real-time phase-balancing in power distribution networks. We also propose that the proposed MI-QP model be used to locate and size three-phase distributed generators, finding unbalanced distribution networks for power loss minimization, or enlarge the voltage stability margin.