A Mixed-Integer Convex Model for the Optimal Placement and Sizing of Distributed Generators in Power Distribution Networks

: The optimal placement and sizing of distributed generators is a classical problem in power distribution networks that is usually solved using heuristic algorithms due to its high complexity. This paper proposes a different approach based on a mixed-integer second-order cone programming (MI-SOCP) model that ensures the global optimum of the relaxed optimization model. Second-order cone programming (SOCP) has demonstrated to be an efﬁcient alternative to cope with the non-convexity of the power ﬂow equations in power distribution networks. Of relatively new interest to the power systems community is the extension to MI-SOCP models. The proposed model is an approximation. However, numerical validations in the IEEE 33-bus and IEEE 69-bus test systems for unity and variable power factor conﬁrm that the proposed MI-SOCP ﬁnds the best solutions reported in the literature. Being an exact technique, the proposed model allows minimum processing times and zero standard deviation, i.e., the same optimum is guaranteed at each time that the MI-SOCP model is solved (a signiﬁcant advantage in comparison to metaheuristics). Additionally, load and photovoltaic generation curves for the IEEE 69-node test system are included to demonstrate the applicability of the proposed MI-SOCP to solve the problem of the optimal location and sizing of renewable generators using the multi-period optimal power ﬂow formulation. Therefore, the proposed MI-SOCP also guarantees the global optimum ﬁnding, in contrast to local solutions achieved with mixed-integer nonlinear programming solvers available in the GAMS optimization software. All the simulations were carried out via MATLAB software with the CVX package and Gurobi solver.


Introduction
Distributed energy resources (DERs) are an attractive technology to install in power distribution networks due to their technical, economic, and environmental benefits [1]. In general, a DER is a small power plant that can employ different kinds of energy source such as small-hydro, micro-turbine, fuel cell, methane-generating bioreactors, gas-turbine, and renewable distributed generation (wind and solar power). Typically, DERs are located close to the customers in a distribution network, which has allowed decreasing power losses associated with transmission lines and improving the voltage profiles [2,3]. However, DER can generate adverse effects in electrical distribution networks, such as dis-coordination of protection schemes, increased current fault level, bidirectional power flows, transmission line overload, and voltage profile deterioration [4] if these are not correctly placed and sized. An optimization model for the optimal location and sizing of DERs is required to avoid these possible problems.
This problem has been extensively studied in the literature with several mathematical formulations and methodologies. The problem is complicated to solve since its formulation combines discrete and continuous variables with nonlinear/non-convex constraints. Therefore, the global optimum cannot be guaranteed [5]. The problem has been typically solved using master-slave strategies [6]. The master strategy chooses the DERs placement and size, while the slave strategy solves an optimal power flow problem. For the master strategy, metaheuristic optimization techniques have been commonly implemented, such as Chu-Beasley genetic algorithms [7], tabu search algorithms [8], simulated annealing methods [9], ant-lion optimizers [10], whale optimization algorithm [11], krill herd algorithms [12], teaching-based learning optimizers [13], population-based incremental learning [6], imperialist competitive algorithms [14], harmonic search algorithms [15], and bat and firefly algorithms [16], among others. Although these optimization techniques ensure good solutions, none of them can guarantee the optimal solution, and some are not well justified mathematically [17]. These techniques include some degree of randomness, which leads to slightly different solutions from one execution to another. In addition, they contain many adjustment parameters that make them highly dependent on the programmer, adding to the fact that statistical analysis must be required to determine their reproducible capacities [18][19][20].
Recent approaches have combined heuristic methods with second-order cone programming (SOCP) models to solve the problem of the optimal location and sizing of DERs in distribution networks as the case of the discrete version of the sine-cosine algorithm reported in [21] and the genetic algorithm proposed in [22] however, these approaches only ensures the optimal solution for the continuous part of the problem, i.e., the DERs' sizing problem, while for the integer part, i.e., the DERs location problem, it is not possible due to the heuristic nature of the discrete algorithms. The heuristic nature of the the integer part does not allow for ensuring the global optimum reaching, which implies that statistical test are required to verify their efficiency [19].
Unlike heuristics, this paper proposes an exact formulation based on recent advances in convex optimization for power systems applications [23]. The optimal placement and sizing of DERs are transformed into a mixed-integer second-order cone programming (MI-SOCP) model. The problem is solved via the Branch & Bound method, which guarantees optimal global optimization if each sub-problem of optimization is convex. For this purpose, the optimal power flow becomes a convex problem through a relation. Thus, the main contributions of this work are summarized below: • A reformulation of the classical model represents the problem of the optimal location and sizing of DERs in power distribution networks by using an MI-SOCP model. This model guarantees the global optimum of the relaxed problem, being indeed a suitable solution to the problem for practical purposes; • The application of the proposed MI-SOCP model to the optimal sizing and location of DERs considering the unity power factor and variable one is analyzed. Its numerical results confirm that for the unity power factor found in the best solutions reported in scientific literature and the case variable power factor, better solutions are found when compared to the literature reports. In both simulation cases, with minimal computation efforts and reproducible results, i.e., with the assurance of finding the same optimum solution at each simulation run.
The remainder of the paper is organized as follows: Section 2 describes the mathematical model for the optimal placement and sizing of DERs in power distribution networks. Section 3 presents the optimal power flow convex model based on the SOCP approximation. Section 4 describes the MI-SOCP convex model based on the Branch & Bound method and SOCP approximation. In Section 5, the configurations of the IEEE 33-bus and IEEE 69-bus system and simulation cases are presented. Section 6 shows the MI-SOCP convex model's numerical results compared to multiples strategies previously reported. Finally, some conclusions are drawn in Section 7.

