1. Introduction

Energy efficiency and optimization in the operation of electrical energy distribution systems have always been a topic of great importance (^{Claeys et al., 2021}). Furthermore, with the continuous growth of the demand for electrical energy and the need to reduce electrical losses in distribution networks (^{Wilms et al., 2017}), it has become imperative to develop strategies and technologies that allow for improving the operation of these networks by reducing their losses while maintaining their operating limits (^{Marini et al., 2019}; ^{Soltani et al., 2017}).

In this context, the balanced operation of electrical energy distribution systems plays a crucial role in the efficient functioning of electrical systems (^{Bina & Kashefi, 2011}). Therefore, carrying out load balancing at the terminals of a substation is important because various elements and circuits are interconnected. In unbalanced three-phase distribution systems, load imbalance between phases can cause a series of problems, such as overloads on certain equipment, additional losses of electrical energy, and degradation of the quality of the electrical supply (^{Shen et al., 2018}; ^{Garcés-Ruiz, 2022}). Therefore, these challenges motivate us to propose innovative solutions that allow for mitigating load imbalance at substation terminals (^{Jimenez et al., 2022}). In the quest to enhance operational efficiency and minimize power losses in electrical networks, it is essential to investigate novel methodologies and technologies that enable the optimization of load distribution among phases intelligently and effectively (^{Huangfu et al., 2024}).

In the specialized literature, various methodologies have been proposed to address the issue of phase balance or load redistribution. The work by ^{Al-Kharsan et al (2020)}, applied a particle swarm optimization (PSO) to optimize the balance of phases in unbalanced three-phase distribution systems. The PSO algorithm demonstrated its capability to explore optimal phase configurations to minimize the imbalance of the unbalanced three-phase system. The above was achieved by adjusting the inertia weight to identify the most efficient configuration and, ultimately, enhance the performance of the PSO algorithm.

The study by ^{Cortés-Caicedo et al. (2021)}, implemented a methodology for optimal balancing in unbalanced three-phase distribution systems using a discretized vortex search algorithm (DVSA). The DVSA aimed to minimize power losses in a three-phase system by determining the optimal configuration of phase connections for the system's demand nodes. This algorithm was evaluated in IEEE 8-, 25-, and 37-node test systems and compared with the Chu-Beasley genetic algorithm. The numerical results showed the efficiency of the DVSA in reducing power losses in unbalanced three-phase distribution systems.

The work by ^{Montoya et al. (2021)}, proposed an optimal balancing of phases in three-phase electrical distribution systems using a mixed-integer quadratic convex formulation, which minimized the sum of the squared currents through the distribution lines. This work was tested in three radial test systems and compared with three metaheuristic techniques. The mixed-integer quadratic convex formulation exhibited superior performance compared to the other techniques.

The authors of ^{Bohórquez-Álvarez et al. (2023)}, addressed the challenge of minimizing power losses in unbalanced distribution networks using a convex approximation. The proposed optimization approach aimed to minimize total power losses in a three-phase network by utilizing the concept of electric momentum. The above was achieved by simplifying the power balance restrictions and transforming the objective function into a strictly convex form, resulting in efficient loss minimization. The GAMS software, with its solvers CPLEX, SBB, and XPRESS, was used to compare the performance of the proposed approximation.

Other studies have also explored the optimal phase balance in three-phase distribution networks with asymmetric loads using methods such as specialized genetic algorithm (^{Granada-Echeverri et al., 2012}), smart meter data-based algorithm (^{Grigoraș et al., 2020}), modified PSO algorithm (^{Tuppadung & Kurutach, 2006}), a heuristic algorithm that uses pole measurement (^{El Hassan et al., 2022}), Derivative-Free algorithm (^{Montoya et al., 2021}), and genetic algorithms incorporating group theory (^{Garcés et al., 2020}), among several other proposed approaches.

In light of the information provided, this research compares three mixed-integer quadratic programming models designed to address the expected imbalance in active, reactive, or apparent power at substation terminals within three-phase asymmetric networks. These optimization models, which belong to the domain of mixed-integer convex models, guarantee solution repeatability in each execution by utilizing the Gurobi solver in the Julia programming language, supplemented by the HiGHS solver in the JuMP optimization framework. These contributions offer an innovative approach to addressing power imbalances in asymmetric networks, providing a robust and reliable tool for the efficient planning and operation of electrical systems.

This paper is organized as follows: Section 2 describes the proposed methodology for optimal load balancing in terminals of the substation. Section 3 presents two test systems to assess the proposed methodology, as well as their main results. Section 4 presents the main concluding remarks of the study.

2. Methodology

In the proposed methodology for the optimal load balancing in terminals of the substation and the power loss minimization, two main elements are considered. The first element corresponds to the mixed-integer quadratic models and the second element is the three-phase power flow formulation via Broyden’s formulation to evaluate the final load redistribution. Each one of these elements is widely described below.