General Mathematical Model
The mathematical model for the optimal placement and sizing of DERs in electrical distribution networks is a mixed-integer nonlinear problem. This implies that it is composed of integers and continuous variables. The grid is represented as an oriented graph G = {N , E }, where N = {1, 2, . . . , n} represents the set of nodes of the power distribution network and E ⊂ N × N the set of branches. Without loss of generality, the slack node is assumed to be k = 1. For a compact representation, all variables are presented in the set complex, considering that they can be separated into real and imaginary parts. The location or not of the DERs in the network shape the integer (binary) variables, while the continuous variables are defined in the optimal power flow formulation. The optimal problem tackled in this paper is described in Nomenclature.

Objective Function
The main objective of the optimal placement and sizing of DERs in electrical distribution network is to minimize the active power loss function, which is given by (1), where p L represents the power loss of the network, y km is entry km of the line admittance matrix Y bus , and (·) * denotes the complex conjugate; and v k and v k are the voltages phasor at nodes k and m, respectively. This objective function constitutes a convex quadratic form as will be demonstrated in Section 2.3.

Constraints
The main constraint of the model is given by the set of power flow equations that are given by (2): where s k , is the apparent power generated at node k by a existing generators, s new k is the apparent power generated by new DERs at node k, and d k is the apparent power demanded at node k. Notice this constraint constitutes a set of 2n non-affine equality equations (n for the real part and n for the imaginary part) which make the problem non-convex.
The voltage in the slack node (substation) is supposed to be v 1 = 1 + j0; the rest of voltages are limited as follows: where V nom is rated voltage of the system and δ is the maximum deviation given by the grid code (usually δ is a constant between 0.05pu to 0.1pu). The left-hand side of this constraint is non-convex, however, the right hand side is. Flows in the distribution lines are limited by the following box constrain which is, of course, convex: where s km is the power flow through line km and s km is the maximum apparent power flow allowed. Capacity of generators, both existing and new are considered as follows: where s k is the maximum apparent power generated at node k by a existing generators, s new is the maximum apparent power generated by new DERs at node k, and z k ∈ {0, 1}, ∀k ∈ N (7) which denotes the binary variable of the problem, which has a value of 1 if a DER is installed at node k or 0 otherwise. There is a limit of DERs to be installed in the system, given by (8), For the sake of completeness, the entire model is presented below: Model (9) is a non-convex non-linear integer programming problem. These problems are highly complex, hence the need to use approximate techniques for their solution. However, some constraints are convex and may be used to generate a tractable problem.