A. Mixed-integer quadratic models

The problem of the effective load balancing in terminals of the substations assumes that in a three-phase distribution network, all the loads can be summarized as an equivalent load per phase, as presented in Figure 1.

Figure 1 Equivalent load connections in terminals of the substation. 

In Figure 1 it is observed that the voltage profiles in the substation are measured with respect to the neutral point n, which are defined as Van, Vbn, and Vcn, respectively, in addition, San, Sbn and Scn represent the equivalent apparent power load connections in terminals of the substation (i.e., complex values that can be separated in their real and imaginary parts called active and reactive powers). Note that for a three-phase, load, there exists six possible connections in terminals of a particular k node. Table 1 summarizes these possibilities, where Mkj represents the possible load rotation at node k, which have six options, i.e., j=1,2,…,6, being Sk,0abc the initial load connection and Sk,fabc the final load connection at node k.

Table 1 Possible load redistribution in a particular k node 

From Table 1 it is observed for any k node it is possible to obtain a load rotation considering the possible load rotation Mj related with the binary variable xkj, using a vector representation, as defined in Equation

Sk,fabc=∑j=16xkjMkjSk,0abc, ∀k∈N (1)

where N is the set containing all the nodes of the network, in addition, from Equation (1) the real and imaginary parts of the apparent power (i.e., Pk,fabc and Qk,fabc) can be split into Equations (2) and (3).

Pk,fabc=∑j=16xkjMkjPk,0abc, ∀k∈N (2)

Qk,fabc=∑j=16xkjMkjQk,0abc, ∀k∈N (3)

To ensure that at each k node only and only one possible load distribution (as shown in Table 1), the following constraint is imposed to the binary variables, as defined in Equation (2).

∑j=16xkj=1, ∀k∈N (4)

Regarding the objective function, the problem of the optimal load redistribution in terminals of the substation can be done with respect to the active, reactive power, apparent powers. The proposed three quadratic objective functions are defined from (5) to (7).

min⁡F1=100%3Pav2∑p∈{a,b,c}Pk,fp-Pav2 (5)

min⁡F2=100%3Qav2∑p∈{a,b,c}Qk,fp-Qav2 (6)

min⁡F3=12F1+F2 (7)

where F1 is the average active power unbalance, F2 is the average reactive power unbalance and F3 represents the average apparent power unbalance, being Pav and Qav the ideal active and reactive power consumption per phase.

3. Results and discussion

In this section is presented the test feeders under analysis, which consist to the 15-bus and 35-bus unbalanced distribution grids, and the numerical validations of the proposed optimization models for optimal load balancing the total power consumptions in terminals of the substation and their effect of the expected power losses level.

3.1 Test feeders

The medium-voltage distribution networks under analysis correspond to the 15- and 35-bus grids, which are unbalanced distribution systems with radial topology.

The 15-bus grid

The 15-bus grid is a distribution network operating with a line-to-line voltage of 13.2 kV composed of 14 distribution lines with 6 type of branches impedances. The total constant power load consumption in phases a, b, and c is 9605+j5226 kVA, 6480+j4940 kVA, and 11977+j8778 kVA, respectively. Table 2 and 3 defines the power consumption per node and the nodal connection (branches information including the impedance type).

Table 2 Load information for each phase in the 15-bus grid 

Table 3 Information of nodal connection and impedance type for the 15-bus grid. 

The set of impedances associated with each type of conductor listed in Table 3 is presented in Table 4.

Table 4 Impedance assignable to each type of conductor listed in Table 3. 

The 35-bus system

The 35-bus grid corresponds to a radial distribution network composed of 35 nodes and 34 distribution lines. The line-to-ground voltage in terminals of the substation corresponds to 15 kV. The total constant power load consumption in phases a, b, and c 4435+j1940 kVA, 2868+j2059 kVA, and 2925+j1370 kVA, respectively. Table 5 and 6 defines the power consumption per node and the nodal connection (branches information including the impedance type).

Table 5 Load information for each phase in the 35-bus system. 

Table 6 Impedance information for each branch in the 35-bus system 

It is worth mentioning that in the case of the impedance information of the 35-bus system, the impedance assignable to each branch is diagonal, i.e., this presents the structure of defined in Equation (15).

Zkm=Rkm+jXkm000Rkm+jXkm000Rkm+jXkm (15)

3.2 Results

In this subsection, the numerical analysis of both test feeders is reported. The benchmark case and the solution of each of the three-proposed optimization models are also reported, added with the final value of the power losses. Note that all the numerical results reported in this section were reached with the usage of the Gurobi solver in the Julia programming language by its combinations with the HiGHS solver in the JuMP optimization environment (^{Bezanson et al., 2017}; ^{Lubin et al., 2022}).