Identification of Model Properties
The model for the optimal placement and sizing of DERs in power distribution networks is highly complex, hence the use of heuristic algorithms. However, some parts of the model exhibit suitable properties to be included in a convex optimization model. Table 1 shows this properties. Table 1. Identification of model properties.

Equation Property
(1) convex objective (quadratic-convex form) convex (box constraint) (6) convex (box constraint) (8) convex (affine inequality constraint) (7) binary constraint Notice the objective function can be represented in a real domain by splitting y km = g km + jb km and v = v real + jv imag , resulting in the following equivalent expression: ).
Since G = [g km ] ∈ R n×n is positive semidefinite (note that G can be calculated as G = AG p A , where A is the incidence matrix and G p is a diagonal matrix with the resisting effect of each branch), then p L is a convex function. Most of the constraints are also convex. Therefore, the main source of complexity is the power flow equations and the binary nature of variables z. A suitable approximation is required in order to obtain a tractable model, as presented in the following section.

SOCP Optimization
The second order cone programming (SOCP) is sub-field of convex optimization that deals with an optimization model with an affine linear and second order cones constraints [24]. Tailored algorithms exist for these type of problems that efficiently solve large SOCP models in only microseconds. A second order cone is a convex set defined by the following constraint: x ≤ z where x ∈ R n and z ∈ R. x is the norm-2 of the vector x (hence it is also called a Euclidean cone). Figure 1 shows a second order cone in R 3 , clearly a convex set. For further details about optimization convex, see [25]. x 1 x 2 z Figure 1. Representation of the second order cone Ω = { x ≤ z}. with x ∈ R 2 and z ∈ R.