Numerical results for the 15-bus grid

Table 7 reports the solution reached with the solution of each objective function defined in the optimization model (1)-(7) by comparing the initial and final values of each objective function and their corresponding power losses for the 15-bus system simulation case.

Table 7 Numerical results in the 15-bus grid. 

Numerical results in Table 7 shows that:

Each optimization model using the quadratic formulation allows minimizing the total expected power unbalanced aggregated in terminals of the substation since function F1, F2, and F3 is minimized as much as possible, owing the theoretical global optimal solution of each one is zero.

The expected level of active power losses, as defined by formula (14), is reduced with respect the benchmark case. When F1 is minimized, the reduction is about 23.1644 kW, i.e., 17.2550 %. In the case of minimizing F2, the total reduction in power losses is about 19.7167 kW, i.e., 14.6869 %. Finally, in the case of minimizing function F3, the total reduction in power losses is about 12.9624 %, i.e., 17.4016 kW.

Note that in the case of the 15-bus grid, the proposed optimization model (1)-(7) presents a better result when the objective function regarding the active power balance is minimized; however, it is recommended that for each distribution network, all the performance indexes in Equations (5)-(7) must be tested, since the numerical performance regarding power losses can be vary, which make the optimization result for a particular objective function non-generalizable.

Numerical results for the 35-bus grid

Table 8 reports the solution reached with the solution of each objective function defined in the optimization model (1)-(7) by comparing the initial and final values of each objective function and their corresponding power losses for the 35-bus system simulation case.

Table 8 Numerical results in the 35-bus grid 

Numerical results in Table 8 shows that:

The values reached for each objective function are near to zero, i.e., with values lower than 1x10^-05 %, which implies that the proposed quadratic convex models reach the optimal solution for each one of the objective functions analyzed. Note that for many practical applications the expected solution equal to zero is not reachable since it is highly probably that with the discrete movements in loads the ideal power consumption per phase is not practically obtained.

The expected level of active power losses, as defined by formula (14), is reduced with respect the benchmark case. When F1 is minimized, the reduction is about 24.0206 kW, i.e., 5.0771 %. In the case of minimizing F2, the total reduction in power losses is about 36.6142 kW, i.e., 7.7389 %. Finally, in the case of minimizing function F3, the total reduction in power losses is about 5.6401%, i.e., 26.6843 kW.

Results in the case of the 35-bus system shows that the proposed optimization model (1)-(7), allows reaching the lowest possible power losses when function F2 is minimized. These results confirm that it is not possible to correlate the minimum power losses results with a particular objective function, which implies that for each distribution network under analysis, all the quadratic functions must be evaluated to obtain the higher reduction in power losses via Broyden’s power flow approach.

4. Conclusions

This article has presented three mixed-integer quadratic programming models to reduce the expected active, reactive, or apparent power unbalance in terminals of the substations in three-phase asymmetric networks. The proposed optimization models belong to the family of mixed-integer convex models and allows ensuring the solution repeatability at each execution using the Gurobi solver in the Julia programming language with the help of the HiGHS solver in the JuMP optimization environment.

A two-stage optimization approach was developed to evaluate the initial and final expected power losses in the three-phase network via the power flow approach based on Broyden’s numerical method. The numerical results in the 15- and 3 bus grids demonstrated that:

In the case of the 15-bus grid, the best optimization model corresponds to the minimization of objective function F1 with an expected reduction in the power losses of about 17.2550 %. However, in the case of the 35-bus grid, the best objective function was F2, since it allowed reduction in power losses of about 7.7389 %.

The aforementioned results showed that is not possible to relate the maximum power losses reduction calculated with a power flow approach (i.e., Broyden’s numerical method) with a particular objective function (see Equations (5)-(7)), since for each distribution network is mandatory to evaluate the proposed solution methodology in order to identify which is the better objective function indicator when the final objective is minimizing the expected grid power losses.

To improve the numerical performance of the proposed solution methodology it is recommended to include in the mixed-integer quadratic programming formulation the effect of the distribution branches. This effect could be modeled using a linear approximation of the power flow equations, which can help to have better numerical values for the expected power losses when compared to the proposed two-stage optimization approach.

Possible future works derived from this research can be the following: (a) to extend the proposed two-stage optimization approach to bipolar direct current networks with asymmetric loading; (b) the addition in the mixed-integer quadratic formulation of the current flow effect in branches using a linear approximation for the three-phase power flow problem, and (c) the application of hybrid optimization algorithms based on the combination of metaheuristic optimizers with Broyden numerical method. The metaheuristic will be entrusted with determining the load connection per node, while the power flow method (Broyden’s numerical method) will be encharged with determining the power losses level.