SOCP Approximation for the Power Flow Equations
The main idea behind the SOCP of the power flow is to transform non-convex (bilinear) constraints into second order cone constraints. The formulation allows transforming the power flow equations, that are a nonlinear/non-convex problem, into a series of convex constraints. The proposed SOCP model is built on the approximation proposed by Low in [26]. A new matrix W = [w km ] ∈ C n×n is defined, as follows: Therefore, power flow equations (2) can be represented as the following equivalent equation: Note (11) is affine and the non-convexity appears in (10), which in turns implies the following expression: A new vector U ∈ R n is defined, with entries u k = v k 2 , then (12) is transformed into (13): Being a complex number, w km stores information related to the angle of the voltages whereas u k , u m stores information related to the magnitude. We use the following hyperbolic representation for the right hand side of (13): Resulting in the following equivalent expression: At this point, the problem is still not-convex, however, (14) relaxed and transformed into the following constraint: Notice that (15) represents a second order cone, namely: where the term inside the norm in the left hand side of the inequality is a vector in C 2 (square (16) to demonstrate it is equivalent to (15)). The problem is now convex and tractable, except for the binary constraints given by (7). The problem must be solved in the framework of MI-SOCP.
The new decision variables u k allow also to transform (3) into the following convex constraint: ( This constraint is equivalent to (3) since 1 − δ ≥ 0 and therefore the inequality holds even if it is squared (i.e., (17) is not an approximation but an equivalent constraint).
The objective function (1) although already convex, is transformed into a simpler affine equation using the new variables w km : min p L = real(y km w km ) remember that both y km and w km are represented as complex variables but the objective p L is a real function. In theory, these variables are separated in real and imaginary parts. However, the convex representation is more compact and can be directly implemented in practice using languages such as Matlab and/or Python.
Combining all the aforementioned considerations, the MI-SOCP model is obtained: Notice the model is now a SOCP except for the last constraint that is binary, leading to a MI-SOCP model.

Return to the Original Variables
The model for the SOC approximation is given as a function of w km and u k instead of the voltages v k . However, it is possible to approximate the original voltages by the following two-step procedure: First, the magnitude of the voltages is computed as v k = √ u k . This value exists and is real since u k ≥ 0 and second, the angle of the voltages is calculated from θ km = ang(w km ) in a forward iteration, starting from θ 1 = 0. Therefore, a power flow calculation is not required after the optimization problem is solved.

Solution of the MI-SOCP Model
Mixed-integer second-order cone programming (MI-SOCP) are problems that includes the following type of constraints: where decision variables include continuous x and binary variables z; A k are real matrices, b k , α k , β k are real vectors, and γ k are constants for each constraint k.
Like most of the integer programming problems, an MI-SOCP can be solved using a modified version of Branch & Bound (B&B), as shown in Figure 2. At each iteration, the B&B solves a SOCP problem that uses an interior point method, specially designed for this type of problem [24]. The method is benefited for the properties of the SOCP problems concerning convexity and fast convergence of the interior point methods [27]. There are different versions of the interior point method for convex optimization problems. In general, these methods can be classified into two groups: Primal and primaldual methods. Primal methods are based on the seminal paper of Karmarkar [28], which can be straightforwardly extended to SOCP problems. Primal-dual requires more effort to be applied. However, numerical experiments have demonstrated advantages in terms of convergence [29].
The B&B method departs from a solution of a SOCP problem where the variables z k ∈ [0, 1] are considered continuous. This solution is named as N0 and constitutes a lower bound for the MI-SOCP problem. Next, a series of SOCP optimization problems are solved for integer partitions of z k , as depicted in Figure 2. Each SOCP problem has continuous variables, and hence it can be solved efficiently despite having a high amount of variables. The method continues until a binary solution is obtained. Other approaches, such as outer approximation, cutting-plane methods, and generalized benders decomposition, can be used to solve MI-SOCP (see [27] and the references therein for a complete review of MI-SOCP). However, most of the solvers that allow solving MI-SOCP use B&B. Since our approach is application-oriented, we will rely on the solution given by these solvers.

Test Systems and Simulation Cases
The MI-SOCP formulation is validated in two widely-used test systems for the problem of optimal location and sizing of DERs in power distribution networks. The test systems are the IEEE 33-and 69-bus systems. The most relevant information about these test systems is presented below.

IEEE 33-Bus System
This test system operates at 12.66 kV and consists of 33 nodes, 32 branches, and a unique slack node located at node 1. The configuration of this feeder is depicted in Figure 3a, and its parameters can be found in [3]. Its initial active power losses are equal to 210.9876 kW. The total active and reactive power demands of this system are 3715 kW and 2300 kVAr, respectively. The voltage and power base values considered are 12.66 kV and 1 MW, respectively.

IEEE 69-Bus System
This test system works at 12.66 kV, which is composed of 69 nodes, 68 branches, and a unique slack node located at node 1. The configuration of this feeder is shown in Figure 3b, and its parameters can be consulted in [3]. Its initial active power losses are equal to 224.9520 kW. The total active and reactive power demands of this system are 3890.7 kW and 2693.6 kVAr, respectively. The voltage and power base values chosen are 12.66 kV and 1 MW, respectively.

Simulation Cases
To validate the effectiveness and robustness of the proposed MI-SOCP model, we consider two simulation cases as follows: • Case-1: This case considers the possibility of installing three DERs. For the IEEE 33-bus system, the capacity of DERs is from 0 kW to 1.2 MW for each one. While for the IEEE 69-bus system, the capacity of DERs is from 0 kW to 2 MW for each one. In addition, each DER is assumed to operate with a unit power factor; • Case-2: This case also considers the possibility of installing three DERs with the capacity to generate reactive power, i.e., variable power factor.

Results
The MI-SOCP convex model has been solved using CVX, and Gurobi [30]. It is important to mention that the solution MI-SOCP model shown in the results is a solution recovery in original variables using the exact power flow model presented in Section 2.

Remark 1.
The comparative shown for the distributed generators' optimal locations and sizes have been taken from the original reports. This research's main interest is to present an effective and robust optimization methodology via MI-SOCP to place and size DGs guaranteeing the global optimum finding when this is compared with the most classical and well-known metaheuristic approaches. Table 2 presents the results for the optimal placement and sizing of DERs of the proposed convex model and those reported in the literature.
From Table 2, we can conclude that: • The solution of the proposed MI-SOCP model can ensure that it is the optimal global solution for the IEEE 33-bus and IEEE 69-bus systems. This is because it has evaluated all the possible nodal combinations by using an exhaustive search combined with a classical optimal power flow (OPF) model solved via interior point methods.  Figure 4 depicts the power loss reductions reached for each methodology presented in Table 2. Note that the proposed convex model achieves the greatest reduction in power losses for the IEEE 33-bus system (see Figure 4a). While the IEEE 69-bus system, the MI-SOCP, and MSSA methodologies present the same power loss reduction (see Table 2), even their DERs sizing are different. This indicates that there are multiple global optima for this system, and there may be other sizes of DERs that reach the same power losses. The processing time for the IEEE 33-and 69-bus systems is 6.2 s and 49.7 s, respectively.

Results of Case-2
In Case-2, the proposed convex model is compared to improved analytical (IA) [40], particle swarm optimization (PSO) [41], and MINLP. Table 3 shows the results for IEEE 33-buses and IEEE 69-buses systems for the optimal placement and sizing of DERs.   Table 3 that the MI-SOCP convex model proposed has a superior performance for placement and sizing of multiple DERs for both test systems with real and reactive power capability when compared with the other three methodologies. In addition, it can also be observed that the MI-SOCP and MINLP present the same placement for the DERs in the IEEE 33-bus system, but the proposed convex model achieves a better solution than the MINLP model. This is because the MINLP model can not guarantee the global optimum. While for the IEEE 69-bus system, none of the other methodologies could find the combination of buses for the location of the DERs. The power loss reduction in Table 3 confirms that the MI-SOCP convex model gives a maximum reduction of power losses. However, the total DERs power for the proposed convex model in the IEEE 33-bus system is slightly higher than the MINLP model. While for the IEEE 33-bus system, the total DERs power is slightly lower in the MI-SOCP convex model. Finally, the IEEE 33-and 69-bus systems' processing time is 11.8 s and 75.1 s, respectively.

Comparison with a Hybrid Optimization Approach
To demonstrate the effectiveness and robustness of the proposed MI-SOCP model to determine the optimal location and sizing of renewable DGs in AC distribution networks, here, we present a comparison with hybrid optimization approaches that combines a discrete metaheuristic with a SOCP formulation to locate and size DGs in electrical distribution grids. This approach was initially presented in [21], and it works with the discrete version of the sine-cosine algorithm to determine the nodes where DGs will be located, while the SOCP is used to determine their optimal sizes (this approach is named DSCA-SOCP). In addition, we include comparisons with MINLP solvers available in the GAMS optimization tool. The renewable energy resources correspond to photovoltaic (PV) sources considering a typically daily generation curve and a daily load curve. These curves are presented in Figure 5. This comparison is considered the IEEE 69 node test feeder since these are more complex than the IEEE 33 node test feeder. Furthermore, the test feeder has about 50,116 possible combinations to locate three PV sources. In the case of MINLP solvers, we consider the following optimization tools available in GAMS: Solvers BONMIN, COUENNE, and DICOPT, respectively. The complete comparison of the proposed MI-SOCP model with the DSCA-SOCP and the GAMS solvers is reported in Table 4.
The results in Table 4 demonstrate that: (i) With the location of three photovoltaic sources, the proposed approach, i.e., the proposed MI-SOCP and the DSCA-SOCP approaches, reduces the daily energy loss per day to about 638.1187 kWh/day, i.e., 23.93%; while the best GAMS approach using the COUENNE solver reach a reduction about 23.84%.The solutions found by the different MINLP solvers in GAMS demonstrate that these stuck in optimal local solutions when compared with the optimal solution found by the MI-SOCP and the DSCA-SOCP methodologies, and (ii) in all the solutions presented in Table 4 nodes higher than 60 show the high power injection regarding PV penetration, which allows observing that these nodes are more sensitive to power injections when compared with the remainder of buses. It is worth motioning that even if the DSCA-SOCP reach the same optimal solution of the proposed MI-SOCP approach, these methodologies are not comparable due to: The DSCA-SOCP approach being a metaheuristic approach that requires a statistical test based on multiple evaluations to determine its capability of finding the global optimum. In this sense, after 100 consecutive evaluations, the DSCA-SOCP methodol-ogy reach the solution reported in Table 4 only 5 times, i.e., an efficiency of 5%, while the proposed MI-SOCP finds the same solution at each running, due to the convex structure of the solution space that guarantees the global optimum finding without recurring to statistical tests; Each evaluation of the DSCA-SOCP approach (i.e., one complete evaluation) takes about 16,747.35 s with the main problem being that it is not possible to ensure the finding of the optimal solution while the proposed MI-SOCP only takes 315.60 s to determine the global optimal solution. These processing times show that our proposal is about 53 times more efficient than the DSCA-SOCP approach.

Conclusions
The problem in the optimal location and sizing of DERs in electrical AC distribution networks was addressed in this paper, from the exact mathematical optimization point of view. To solve the location problem, the classical and well-known Branch & Bound method was employed and in the case of the sizing problem, it was proposed as second-order cone relaxation of the optimal power flow problem. The combination of both approaches generated the MI-SOCP method proposed. The main advantages of our optimization proposal regarding current literature were: (i) The guarantee of finding the global optimum for each binary variable combination explored in the Branch & Bound stage since each node in this problem has a convex structure. (ii) The non-dependency of parametric adjustments since the MI-SOCP model is an exact approach that only requires an adequate representation of the optimal power flow problem, which allows repeatability in the solution at each simulation run. (iii) The easy extension of the proposed model from the unity power factor case to the variable one allows sizing the apparent power required in the DERs for minimizing the total grid power losses. Numerical results in the IEEE 33and 69-bus systems demonstrated that for the unity power factor case were found too be the best solutions reported in the literature. While in the case of the variable power factor being improved, the best solutions in scientific research, in both cases with minimum computation efforts, i.e., low-processing times.
Numerical validations considering daily load and generation curves for photovoltaic sources demonstrated that the proposed MI-SOCP has the ability to achieve the global optimum, while hybrid metaheuristic approaches, i.e., the DSCA-SOCP, requires a statistical test to find the optimal solution without guaranteeing that these can be found. In addition, comparisons with MINLP solvers available for the GAMS optimization package demonstrates that the non-convex structure of the problem create an issue of having them stuck in local optimal solutions.
In future works, we suggest applying the proposed MI-SOCP problem to the following issues: (i) Optimal location of fixed-step capacitor banks with continuous and discrete capacities. (ii) The extension of the convex power flow to the optimal operation of battery packages in AC distribution networks, as well as in exploring their optimal locations and sizes, and (iii) the reformulation of the phase-balancing problem for AC networks from the convex point of view to minimize grid power losses with the guarantee of global convergence. Funding: This research was funded by the Agencia Estatal de Investigación, Spain (AEI) and the Fondo Europeo de Desarrollo Regional (FEDER) aimed at the Challenges of Society (grant No. ENE 2017-83860-R "Nuevos servicios de red para microredes renovables inteligentes. Contribución a la generación distribuida residencial").

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.
Acknowledgments: This work was supported in part by the Centro de Investigación y Desarrollo Científico de la Universidad Distrital Francisco José de Caldas under grant 1643-12-2020 associated with the project: "Desarrollo de una metodología de optimización para la gestión óptima de recursos energéticos distribuidos en redes de distribución de energía eléctrica." and in part by the Dirección de Investigaciones de la Universidad Tecnológica de Bolívar under grant PS2020002 associated with the project: "Ubicación óptima de bancos de capacitores de paso fijo en redes eléctricas de distribución para reducción de costos y pérdidas de energía: Aplicación de métodos exactos y metaheurísticos."

Conflicts of Interest:
The authors declare no conflicts of interest.
Nomenclature δ maximum deviation of the voltages (V). B set of binary numbers. C set of complex numbers. E set of branches. N set of nodes.
x, x limits minimum and maximum of the corresponding variable x (VA). d k load at node k (VA). n T number of total DERs to be placed in the grid. p L total power loss (W). s k apparent power of existing generators at node k (VA). s new k apparent power of new distributed resources at node k (VA). s km load flow in the line km (VA). u k , u m magnitude of voltages squared (V 2 ). V nom rated voltage of network (V). v k , v n nodal voltages (V). w km entries of auxiliary matrix (V 2 ). y km entries of the Y bus (Ω). z k binary variable that indicates placement of DERs in node k.