# Integrated Planning of MV/LV Distribution Systems with DG Using Single Solution-Based Metaheuristics with a Novel Neighborhood Search Method Based on the Zbus Matrix

Authors: Alejandro Valencia; Ricardo A. Hincapie; Ramon A. Gallego
Journal: Journal of Electrical and Computer Engineering
 (2022)

---

## Abstract

This paper presents a new methodology for the optimal integrated planning of medium- and low-voltage distribution systems, considering the location and sizing of distributed generation. The integrated problem is formulated as a mixed-integer nonlinear model and, to solve it, two well-known optimization algorithms (simulated annealing and iterated local search) are used. The intensification and diversification processes are usually the bottleneck of metaheuristic techniques for solving complex problems. To overcome such complexity, a new neighborhood search method based on the Zbus matrix (NSZM) is proposed, to explore the solution space more efficiently and effectively for both algorithms. The proposed methodology is validated and tested on a real distribution system taken from the literature. The results obtained are better than those reported in the literature. To verify the efficiency of the new NSZM method, the Wilcoxon signed ranks test is used to measure the performance behavior of the NSZM method in the two optimization algorithms used. The numerical results demonstrate that the NSZM method enhances both algorithms equally.


---

## Body

## 1. Introduction

Distribution system planning (DSP) is a set of strategies that allows one to determine how many elements can be installed—and also where and when—in an electric network, to satisfy the growing demand for a specific time horizon [1, 2]. Traditionally, the DSP problem considers the installation of new electric elements and the upgrading of existing elements in medium- and low-voltage (MV/LV) networks (primary feeders, secondary circuits, substations, and distribution transformers) [3]. Currently, distributed generation (DG) has been integrated into distribution systems to improve network efficiency. Energy distribution companies can obtain major advantages through the integration of DG into distribution networks. These include network expansion deferral, reduction of operating costs (technical power losses of the networks), improved voltage levels, improvements in network reliability, and a decrease in environmental contamination [4–6].Distribution companies are facing new challenges in the DSP problem due to the integration of DGs, as wrong decisions in the planning process can lead to bad operating states in the network [7, 8]. Distribution companies therefore need reliable methodologies that can integrate the allocation and sizing of DGs in the DSP problem.From a mathematical point of view, the model of the DSP problem is a mixed-integer, nonlinear model (MINLP) because of the nonlinearity of the power flow equations and the integer variables that represent the investment decisions. Due to the combinatorial nature (NP-hard) of the problem, both voltage levels (MV/LV) have been solved separately, allowing for a reduced search space but not leading to a joint global solution [3, 9]. The number of papers regarding DSP for a medium voltage [10–15] is higher than for low voltage [16–19].Moreover, only a few papers consider both levels [20–23]. In [20], a heuristic method based on a biogeographical technique is proposed in order to carry out the integrated planning of MV/LV networks through the iterative construction of both networks. An MINLP mathematical model is proposed, which considers the minimization of the costs of MV/LV systems with regard to maintenance, reliability, operation, and investment (conductors, transformers, and substations). In [21], a general variable neighborhood search metaheuristic method is used to solve the integrated planning of MV/LV networks. An MINLP model is proposed, which aims at minimizing investment (substations, transformers, conductors, poles, and structures), operating, and reliability costs. The coupling constraints between MV/LV networks are represented by MV/LV modeling transformers through the power flow equations. In [22], an imperialist competitive algorithm is used to solve the integrated DSP problem of MV/LV networks, considering the location of DGs in the primary network. The authors proposed a single cost function, which comprises costs of MV/LV feeders, MV substations, DGs, and high-voltage substations. In [23], a decomposition approach and a Tabu search algorithm are proposed to solve the integrated DSP problem of MV/LV networks, considering the location of DGs in the LV network. The authors proposed an MINLP model, which considers the minimization of two objective functions that involve the costs of installing and upgrading the new and existing elements (MV/LV branches, distribution transformers, substations, and DGs) and the cost of the technical energy losses.The foregoing literature survey shows that only two papers consider the location and sizing of DG in the integrated DSP problem of MV/LV networks [22, 23]. Unlike the studies mentioned, this article presents a new methodology that solves the integrated DSP problem of MV/LV networks as a single network, where the location and sizing of DG are considered in LV networks. The installation of renewable DGs is considered in LV networks because some countries (such as Colombia) have regulated the DG connections at this voltage level. Moreover, distribution transformer (DT) overloads are common in LV networks with high penetration of DG [24]. This paper therefore models the electrical (transformer modeling) and the physical (route sharing between MV/LV networks and MV/LV optimal transformer location) relationships existing between MV/LV systems. Furthermore, considering both systems as a single network can lead easily to a joint, global solution.Classic optimization and metaheuristic methods have been used to solve the DSP problem [2]. Classic optimization guarantees an optimal global solution but requires excessive computational efforts for medium and large real systems because of the complexity of the DSP problem [3, 23]. Metaheuristics are powerful tools for solving complex optimization problems and are widely used in distribution planning problems [3, 21, 22]. As a consequence, metaheuristic techniques are used to solve the integrated DSP problem of MV/LV networks, considering the integration of DG in LV networks. The intensification and diversification processes are usually the bottleneck of metaheuristic techniques for solving complex problems, such as the proposed problem [9, 20]. To overcome such complexity, a new neighborhood search method based on the Zbus matrix (NSZM) is proposed to enhance the metaheuristic techniques used. This method uses the Zbus matrix to model MV/LV systems as a single network because it allows for the establishment of an electrical relationship between the primary and secondary nodes. As a consequence, a change in one network is reflected in the electrical characteristics of the other. Moreover, the Zbus matrix represents the relations between voltages and current injections in a system [25]. This electrical information can be used through sensitivity factors to reduce the solution space obtained by the neighborhood structure (branch exchange, upgrading existing elements, and the location of DGs and DTs).The Zbus matrix can be modified to reflect the changes produced by the neighborhood structure without any need to build it again [26]. Consequently, the NSZM method performs a more efficient and effective exploration of the solution space in order to find a joint, global solution that establishes a balanced benefit for integrated planning. In the specialized literature, the Zbus matrix has been used for different applications, such as load flow [27], reconfiguration of networks [28], and location of DGs [29, 30]. To the authors’ best knowledge, however, the Zbus matrix has not been used to solve the DSP problem.To verify the efficiency of the new NSZM method proposed, two well-known optimization algorithms (simulated annealing and iterated local search) were used. The algorithms SA and ILS, in conjunction with the NSZM method, are called in this paper SA-NSZM and ILS-NSZM, respectively. A decision was made to use these algorithms to solve the proposed integrated model due to the dimension and nonconvexity of the problem, the nature of some variables, and the mathematical complexity of the problem (NP-hard) and because these algorithms have been widely used, with satisfactory results, in solving planning problems with similar characteristics. Nevertheless, any other optimization algorithm could be used [2, 31, 32]. Moreover, there is a need to measure the performance behavior of the NSZM method in the two optimization algorithms used. A nonparametric statistical test, called the Wilcoxon signed ranks, is therefore used [33]. This test is a pairwise comparison, used in hypothesis-testing situations. The Wilcoxon signed ranks test is chosen because it is less influenced by sample size (number of case problems) than other nonparametric tests [33].The main contributions of this work are as follows:(i)
A new methodology, based on the impedance matrix, is proposed for solving the integrated DSP problem of MV/LV networks, considering the location and sizing of DG, where the problem is formulated as a single mathematical model. The Zbus matrix is used to represent both systems as a single network because it can lead easily to a joint, global solution. Furthermore, the methodology can be applied for real, large-size MV/LV distribution systems.(ii)
A new neighborhood search method based on the Zbus matrix (NSZM) is proposed for exploring the search space to find good-quality solutions. This method uses a defined neighborhood structure combined with sensitivity factors, based on the impedance matrix, in order to reduce the solution space for the SA and ILS algorithms. The Zbus matrix can be modified to reflect changes in the network produced by the neighborhood structure, without any need to rebuild it. As a consequence, the NSZM method enhances the search process of both algorithms. In addition, this method can be used by any metaheuristic technique that considers exploration and diversification of the search space to solve the DSP problem, considering the location and sizing of DG.(iii)
A new codification scheme is used to represent the proposed model. The difference from other papers within the literature on the subject is that all the installing decision variables are represented in a single vector. This allows for easy manipulation of the actions proposed by the SA-NSZM and the ILS-NSZM methodologies proposed.(iv)
To solve the allocation and sizing of DGs in the integrated DSP problem, the NSZM method uses a sensitivity analysis based on the impedance matrix. The analysis proposed in the NSZM method uses an equation, based on the impedance matrix, in order to measure easily the impact of the DG on the technical losses of both networks. The impedance matrix has been used for solving the allocation and sizing problems of DGs in distribution systems with good results [29, 30], but the Zbus matrix has not been used previously to solve the location and sizing of DGs together with the DSP problem.This paper is structured as follows: In the next section, the general mathematical formulation for the integrated DSP problem, considering allocation and sizing of DGs, is presented. Next, the proposed SA-NSZM and ILS-NSZM methodologies are described in Section 3. The proposed methodologies are tested on a distribution system, with the results being presented in Section 4. The results are analyzed and discussed in that section. Lastly, conclusions and future work are discussed in the final section.


## 2. Mathematical Formulation of the Problem

The impedance matrix is used in this mathematical model to represent both networks as a single system [34]. The problem is formulated as an MINLP model and both networks are represented by a single-phase equivalent. Equation (1) is the objective function that minimizes investment costs, and equations (2)–(22) provide the set of constraints. Equation (1) is the present value of nine terms.Terms 1 and 2 are the installation and upgrading costs of both new and existing primary feeders, respectively. Terms 3 and 4 are the installation and upgrading costs of both new and existing substations, respectively. Term 5 is the installation cost of new DTs. Terms 6 and 7 are the installation and upgrading costs of both new and existing secondary circuits, respectively. Term 8 is the installation cost of new DGs. Term 9 is the operating cost of the networks (technical energy losses in all elements).(1)min∑ij∈ΩNP∑p∈ΩTPCij,pNPγij,pNP+∑ij∈ΩEP∑p∈ΩTPCij,pEPγij,pEP+∑i∈ΩNSS∑s∈ΩTSSCsNSSγi,sNSS+∑i∈ΩESS∑s∈ΩTSSCsESSγi,sESS+∑i∈ΩNDT∑d∈ΩTDTCdNDTγi,dNDT+∑ij∈ΩNS∑c∈ΩTSCij,cNSγij,cNS+∑ij∈ΩES∑c∈ΩTSCij,cESγij,cES+∑i∈ΩNDG∑g∈ΩTDGCgNDGγi,gNDG+k1∑l=1NLk2lℝe∑i∈ΩSSSi,lSS−∑i∈ΩPNSi,lPDZIPi,lP+∑i∈ΩSNSi,lSDZIPi,lSΔTl.Equations (2) and (3) represent the power nodal balance given by Kirchhoff’s laws for MV and LV networks. Equation (4) and (5) represent the ZIP model for demands in MV and LV networks. Equations (6) and (7) use Ohm’s law to calculate the bus voltage using the Zbus matrix. Equation (8) uses Kirchhoff’s first law to calculate the bus currents in substations. Equations (9) and (10) calculate the currents in the primary feeders and secondary circuits, respectively. Equation (11) determines the injected power in each DT. Equation (12) determines the power generated from the installed DG.(2)Si,lSS=Si,lDT+Si,lPDZIPi,lP+Ui,lBUSIi,lBUS∗∀i∈ΩPN;∀l∈ΩNL,(3)Si,lDT+Si,lDG=Si,lSDZIPi,lS+Ui,lBUSIi,lBUS∗∀i∈ΩSN;∀l∈ΩNL,(4)ZIPi,lP=a0+a1Ui,lBUSUnomP+a2Ui,lBUSUnomP2∀i∈ΩPN;∀l∈ΩNL,(5)ZIPi,lS=a0+a1Ui,lBUSUnomS+a2Ui,lBUSUnomS2∀i∈ΩSN;∀l∈ΩNL,(6)Ui,lBUS=UnomP+∑j∈ΩNZijBUSIi,lBUS∀i∈ΩPN;∀l∈ΩNL,(7)Ui,lBUS=UnomS+∑j∈ΩNZijBUSIi,lBUS∀i∈ΩSN;∀l∈ΩNL,(8)Ii,lBUS=∑ij∈ΩPFIij,lP∀i∈ΩSS;∀l∈ΩNL,(9)Iij,lP=Ui,lBUS−Uj,lBUSZij,pP∑p∈ΩTPγij,pNP+γij,pEP∀ij∈ΩPF;∀l∈ΩNL,(10)Iij,lS=Ui,lBUS−Uj,lBUSZij,cS∑c∈ΩTSγij,cNS+γij,cES∀ij∈ΩSC;∀l∈ΩNL,(11)Si,lDT=Ui,lBUS∑ki∈ΩPFIki,lP−∑im∈ΩPFIim,lP∗∑d∈ΩTDTγi,dNDT∀i∈ΩNDT;∀l∈ΩNL,(12)Si,lDG=∑g∈ΩTDGγi,gNDGS¯g,lDG∀i∈ΩNDG;∀l∈ΩNL.Equations (13)–(17) represent the operative limits of the system. Equations (13) and (14) are the operative limits of the primary feeders and secondary circuits, respectively. Equations (15) and (16) are the operative limits of the substations and DTs, respectively. Equation (17) is the voltage limit of the nodes of the MV/LV networks. Equations (18)–(22) ensure, respectively, that only one type of wire, substation, DT, and DG can be installed in the same place.(13)0≤Iij,lP≤∑p∈ΩTPγij,pNP+γij,pEPI¯pP∀ij∈ΩPF;∀l∈ΩNL,(14)0≤Iij,lS≤∑c∈ΩTSγij,cNS+γij,cESI¯cS∀ij∈ΩSC;∀l∈ΩNL,(15)0≤Si,lSS≤∑s∈ΩTSSγi,sNSS+γi,sESSS¯sSS∀i∈ΩSS;∀l∈ΩNL,(16)0≤Si,lDT≤∑d∈ΩTDTγi,dNDTS¯dDT∀i∈ΩNDT;∀l∈ΩNL,(17)Uimin≤Ui,lBUS≤Uimax∀i∈ΩN;∀l∈ΩNL,(18)∑p∈ΩTPγij,pNP+γij,pEP≤1∀ij∈ΩPF,(19)∑c∈ΩTSγij,cNS+γij,cES≤1∀ij∈ΩSC,(20)∑s∈ΩTSSγi,sNSS+γi,sESS≤1∀i∈ΩSS,(21)∑d∈ΩTDTγi,dNDT≤1∀i∈ΩNDT,(22)∑g∈ΩTDGγi,gNDG≤1∀i∈ΩNDG.


## 3. Methodology Proposed

The MINLP model presented here is solved using two well-known, single solution-based metaheuristics (simulated annealing and iterative local search), combined with a new neighborhood search method based on the impedance matrix (NSZM). The proposed method uses the impedance matrix and the per-unit representation to model both networks as a single system. Furthermore, the NSZM method uses a defined neighborhood structure and equations, based on the impedance matrix, to reduce the search space in order to find good-quality solutions for the SA and ILS algorithms. Thus, the NSZM method enhances the search process of both metaheuristics.The proposed methodology starts with an initial, feasible solution for both networks (medium and low voltage), obtained from a constructive heuristic algorithm. Next, the Zbus matrix of the initial solution is built using a constructive algorithm. For more detailed information, the algorithm is explained in [26]. The NSZM method modifies the Zbus matrix of the current solution to reflect the changes made by the neighborhood structure; thus, the new solution obtained is represented by a new Zbus matrix without the need to build it again. The transition between solutions is controlled by each metaheuristic.The main aspects of the new NSZM method and the SA and ILS algorithms used in this paper will be discussed in the following.


### 3.1. Codification

The integrated DSP problem considering DG is represented using the new codification vector of Figure1. This codification vector uses integer numbers to represent the installing decision variables—where a number different to zero represents the capacity of an installed element, and a zero value shows that an element is not installed.Figure 1
Codification scheme for the integrated DSP problem considering DG.The vector used can be partitioned into five areas. The first area contains the locations and sizes of the existing and new primary feeders (sizen1+n2). The second area includes the locations and sizes of the existing and new secondary circuits (size n3+n4). The third area relates to the locations and capacities of the DTs (size n5). The fourth area relates to the locations and sizes of the existing and new substations (size n6+n7). The fifth area relates to the locations and capacities of the DGs (size n8).


### 3.2. Neighborhood Search Method Based on the Zbus Matrix (NSZM)

This method proposes a defined neighborhood structure and uses the Zbus matrix to obtain electrical and physical information about the solutions network. The information obtained from the Zbus matrix allows one to explore the solution space more efficiently and effectively. The criteria defined for the neighborhood structure use equation (23) or (25) to reduce the solution space. Equation (23) is used to obtain the bus voltages of the n nodes of the system, and, once these are calculated, other variables, such as the electric current in the branches and the power in the DTs, can be obtained easily. The symbols ° and ⊘ are the Hadamard product and division, respectively.(23)Un×1BUS=Un×1Slack+ℤn×nBUS⋅Sn×1D°ZIPn×1⊘Un×1Slack+Sn×1DG⊘Un×1Slack∗.Equation (25) presents the sensitivity factor used to reduce the search space for the allocation and sizing problem of DGs in the neighborhood structure. A demonstration of equation (25) is presented next. Equation (24) calculates the active and reactive losses of the system, using the impedance and current matrices, where ℤn×nBUS=Rn×nBUS+jXn×nBUS and In×1BUS=In×1BUS∠Θn×1 are complex numbers. Taking the partial derivative of the real part of equation (24) with respect to the bus current of node k yields equation (25):(24)SLosses=In×1BUST⋅ℤn×nBUS⋅In×1BUS∗.Equation (25) calculates the change in the technical losses of the system when the bus current injected in any bus k changes by the installation of different types of DG, where R⟶∗,kBUS is the column k of the real part of the impedance matrix, Θ⟶k is a vector containing the angle of the bus current in the node k, and I⟶BUS and Θ⟶ are vectors containing the bus current magnitude and angle in the n nodes of the system. Moreover, equation (25) can also be used to reduce the search space, when all buses from the system are considered as candidate buses for installing DGs.(25)∂PLosses∂IkBUS=2⟶1×n⋅I⟶n×1BUS°R⟶∗,kBUSn×1°cosΘ⟶kn×1−Θ⟶n×1.The impedance matrix is also used to obtain physical information about the current solution, such as the identification of terminal nodes, nodes upstream or downstream from a specific node, and send and receipt nodes. This information is obtained from equations (26)–(29).Equation (26) is used to find out if a specific node is a terminal node. The node k is a terminal node, if equation (26) is true for all i=1,2,…,n.(26)ZkkBUS≠ZkiBUSi=1,2,…,n;i≠k.Equation (27) is used in order to determine which nodes are upstream from a specific node. If this equation is true, the node i is upstream from the node k.(27)ZiiBUS=ZkiBUS.Equation (28) is used to know which nodes are downstream from a specific node. If this equation is true, the node i is downstream from the node k.(28)ZkkBUS=ZkiBUS.Equation (29) is used to identify which nodes are send and receipt nodes. If this equation is true, the node p is the receipt node and the node q is the send node. Otherwise, the node p is the send node and the node q is the receipt node.(29)ZppBUS>ZqqBUS.The physical information obtained from equations (26)–(29) is used in the neighborhood structure to make the changes proposed by the criteria. These equations are obtained by analyzing the constructive algorithm of the Zbus matrix proposed in [26]. Moreover, equations (26)–(29) can only be used in radial and single-phase networks that do not consider the shunt impedance of lines or buses in the Zbus matrix.


### 3.3. Neighborhood Structure

The neighborhood structure of the NSZM method uses the following criteria: the upgrading of existing substations, primary feeders and secondary circuits, installation of new substations, DTs and DGs, and branch exchange in both networks. Equations (23) and (25) are used to reduce the search space of the neighborhood. The neighborhood structure considered is explained in the following:(i)
Branch Exchange. The branch with maxΔUij is selected to connect the link. Next, the branch with minIij of the loop formed is selected to be removed. The currents of the loop are calculated using equation (23).(ii)
DT Installation. A DT is chosen randomly from the three DTs with the largest installed capacities. If there are DTs disconnected nearby, one is chosen at random to be installed, and some loads are transferred from the first DT. Otherwise, the DT to be installed is chosen at random from the disconnected DTs. Some loads need to be transferred to the new DT to comply with radiality in the secondary circuits. The loads are transferred from the nearest existing DT to the new one, by removing the circuit with minIij of the loop that is formed. These currents are calculated using equation (23). Finally, the operative limits of the DTs involved in the process are checked and corrected.(iii)
DT Removal. The DT to be removed is chosen randomly from the three DTs with the smallest installed capacities. The loads of this DT need to be transferred to another existing DT. The DT chosen for receiving these loads is the nearest one with less capacity. If there is more than one such, the DT that is less overloaded is chosen. If there is more than one DT that fulfills these conditions, the DT is chosen randomly from among them. Loads are transferred by removing the circuit with minIij of the loop formed. These currents are calculated using equation (23). Finally, the operative limits of the DTs involved in the process are checked and corrected.(iv)
Upgrading of Existing Branches and DTs. Two alternatives are considered. In the first case, an element is chosen randomly from the three most overloaded elements, and that element is replaced by one with a bigger size. In the second case, an element is chosen randomly from the three least overloaded elements, and that element is replaced by a smaller size. One of these alternatives is chosen randomly. Primary feeders, secondary circuits, and DTs are chosen for upgrading capacity, and one of these elements is selected randomly. Applying these criteria should ensure that the sections upstream do not present a smaller-size wire in the primary feeders and secondary circuits.(v)
Installation and Sizing of DGs. Four alternatives are considered: installation of DGs, removal of installed DGs, increased installed capacity of DGs, and decreased installed capacity of DGs. One of these alternatives is chosen randomly. For each alternative, equation (25) is used, and the bus that produces the minimum value is selected. For installing a DG, a random capacity of DG is selected from among the different types of DG, and equation (25) is used. For removing an installed DG, equation (25) is applied to all installed DGs. For increasing or decreasing the installed capacity of a DG, equation (25) is applied to all installed DGs, considering the new increase or decrease in the injected power of the DGs.The neighborhood structure uses equations (26)–(29) through the process to determine the physical network information of the solutions.


### 3.4. Modifications of the Zbus Matrix to Reflect Changes Made by the Neighborhood Structure

The Zbus matrix can be modified to reflect changes that are made to a network without the need to rebuild it. These changes may be the addition or removal of elements (feeders or circuits) for branch exchange, the addition or removal of elements (DTs or circuits) for DTs location, and changes in the impedance of elements (primary feeders, DTs, and secondary circuits) for upgrading or degrading the capacity. As a result, the Zbus matrix can be modified to reflect changes caused by the neighborhood structure (see Figure2).Figure 2
Zbus modification due to neighborhood criteria. (a) Branch exchange. (b) Upgrading or degrading capacity. (c) DT installation. (d) DT removal.
(a)(b)(c)(d)These changes, produced by the neighborhood structure, are reflected in the Zbus matrix by the addition of links in the network (see Figure2). To include the link in the Zbus matrix, some equations based on the relations between the current injections and the voltages in the system are used [26]. When these equations are applied, a new row and column are added to the Zbus matrix, and a fictitious node l is added. When a link is added, no new bus is added to the network. The dimensions of the Zbus matrix are therefore unchanged, but all the elements in the matrix must be recalculated to include the effect of the added link.For recalculating the elements of the Zbus matrix, Kron’s reduction is used. This procedure is based on the elimination of a node in a matrix where the independent variable of this node is equal to zero. The procedures to calculate the new Zbus matrix produced by the neighborhood structure are explained in the following.For branch exchange, the new element is added as a link between nodesp−q (see Figure 2(a)), and when it is added, Kron’s reduction is used. After this, the element k−m is removed. Between k−m, a link is added in parallel, the impedance of which is equal to the negative of the impedance of the element k−m to be removed (see Figure 2(a)). Kron’s reduction is then used.In per-unit representation, the DTs are treated as lines. As a consequence, the procedure described for branch exchange is the same for the DTs’ location (adding a DT or removing a DT). For the DTs’ installation, the new DT is added as a link between nodesB−e, and the secondary circuit c−d is removed (see Figure 2(c)). For the DTs’ removal, the new secondary circuit c−d is added as a link, and the DT between B−e is removed (see Figure 2(d)).The procedure for upgrading or degrading the capacity of an element is explained next. Betweenk−m, a link is added in parallel, the impedance of which is such that the equivalent impedance of the two elements in k−m is the desired value (see Figure 2(b)). The elements of the Zbus matrix are then recalculated. Equation (30) calculates the value of the impedance added in the link k−m, in order to obtain the new value of impedance from the old value.(30)zkmchange=zkmoldzkmnewzkmold−zkmnew.


### 3.5. Evaluation of the Configurations

The load flow is essential for knowing the operative conditions of each generated configuration. Hence, a load flow based on the Zbus Gauss method is used. For more detailed information, the load flow is explained in [26]. The load flow is used for calculating technical energy losses, to be summed in the objective function of equation (1) and for checking the feasibility of solutions.Unfeasible solutions are penalized in the objective function through the penalty costs of the respective violated constraint. The combination of objective function plus penalty costs is called the “fitness function”Ffit and has been widely used in [9, 10, 21, 23, 30]. Penalty costs are expressed in monetary units.


### 3.6. Single Solution-Based Metaheuristics

A single solution-based metaheuristic algorithm starts with an initial solution and moves away from it, describing a trajectory in the search space. Common single solution-based metaheuristic algorithms include simulated annealing, iterated greedy algorithm, Tabu search, GRASP, gravitational search, variable neighborhood search, variable neighborhood descent method, guided local search, iterated local search, and their variants.In this document, the SA and ILS algorithms are used and explained in the following.


#### 3.6.1. Simulated Annealing (SA) Algorithm

The SA algorithm is a stochastic optimization procedure that models the process of annealing in solids [35]. This technique simulates the heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The process ends when thermal equilibrium is achieved. The SA algorithm avoids getting trapped in a local minimum by accepting deteriorations of current solutions throughout the search process. Starting from an initial solution, the SA algorithm generates an amount of solution that will be evaluated to find the next new solution. A stochastic mechanism controls the transition between two solutions for the SA algorithm. The stochastic mechanism accepts a new solution depending on the values of the temperature level Tk and the objective function for the current and new solutions (Ffitcurrent and Ffitnew). If Ffitnew≤Ffitcurrent the new solution is accepted, and it replaces the old one. If Ffitnew≥Ffitcurrent, the new solution can be accepted or not with a certain probability given by equation (31). If this equation is true, then the new solution is accepted. Otherwise, the new solution is rejected.(31)expFfitcurrent−FfitnewTk>random0,1.Thermal equilibrium can be achieved by defining an adequate cooling scheme. The computational time, and the quality of the final solution, will depend on the selection of the cooling scheme and the tuning of its parameters. The cooling scheme employed is presented in equations (32) and (33). Equation (32) calculates the initial temperature T0, equation (33) calculates the new value for temperature Tk+1 at level k, and equation (34) calculates the new value for the number of transitions Nk+1 at level k. For more detailed information, the cooling scheme used is explained in [36].(32)T0=η−lnλFx0,(33)Tk+1=βTk,(34)Nk+1=αNk.The algorithm stops when the incumbent solution is not improved for a defined number of temperature levelsk, or when a fixed number of temperature levels k for the cooling process is achieved [32, 36].The search process for the SA algorithm is repeated until a certain number of transitionsNk is reached. Next, the temperature Tk is decreased, and the search process is repeated for a new temperature level. Hence, the probability of accepting deteriorating moves is high at the beginning of the search and gradually decreases. The cooling process ends when the stop criterion is met.


#### 3.6.2. Iterated Local Search (ILS) Algorithm

ILS is a metaheuristic, based on a simple idea [37]. The algorithm starts a local search from an initial solution. After a local optimum is discovered, the ILS algorithm generates a new starting solution for the next local search, by perturbing the local optimum found in the previous search. This is done in the expectation that the perturbation mechanism makes it possible to search for another, better, local minimum. Too weak perturbation might be insufficient to escape from the current local optimum. Conversely, too strong perturbation would render ILS akin to a multistart, local search with randomly generated starting solutions. The local search, combined with the perturbation mechanism, enables control over the trade-off between intensification and diversification.In the local search, solutions with lower costs are always accepted, while solutions with higher costs are always rejected. The local search process is repeated until a certain number of local iterations is completed. Afterwards, the new starting solution for the next local search is obtained by perturbing the local optimum discovered in the previous search. This process is repeated until a certain number of global iterations is completed.


### 3.7. General Methodology of the SA-NSZM

The procedure begins by obtaining the initial configuration, using a constructive heuristic algorithm. Next, the Zbus matrix is built for the initial solution, using the constructive algorithm proposed in [26]. Afterwards, the incumbent solution and codification vector are initialized.Next, the cooling scheme is calculated, as explained in Subsection3.6 The NSZM method is then used to obtain a new solution for the SA algorithm, as explained in Subsection 3.2. The new configuration represented by the Zbus matrix is evaluated, using the load flow and the fitness function. Subsequently, the stochastic mechanism of the SA algorithm controls the transition between solutions.The search process is repeated untilNk transitions are completed. For the next temperature level, Tk and Nk are recalculated by using equations (33) and (34). If the incumbent solution is improved in the search process, the incumbent is updated. Once a temperature level is completed, the initial Zbus matrix for the next temperature level is the incumbent Zbus matrix.The cooling process ends when the stop criterion is achieved. Figure3 shows the proposed SA-NSZM methodology.Figure 3
Flowchart of the SA-NSZM.


### 3.8. General Methodology of the ILS-NSZM

The procedure begins by obtaining the initial configuration, using a constructive heuristic algorithm. Next, the Zbus matrix is built for the initial solution, using the constructive algorithm proposed in [26]. Subsequently, the incumbent solution, the codification vector, and the factors (nd, ilocal, and jglobal) are initialized. Following this, the NSZM method is used to obtain a new solution for the ILS algorithm, as explained in Subsection. The new configuration represented by the Zbus matrix is evaluated, using the load flow and the fitness function.IfFfitnew≤Ffitcurrent, the new solution is accepted. Otherwise, the new solution is rejected. This local search process is repeated until ilocal local iterations are finished. Moreover, if the incumbent solution is improved in the local search process, the incumbent is updated. After a local search is finished, the incumbent solution is perturbed using the NSZM method nd times. These perturbations are repeated until jglobal global iterations are finished. The ILS-NSZM algorithm ends when the global search process has been repeated to the point at which jglobal global iterations are finished. Figure 4 shows the proposed ILS-NSZM methodology.Figure 4
Flowchart of the ILS-NSZM.


## 3.1. Codification

The integrated DSP problem considering DG is represented using the new codification vector of Figure1. This codification vector uses integer numbers to represent the installing decision variables—where a number different to zero represents the capacity of an installed element, and a zero value shows that an element is not installed.Figure 1
Codification scheme for the integrated DSP problem considering DG.The vector used can be partitioned into five areas. The first area contains the locations and sizes of the existing and new primary feeders (sizen1+n2). The second area includes the locations and sizes of the existing and new secondary circuits (size n3+n4). The third area relates to the locations and capacities of the DTs (size n5). The fourth area relates to the locations and sizes of the existing and new substations (size n6+n7). The fifth area relates to the locations and capacities of the DGs (size n8).


## 3.2. Neighborhood Search Method Based on the Zbus Matrix (NSZM)

This method proposes a defined neighborhood structure and uses the Zbus matrix to obtain electrical and physical information about the solutions network. The information obtained from the Zbus matrix allows one to explore the solution space more efficiently and effectively. The criteria defined for the neighborhood structure use equation (23) or (25) to reduce the solution space. Equation (23) is used to obtain the bus voltages of the n nodes of the system, and, once these are calculated, other variables, such as the electric current in the branches and the power in the DTs, can be obtained easily. The symbols ° and ⊘ are the Hadamard product and division, respectively.(23)Un×1BUS=Un×1Slack+ℤn×nBUS⋅Sn×1D°ZIPn×1⊘Un×1Slack+Sn×1DG⊘Un×1Slack∗.Equation (25) presents the sensitivity factor used to reduce the search space for the allocation and sizing problem of DGs in the neighborhood structure. A demonstration of equation (25) is presented next. Equation (24) calculates the active and reactive losses of the system, using the impedance and current matrices, where ℤn×nBUS=Rn×nBUS+jXn×nBUS and In×1BUS=In×1BUS∠Θn×1 are complex numbers. Taking the partial derivative of the real part of equation (24) with respect to the bus current of node k yields equation (25):(24)SLosses=In×1BUST⋅ℤn×nBUS⋅In×1BUS∗.Equation (25) calculates the change in the technical losses of the system when the bus current injected in any bus k changes by the installation of different types of DG, where R⟶∗,kBUS is the column k of the real part of the impedance matrix, Θ⟶k is a vector containing the angle of the bus current in the node k, and I⟶BUS and Θ⟶ are vectors containing the bus current magnitude and angle in the n nodes of the system. Moreover, equation (25) can also be used to reduce the search space, when all buses from the system are considered as candidate buses for installing DGs.(25)∂PLosses∂IkBUS=2⟶1×n⋅I⟶n×1BUS°R⟶∗,kBUSn×1°cosΘ⟶kn×1−Θ⟶n×1.The impedance matrix is also used to obtain physical information about the current solution, such as the identification of terminal nodes, nodes upstream or downstream from a specific node, and send and receipt nodes. This information is obtained from equations (26)–(29).Equation (26) is used to find out if a specific node is a terminal node. The node k is a terminal node, if equation (26) is true for all i=1,2,…,n.(26)ZkkBUS≠ZkiBUSi=1,2,…,n;i≠k.Equation (27) is used in order to determine which nodes are upstream from a specific node. If this equation is true, the node i is upstream from the node k.(27)ZiiBUS=ZkiBUS.Equation (28) is used to know which nodes are downstream from a specific node. If this equation is true, the node i is downstream from the node k.(28)ZkkBUS=ZkiBUS.Equation (29) is used to identify which nodes are send and receipt nodes. If this equation is true, the node p is the receipt node and the node q is the send node. Otherwise, the node p is the send node and the node q is the receipt node.(29)ZppBUS>ZqqBUS.The physical information obtained from equations (26)–(29) is used in the neighborhood structure to make the changes proposed by the criteria. These equations are obtained by analyzing the constructive algorithm of the Zbus matrix proposed in [26]. Moreover, equations (26)–(29) can only be used in radial and single-phase networks that do not consider the shunt impedance of lines or buses in the Zbus matrix.


## 3.3. Neighborhood Structure

The neighborhood structure of the NSZM method uses the following criteria: the upgrading of existing substations, primary feeders and secondary circuits, installation of new substations, DTs and DGs, and branch exchange in both networks. Equations (23) and (25) are used to reduce the search space of the neighborhood. The neighborhood structure considered is explained in the following:(i)
Branch Exchange. The branch with maxΔUij is selected to connect the link. Next, the branch with minIij of the loop formed is selected to be removed. The currents of the loop are calculated using equation (23).(ii)
DT Installation. A DT is chosen randomly from the three DTs with the largest installed capacities. If there are DTs disconnected nearby, one is chosen at random to be installed, and some loads are transferred from the first DT. Otherwise, the DT to be installed is chosen at random from the disconnected DTs. Some loads need to be transferred to the new DT to comply with radiality in the secondary circuits. The loads are transferred from the nearest existing DT to the new one, by removing the circuit with minIij of the loop that is formed. These currents are calculated using equation (23). Finally, the operative limits of the DTs involved in the process are checked and corrected.(iii)
DT Removal. The DT to be removed is chosen randomly from the three DTs with the smallest installed capacities. The loads of this DT need to be transferred to another existing DT. The DT chosen for receiving these loads is the nearest one with less capacity. If there is more than one such, the DT that is less overloaded is chosen. If there is more than one DT that fulfills these conditions, the DT is chosen randomly from among them. Loads are transferred by removing the circuit with minIij of the loop formed. These currents are calculated using equation (23). Finally, the operative limits of the DTs involved in the process are checked and corrected.(iv)
Upgrading of Existing Branches and DTs. Two alternatives are considered. In the first case, an element is chosen randomly from the three most overloaded elements, and that element is replaced by one with a bigger size. In the second case, an element is chosen randomly from the three least overloaded elements, and that element is replaced by a smaller size. One of these alternatives is chosen randomly. Primary feeders, secondary circuits, and DTs are chosen for upgrading capacity, and one of these elements is selected randomly. Applying these criteria should ensure that the sections upstream do not present a smaller-size wire in the primary feeders and secondary circuits.(v)
Installation and Sizing of DGs. Four alternatives are considered: installation of DGs, removal of installed DGs, increased installed capacity of DGs, and decreased installed capacity of DGs. One of these alternatives is chosen randomly. For each alternative, equation (25) is used, and the bus that produces the minimum value is selected. For installing a DG, a random capacity of DG is selected from among the different types of DG, and equation (25) is used. For removing an installed DG, equation (25) is applied to all installed DGs. For increasing or decreasing the installed capacity of a DG, equation (25) is applied to all installed DGs, considering the new increase or decrease in the injected power of the DGs.The neighborhood structure uses equations (26)–(29) through the process to determine the physical network information of the solutions.


## 3.4. Modifications of the Zbus Matrix to Reflect Changes Made by the Neighborhood Structure

The Zbus matrix can be modified to reflect changes that are made to a network without the need to rebuild it. These changes may be the addition or removal of elements (feeders or circuits) for branch exchange, the addition or removal of elements (DTs or circuits) for DTs location, and changes in the impedance of elements (primary feeders, DTs, and secondary circuits) for upgrading or degrading the capacity. As a result, the Zbus matrix can be modified to reflect changes caused by the neighborhood structure (see Figure2).Figure 2
Zbus modification due to neighborhood criteria. (a) Branch exchange. (b) Upgrading or degrading capacity. (c) DT installation. (d) DT removal.
(a)(b)(c)(d)These changes, produced by the neighborhood structure, are reflected in the Zbus matrix by the addition of links in the network (see Figure2). To include the link in the Zbus matrix, some equations based on the relations between the current injections and the voltages in the system are used [26]. When these equations are applied, a new row and column are added to the Zbus matrix, and a fictitious node l is added. When a link is added, no new bus is added to the network. The dimensions of the Zbus matrix are therefore unchanged, but all the elements in the matrix must be recalculated to include the effect of the added link.For recalculating the elements of the Zbus matrix, Kron’s reduction is used. This procedure is based on the elimination of a node in a matrix where the independent variable of this node is equal to zero. The procedures to calculate the new Zbus matrix produced by the neighborhood structure are explained in the following.For branch exchange, the new element is added as a link between nodesp−q (see Figure 2(a)), and when it is added, Kron’s reduction is used. After this, the element k−m is removed. Between k−m, a link is added in parallel, the impedance of which is equal to the negative of the impedance of the element k−m to be removed (see Figure 2(a)). Kron’s reduction is then used.In per-unit representation, the DTs are treated as lines. As a consequence, the procedure described for branch exchange is the same for the DTs’ location (adding a DT or removing a DT). For the DTs’ installation, the new DT is added as a link between nodesB−e, and the secondary circuit c−d is removed (see Figure 2(c)). For the DTs’ removal, the new secondary circuit c−d is added as a link, and the DT between B−e is removed (see Figure 2(d)).The procedure for upgrading or degrading the capacity of an element is explained next. Betweenk−m, a link is added in parallel, the impedance of which is such that the equivalent impedance of the two elements in k−m is the desired value (see Figure 2(b)). The elements of the Zbus matrix are then recalculated. Equation (30) calculates the value of the impedance added in the link k−m, in order to obtain the new value of impedance from the old value.(30)zkmchange=zkmoldzkmnewzkmold−zkmnew.


## 3.5. Evaluation of the Configurations

The load flow is essential for knowing the operative conditions of each generated configuration. Hence, a load flow based on the Zbus Gauss method is used. For more detailed information, the load flow is explained in [26]. The load flow is used for calculating technical energy losses, to be summed in the objective function of equation (1) and for checking the feasibility of solutions.Unfeasible solutions are penalized in the objective function through the penalty costs of the respective violated constraint. The combination of objective function plus penalty costs is called the “fitness function”Ffit and has been widely used in [9, 10, 21, 23, 30]. Penalty costs are expressed in monetary units.


## 3.6. Single Solution-Based Metaheuristics

A single solution-based metaheuristic algorithm starts with an initial solution and moves away from it, describing a trajectory in the search space. Common single solution-based metaheuristic algorithms include simulated annealing, iterated greedy algorithm, Tabu search, GRASP, gravitational search, variable neighborhood search, variable neighborhood descent method, guided local search, iterated local search, and their variants.In this document, the SA and ILS algorithms are used and explained in the following.


### 3.6.1. Simulated Annealing (SA) Algorithm

The SA algorithm is a stochastic optimization procedure that models the process of annealing in solids [35]. This technique simulates the heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The process ends when thermal equilibrium is achieved. The SA algorithm avoids getting trapped in a local minimum by accepting deteriorations of current solutions throughout the search process. Starting from an initial solution, the SA algorithm generates an amount of solution that will be evaluated to find the next new solution. A stochastic mechanism controls the transition between two solutions for the SA algorithm. The stochastic mechanism accepts a new solution depending on the values of the temperature level Tk and the objective function for the current and new solutions (Ffitcurrent and Ffitnew). If Ffitnew≤Ffitcurrent the new solution is accepted, and it replaces the old one. If Ffitnew≥Ffitcurrent, the new solution can be accepted or not with a certain probability given by equation (31). If this equation is true, then the new solution is accepted. Otherwise, the new solution is rejected.(31)expFfitcurrent−FfitnewTk>random0,1.Thermal equilibrium can be achieved by defining an adequate cooling scheme. The computational time, and the quality of the final solution, will depend on the selection of the cooling scheme and the tuning of its parameters. The cooling scheme employed is presented in equations (32) and (33). Equation (32) calculates the initial temperature T0, equation (33) calculates the new value for temperature Tk+1 at level k, and equation (34) calculates the new value for the number of transitions Nk+1 at level k. For more detailed information, the cooling scheme used is explained in [36].(32)T0=η−lnλFx0,(33)Tk+1=βTk,(34)Nk+1=αNk.The algorithm stops when the incumbent solution is not improved for a defined number of temperature levelsk, or when a fixed number of temperature levels k for the cooling process is achieved [32, 36].The search process for the SA algorithm is repeated until a certain number of transitionsNk is reached. Next, the temperature Tk is decreased, and the search process is repeated for a new temperature level. Hence, the probability of accepting deteriorating moves is high at the beginning of the search and gradually decreases. The cooling process ends when the stop criterion is met.


### 3.6.2. Iterated Local Search (ILS) Algorithm

ILS is a metaheuristic, based on a simple idea [37]. The algorithm starts a local search from an initial solution. After a local optimum is discovered, the ILS algorithm generates a new starting solution for the next local search, by perturbing the local optimum found in the previous search. This is done in the expectation that the perturbation mechanism makes it possible to search for another, better, local minimum. Too weak perturbation might be insufficient to escape from the current local optimum. Conversely, too strong perturbation would render ILS akin to a multistart, local search with randomly generated starting solutions. The local search, combined with the perturbation mechanism, enables control over the trade-off between intensification and diversification.In the local search, solutions with lower costs are always accepted, while solutions with higher costs are always rejected. The local search process is repeated until a certain number of local iterations is completed. Afterwards, the new starting solution for the next local search is obtained by perturbing the local optimum discovered in the previous search. This process is repeated until a certain number of global iterations is completed.


## 3.6.1. Simulated Annealing (SA) Algorithm

The SA algorithm is a stochastic optimization procedure that models the process of annealing in solids [35]. This technique simulates the heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The process ends when thermal equilibrium is achieved. The SA algorithm avoids getting trapped in a local minimum by accepting deteriorations of current solutions throughout the search process. Starting from an initial solution, the SA algorithm generates an amount of solution that will be evaluated to find the next new solution. A stochastic mechanism controls the transition between two solutions for the SA algorithm. The stochastic mechanism accepts a new solution depending on the values of the temperature level Tk and the objective function for the current and new solutions (Ffitcurrent and Ffitnew). If Ffitnew≤Ffitcurrent the new solution is accepted, and it replaces the old one. If Ffitnew≥Ffitcurrent, the new solution can be accepted or not with a certain probability given by equation (31). If this equation is true, then the new solution is accepted. Otherwise, the new solution is rejected.(31)expFfitcurrent−FfitnewTk>random0,1.Thermal equilibrium can be achieved by defining an adequate cooling scheme. The computational time, and the quality of the final solution, will depend on the selection of the cooling scheme and the tuning of its parameters. The cooling scheme employed is presented in equations (32) and (33). Equation (32) calculates the initial temperature T0, equation (33) calculates the new value for temperature Tk+1 at level k, and equation (34) calculates the new value for the number of transitions Nk+1 at level k. For more detailed information, the cooling scheme used is explained in [36].(32)T0=η−lnλFx0,(33)Tk+1=βTk,(34)Nk+1=αNk.The algorithm stops when the incumbent solution is not improved for a defined number of temperature levelsk, or when a fixed number of temperature levels k for the cooling process is achieved [32, 36].The search process for the SA algorithm is repeated until a certain number of transitionsNk is reached. Next, the temperature Tk is decreased, and the search process is repeated for a new temperature level. Hence, the probability of accepting deteriorating moves is high at the beginning of the search and gradually decreases. The cooling process ends when the stop criterion is met.


## 3.6.2. Iterated Local Search (ILS) Algorithm

ILS is a metaheuristic, based on a simple idea [37]. The algorithm starts a local search from an initial solution. After a local optimum is discovered, the ILS algorithm generates a new starting solution for the next local search, by perturbing the local optimum found in the previous search. This is done in the expectation that the perturbation mechanism makes it possible to search for another, better, local minimum. Too weak perturbation might be insufficient to escape from the current local optimum. Conversely, too strong perturbation would render ILS akin to a multistart, local search with randomly generated starting solutions. The local search, combined with the perturbation mechanism, enables control over the trade-off between intensification and diversification.In the local search, solutions with lower costs are always accepted, while solutions with higher costs are always rejected. The local search process is repeated until a certain number of local iterations is completed. Afterwards, the new starting solution for the next local search is obtained by perturbing the local optimum discovered in the previous search. This process is repeated until a certain number of global iterations is completed.


## 3.7. General Methodology of the SA-NSZM

The procedure begins by obtaining the initial configuration, using a constructive heuristic algorithm. Next, the Zbus matrix is built for the initial solution, using the constructive algorithm proposed in [26]. Afterwards, the incumbent solution and codification vector are initialized.Next, the cooling scheme is calculated, as explained in Subsection3.6 The NSZM method is then used to obtain a new solution for the SA algorithm, as explained in Subsection 3.2. The new configuration represented by the Zbus matrix is evaluated, using the load flow and the fitness function. Subsequently, the stochastic mechanism of the SA algorithm controls the transition between solutions.The search process is repeated untilNk transitions are completed. For the next temperature level, Tk and Nk are recalculated by using equations (33) and (34). If the incumbent solution is improved in the search process, the incumbent is updated. Once a temperature level is completed, the initial Zbus matrix for the next temperature level is the incumbent Zbus matrix.The cooling process ends when the stop criterion is achieved. Figure3 shows the proposed SA-NSZM methodology.Figure 3
Flowchart of the SA-NSZM.


## 3.8. General Methodology of the ILS-NSZM

The procedure begins by obtaining the initial configuration, using a constructive heuristic algorithm. Next, the Zbus matrix is built for the initial solution, using the constructive algorithm proposed in [26]. Subsequently, the incumbent solution, the codification vector, and the factors (nd, ilocal, and jglobal) are initialized. Following this, the NSZM method is used to obtain a new solution for the ILS algorithm, as explained in Subsection. The new configuration represented by the Zbus matrix is evaluated, using the load flow and the fitness function.IfFfitnew≤Ffitcurrent, the new solution is accepted. Otherwise, the new solution is rejected. This local search process is repeated until ilocal local iterations are finished. Moreover, if the incumbent solution is improved in the local search process, the incumbent is updated. After a local search is finished, the incumbent solution is perturbed using the NSZM method nd times. These perturbations are repeated until jglobal global iterations are finished. The ILS-NSZM algorithm ends when the global search process has been repeated to the point at which jglobal global iterations are finished. Figure 4 shows the proposed ILS-NSZM methodology.Figure 4
Flowchart of the ILS-NSZM.


## 4. Results and Methodology Validation

### 4.1. Description of the System

The SA-NSZM and ILS-NSZM methodologies were tested on the real, medium-low voltage distribution system presented in [23]. The primary distribution network has 48 existing buses, 1 existing substation (type 2), and 51 existing feeders (type 1). In the MV network, 60 new feeders and 1 new substation can be installed to supply new demands. In the LV network, 33 new DTs, 15 new DGs, and 147 new secondary circuits are proposed to supply the 138 new secondary-demand nodes. In addition, 5 types of substations, 6 types of wires for primary and secondary, 8 types of DTs, and 4 types of DGs are considered. The full system database is presented in [38].


### 4.2. Parameters of the SA-NSZM and ILS-NSZM Algorithms

The parameters of both algorithms were adjusted by a trial and error process [35, 37]. The factors used for the cooling scheme of the SA-NSZM algorithm were η=1%, λ=5%, α=1.05, and β=0.95. The cost of the initial solution Fx0 for Cases A and B was 1.85 million USD, and for Cases C and D it was 2.03 million USD. The initial number of transitions N0 for Cases A and C was 293, and for Cases B and D it was 308. N0 was defined as the number of variables of each problem. The initial temperature T0, according to equation (32), becomes 6175.6 for Cases A and B and 6798.3 for Cases C and D. The stopping criteria used were 15 temperature levels if the incumbent solution is not improved and 60 temperature levels for the cooling process. According to equation (33), the final temperature Tf for 60 levels becomes 284.50 for Cases A and B and 313.20 for Cases C and D.The factors used for the ILS-NSZM algorithm werend=3, ilocal=100, and jglobal=100.


### 4.3. Validation of the Single-Phase Equivalent

To validate the proposed methodologies, the results are compared with those reported in [23], using the same design features proposed by them. One must clarify that, in [23], the primary and secondary networks are modeled as a single-phase and three-phase equivalents, respectively. In this paper, both networks are modeled as a single network by a single-phase equivalent. Before testing the proposed methodologies, therefore, the final configuration presented in [23] is evaluated by using the single-phase equivalent proposed in this paper.In [23], three study cases are proposed, with Case 2 involving integrated planning without considering DG. The published results for the objective function, and the cost of the energy losses in the MV and LV networks for that case, are 1.789, 0.337, and 0.184 (million USD), respectively. The results obtained by using the single-phase equivalent proposed in this paper, for that case, are 1.773, 0.335, and 0.169 (million USD). The error between the two objective functions can be seen to be less than 1%. As a consequence, the results obtained by modeling both systems as a single network by a single-phase equivalent are reliable and approximate to those reported in [23].


### 4.4. Results

The SA-NSZM and ILS-NSZM methodologies were implemented using Matlab (2017b), in an Intel® Core i5-4460S 12 GB RAM PC. To demonstrate the effectiveness and robustness of the proposed methodologies, four cases were studied: (A) integrated planning without DG (Case 2 in [23]), (B) integrated planning with DG in the LV network (Case 3 in [23]), (C) integrated planning without DG, considering the copper losses in DTs, and (D) integrated planning with DG in the LV network, considering the copper losses in DTs.Table1 shows the average, best, and worst values found for the SA-NSZM and ILS-NSZM algorithms over 30 independent runs for the four cases. Additionally, this table shows the best solutions reported in [23], for Cases A and B. In [23], a Tabu search (TS) algorithm with a bilevel approach is used to solve the integrated DSP problem, considering DG.Table 1
Objective functions in millions of USD.
CasesRef. [23]SA-NSZMILS-NSZMBest reportedMeanBestWorstMeanBestWorstA1.7891.6561.6271.6961.6561.6261.689B1.7371.6111.5731.6461.6231.5741.659C—1.8461.8081.8841.8511.8161.884D—1.7571.7151.7931.7481.7191.793Table1 shows that even the worst solutions found with the SA-NSZM and ILS-NSZM algorithms, for Cases A and B, are better than the best ones reported in [23]. Based on this premise, it is evident that the NSZM method enhances the SA and ILS metaheuristics used for solving the integrated DSP problem, considering DG.Nevertheless, one needs to evaluate which of the two algorithms is better, or whether the NSZM method enhances both algorithms equally. A nonparametric statistical test called the Wilcoxon signed ranks is therefore used [33]. This test is a pairwise comparison, employed in hypothesis-testing situations. The nonparametric test uses the mean of the runs for each case problem of two algorithms, unlike a parametric test in which the null hypothesis is the equivalent of the median [39]. The Wilcoxon signed ranks test uses the average of the best objective function values of Table 1. The Wilcoxon signed ranks test is chosen as it is less influenced by sample size (number of case problems) than other nonparametric tests [33].The Wilcoxon test uses the difference between the performance scores of the two algorithms oni th out of n problems. The differences are ranked according to their absolute values; in case of ties, average ranks are used. R+ is the sum of ranks for the problems in which the first algorithm outperformed the second, and R− is the sum of ranks for the opposite. The critical value is the smaller of the two sums. If this value is less than or equal to the value of the distribution of Wilcoxon for n degrees of freedom (see Table B.12 in [40]), then the null hypothesis of equality of means is rejected. This will mean that the SA-NSZM algorithm outperforms the ILS-NSZM algorithm, with the associated p value. A more detailed explanation of this test can be found in [33].Table2 shows the R+, R−, and p value computed for the Wilcoxon test concerning the SA-NSZM algorithm (the p value has been computed by using SPSS). The results obtained by way of the Wilcoxon test show that the null hypothesis is not rejected. The differences observed between the algorithms are therefore explicable by chance alone. This means that the NSZM method enhances both algorithms equally and can be applied to enhance any single, solution-based metaheuristic in order to solve problems related to the expansion planning of power systems, considering the installation of DGs.Table 2
Wilcoxon signed ranks test results.
ComparisonR+R−p valueSA-NSZM vs ILS-NSZM640.715A more detailed comparison of the best results obtained is shown in Tables3 and 4 in terms of present value, where the term ETL means the cost of energy losses. Table 3 shows a comparison with the results published in [23] (Cases A and B). Table 4 shows more details of the results obtained for Cases C and D because, unlike Cases A and B, these consider the cost of energy losses in the DTs. Tables 3 and 4 show the best solutions obtained using the proposed algorithms, where the SA-NSZM discovered better solutions for Cases B, C, and D; and the ILS-NSZM discovered a better solution for Case A. Despite SA-NSZM identifying the best solutions for three cases, the errors between the best solutions for both algorithms are less than 0.5%. Table 3 shows that the SA-NSZM and ILS-NSZM methodologies are faster than the TS algorithm with the bilevel approach, used in [23]. Moreover, Table 3 shows that ILS-NSZM is the faster algorithm. Although the Wilcoxon test says that the SA-NSZM and ILS-NSZM algorithms show the same behavior, the CPU times obtained from the ILS-NSZM algorithm make it more attractive. Thus, we suggest that use is made of the ILS-NSZM because of the combination of the high-quality solutions obtained and the lower computational times achieved.Table 3
Comparison with results reported in [23] in millions of USD, Cases A and B.
CostDescriptionReported in [23]SA-NSZMILS-NSZMCase ACase BCase ACase BCase ACase BFixedSubstations0.3360.3360.3360.3360.3360.336MV feeders0.4100.3970.3790.3480.3710.341LV circuits0.2550.2540.2790.2670.2640.262DT0.2650.2250.2540.2160.2380.218DG—0.032—0.026—0.023Total1.2671.2461.2491.1931.2091.179VariableETL in MV0.3370.3070.2090.1970.2090.199ETL in LV0.1840.1830.1690.1830.2070.196Total0.5210.4900.3780.3800.4160.395Total cost1.7891.7371.6271.5731.6261.574CPU time (s)30603360875960114111Table 4
Comparison of best solutions in millions of USD—cases C and D.
CostDescriptionSA-NSZMILS-NSZMCases CCase DCase CCase DFixedSubstations0.3360.3360.3360.336MV feeders0.3790.3420.3850.327LV circuits0.2770.2700.2710.277DT0.2690.2250.2650.223DG—0.033—0.036Total1.2611.2061.2561.199VariableETL in MV0.2090.1910.1970.191ETL in LV0.1690.1910.1960.207ETL in DTs0.1670.1280.1670.122Total0.5460.5100.5600.520Total cost1.8081.7151.8161.719CPU time (s)65783792106Table3 shows that the SA-NSZM and ILS-NSZM algorithms discovered lower operating costs, using less investment cost, for DGs than those reported in [23]. As a result, one can see that it is desirable to use the Zbus matrix in the analysis of the impact of DGs on the technical losses of both networks, in order to solve the location and sizing problems of DGs in the DSP problem. Furthermore, the results obtained from Tables 3 and 4 show that DG penetration in the distribution system decreases investment and operating costs. The development of efficient and effective methodologies that can solve the DSP problem, considering DG integration, is therefore highly desirable and necessary.In Table3, it can be seen that the use of the NSZM method proposed in this paper allows the SA and ILS algorithms to discover lower operating costs in MV and LV networks because the Zbus matrix represents both systems as a single network. As a consequence, the Zbus matrix allows the establishment of an electrical relation between MV and LV networks; thus, the circulation of power flow between both networks can reach a global equilibrium.Finally, Table3 shows that the total costs obtained with the SA-NSZM and ILS-NSZM algorithms are approximately 9% better than those reported in [23]. As a consequence, new global solutions are found in this paper. The best configurations obtained in this paper for the four studied cases are shown in Figures 5–12. For all cases, the solutions found have different configurations and are feasible.Figure 5
Case A: primary network.Figure 6
Case A: secondary networks.Figure 7
Case B: primary network.Figure 8
Case B: secondary networks.Figure 9
Case C: primary network.Figure 10
Case C: secondary networks.Figure 11
Case D: primary network.Figure 12
Case D: secondary networks.Tables3 and 4 show that the cost of the energy losses in DTs has a considerable impact on objective function and it is inadvisable to ignore them. Moreover, Tables 3 and 4 show that considering both systems as a single network is faster than solving both networks separately, as was proposed in [23].The LV buses with DGs installed for Cases B and D are presented in Table5. One can see that the locations and sizes obtained for DGs differ between cases.Table 5
Locations and types of DG installed for Cases B and D.
Case BLV node2338438895103110130131136DG type1211111111Case DLV node233843718895103110130131136DG type22111112121


## 4.1. Description of the System

The SA-NSZM and ILS-NSZM methodologies were tested on the real, medium-low voltage distribution system presented in [23]. The primary distribution network has 48 existing buses, 1 existing substation (type 2), and 51 existing feeders (type 1). In the MV network, 60 new feeders and 1 new substation can be installed to supply new demands. In the LV network, 33 new DTs, 15 new DGs, and 147 new secondary circuits are proposed to supply the 138 new secondary-demand nodes. In addition, 5 types of substations, 6 types of wires for primary and secondary, 8 types of DTs, and 4 types of DGs are considered. The full system database is presented in [38].


## 4.2. Parameters of the SA-NSZM and ILS-NSZM Algorithms

The parameters of both algorithms were adjusted by a trial and error process [35, 37]. The factors used for the cooling scheme of the SA-NSZM algorithm were η=1%, λ=5%, α=1.05, and β=0.95. The cost of the initial solution Fx0 for Cases A and B was 1.85 million USD, and for Cases C and D it was 2.03 million USD. The initial number of transitions N0 for Cases A and C was 293, and for Cases B and D it was 308. N0 was defined as the number of variables of each problem. The initial temperature T0, according to equation (32), becomes 6175.6 for Cases A and B and 6798.3 for Cases C and D. The stopping criteria used were 15 temperature levels if the incumbent solution is not improved and 60 temperature levels for the cooling process. According to equation (33), the final temperature Tf for 60 levels becomes 284.50 for Cases A and B and 313.20 for Cases C and D.The factors used for the ILS-NSZM algorithm werend=3, ilocal=100, and jglobal=100.


## 4.3. Validation of the Single-Phase Equivalent

To validate the proposed methodologies, the results are compared with those reported in [23], using the same design features proposed by them. One must clarify that, in [23], the primary and secondary networks are modeled as a single-phase and three-phase equivalents, respectively. In this paper, both networks are modeled as a single network by a single-phase equivalent. Before testing the proposed methodologies, therefore, the final configuration presented in [23] is evaluated by using the single-phase equivalent proposed in this paper.In [23], three study cases are proposed, with Case 2 involving integrated planning without considering DG. The published results for the objective function, and the cost of the energy losses in the MV and LV networks for that case, are 1.789, 0.337, and 0.184 (million USD), respectively. The results obtained by using the single-phase equivalent proposed in this paper, for that case, are 1.773, 0.335, and 0.169 (million USD). The error between the two objective functions can be seen to be less than 1%. As a consequence, the results obtained by modeling both systems as a single network by a single-phase equivalent are reliable and approximate to those reported in [23].


## 4.4. Results

The SA-NSZM and ILS-NSZM methodologies were implemented using Matlab (2017b), in an Intel® Core i5-4460S 12 GB RAM PC. To demonstrate the effectiveness and robustness of the proposed methodologies, four cases were studied: (A) integrated planning without DG (Case 2 in [23]), (B) integrated planning with DG in the LV network (Case 3 in [23]), (C) integrated planning without DG, considering the copper losses in DTs, and (D) integrated planning with DG in the LV network, considering the copper losses in DTs.Table1 shows the average, best, and worst values found for the SA-NSZM and ILS-NSZM algorithms over 30 independent runs for the four cases. Additionally, this table shows the best solutions reported in [23], for Cases A and B. In [23], a Tabu search (TS) algorithm with a bilevel approach is used to solve the integrated DSP problem, considering DG.Table 1
Objective functions in millions of USD.
CasesRef. [23]SA-NSZMILS-NSZMBest reportedMeanBestWorstMeanBestWorstA1.7891.6561.6271.6961.6561.6261.689B1.7371.6111.5731.6461.6231.5741.659C—1.8461.8081.8841.8511.8161.884D—1.7571.7151.7931.7481.7191.793Table1 shows that even the worst solutions found with the SA-NSZM and ILS-NSZM algorithms, for Cases A and B, are better than the best ones reported in [23]. Based on this premise, it is evident that the NSZM method enhances the SA and ILS metaheuristics used for solving the integrated DSP problem, considering DG.Nevertheless, one needs to evaluate which of the two algorithms is better, or whether the NSZM method enhances both algorithms equally. A nonparametric statistical test called the Wilcoxon signed ranks is therefore used [33]. This test is a pairwise comparison, employed in hypothesis-testing situations. The nonparametric test uses the mean of the runs for each case problem of two algorithms, unlike a parametric test in which the null hypothesis is the equivalent of the median [39]. The Wilcoxon signed ranks test uses the average of the best objective function values of Table 1. The Wilcoxon signed ranks test is chosen as it is less influenced by sample size (number of case problems) than other nonparametric tests [33].The Wilcoxon test uses the difference between the performance scores of the two algorithms oni th out of n problems. The differences are ranked according to their absolute values; in case of ties, average ranks are used. R+ is the sum of ranks for the problems in which the first algorithm outperformed the second, and R− is the sum of ranks for the opposite. The critical value is the smaller of the two sums. If this value is less than or equal to the value of the distribution of Wilcoxon for n degrees of freedom (see Table B.12 in [40]), then the null hypothesis of equality of means is rejected. This will mean that the SA-NSZM algorithm outperforms the ILS-NSZM algorithm, with the associated p value. A more detailed explanation of this test can be found in [33].Table2 shows the R+, R−, and p value computed for the Wilcoxon test concerning the SA-NSZM algorithm (the p value has been computed by using SPSS). The results obtained by way of the Wilcoxon test show that the null hypothesis is not rejected. The differences observed between the algorithms are therefore explicable by chance alone. This means that the NSZM method enhances both algorithms equally and can be applied to enhance any single, solution-based metaheuristic in order to solve problems related to the expansion planning of power systems, considering the installation of DGs.Table 2
Wilcoxon signed ranks test results.
ComparisonR+R−p valueSA-NSZM vs ILS-NSZM640.715A more detailed comparison of the best results obtained is shown in Tables3 and 4 in terms of present value, where the term ETL means the cost of energy losses. Table 3 shows a comparison with the results published in [23] (Cases A and B). Table 4 shows more details of the results obtained for Cases C and D because, unlike Cases A and B, these consider the cost of energy losses in the DTs. Tables 3 and 4 show the best solutions obtained using the proposed algorithms, where the SA-NSZM discovered better solutions for Cases B, C, and D; and the ILS-NSZM discovered a better solution for Case A. Despite SA-NSZM identifying the best solutions for three cases, the errors between the best solutions for both algorithms are less than 0.5%. Table 3 shows that the SA-NSZM and ILS-NSZM methodologies are faster than the TS algorithm with the bilevel approach, used in [23]. Moreover, Table 3 shows that ILS-NSZM is the faster algorithm. Although the Wilcoxon test says that the SA-NSZM and ILS-NSZM algorithms show the same behavior, the CPU times obtained from the ILS-NSZM algorithm make it more attractive. Thus, we suggest that use is made of the ILS-NSZM because of the combination of the high-quality solutions obtained and the lower computational times achieved.Table 3
Comparison with results reported in [23] in millions of USD, Cases A and B.
CostDescriptionReported in [23]SA-NSZMILS-NSZMCase ACase BCase ACase BCase ACase BFixedSubstations0.3360.3360.3360.3360.3360.336MV feeders0.4100.3970.3790.3480.3710.341LV circuits0.2550.2540.2790.2670.2640.262DT0.2650.2250.2540.2160.2380.218DG—0.032—0.026—0.023Total1.2671.2461.2491.1931.2091.179VariableETL in MV0.3370.3070.2090.1970.2090.199ETL in LV0.1840.1830.1690.1830.2070.196Total0.5210.4900.3780.3800.4160.395Total cost1.7891.7371.6271.5731.6261.574CPU time (s)30603360875960114111Table 4
Comparison of best solutions in millions of USD—cases C and D.
CostDescriptionSA-NSZMILS-NSZMCases CCase DCase CCase DFixedSubstations0.3360.3360.3360.336MV feeders0.3790.3420.3850.327LV circuits0.2770.2700.2710.277DT0.2690.2250.2650.223DG—0.033—0.036Total1.2611.2061.2561.199VariableETL in MV0.2090.1910.1970.191ETL in LV0.1690.1910.1960.207ETL in DTs0.1670.1280.1670.122Total0.5460.5100.5600.520Total cost1.8081.7151.8161.719CPU time (s)65783792106Table3 shows that the SA-NSZM and ILS-NSZM algorithms discovered lower operating costs, using less investment cost, for DGs than those reported in [23]. As a result, one can see that it is desirable to use the Zbus matrix in the analysis of the impact of DGs on the technical losses of both networks, in order to solve the location and sizing problems of DGs in the DSP problem. Furthermore, the results obtained from Tables 3 and 4 show that DG penetration in the distribution system decreases investment and operating costs. The development of efficient and effective methodologies that can solve the DSP problem, considering DG integration, is therefore highly desirable and necessary.In Table3, it can be seen that the use of the NSZM method proposed in this paper allows the SA and ILS algorithms to discover lower operating costs in MV and LV networks because the Zbus matrix represents both systems as a single network. As a consequence, the Zbus matrix allows the establishment of an electrical relation between MV and LV networks; thus, the circulation of power flow between both networks can reach a global equilibrium.Finally, Table3 shows that the total costs obtained with the SA-NSZM and ILS-NSZM algorithms are approximately 9% better than those reported in [23]. As a consequence, new global solutions are found in this paper. The best configurations obtained in this paper for the four studied cases are shown in Figures 5–12. For all cases, the solutions found have different configurations and are feasible.Figure 5
Case A: primary network.Figure 6
Case A: secondary networks.Figure 7
Case B: primary network.Figure 8
Case B: secondary networks.Figure 9
Case C: primary network.Figure 10
Case C: secondary networks.Figure 11
Case D: primary network.Figure 12
Case D: secondary networks.Tables3 and 4 show that the cost of the energy losses in DTs has a considerable impact on objective function and it is inadvisable to ignore them. Moreover, Tables 3 and 4 show that considering both systems as a single network is faster than solving both networks separately, as was proposed in [23].The LV buses with DGs installed for Cases B and D are presented in Table5. One can see that the locations and sizes obtained for DGs differ between cases.Table 5
Locations and types of DG installed for Cases B and D.
Case BLV node2338438895103110130131136DG type1211111111Case DLV node233843718895103110130131136DG type22111112121


## 5. Conclusions

This paper presents a new methodology that uses a novel neighborhood search method, based on the Zbus matrix (NSZM) with two well-known, single solution-based metaheuristics (simulated annealing algorithm and iterated local search algorithm), in order to solve the integrated planning of MV/LV distribution systems, considering the allocation and sizing of DGs. The NSZM method uses a defined, neighborhood structure, combined with sensitivity factors based on the impedance matrix, to identify attractive solutions swiftly for the SA and ILS algorithms. The proposed methodologies have been validated and tested using a real MV/LV distribution system, taken from the specialized literature. The results obtained with the SA-NSZM and ILS-NSZM algorithms are better than those reported in previous literature. As a result, we demonstrate that the NSZM method can be used to explore the solution space more efficiently and effectively for the SA and ILS algorithms, enhancing both algorithms.In this paper, the Wilcoxon signed ranks test is used to compare the performance behavior of the SA-NSZM and ILS-NSZM algorithms. The results obtained with the test show that the differences observed between the two algorithms are explicable by chance alone. This means that the NSZM method enhances both algorithms equally and can be applied to enhance any single, solution-based metaheuristic to solve the DSP problem, considering the installation of DGs. Although the behavior of the SA-NSZM and ILS-NSZM algorithms is almost identical as regards the quality of the solutions obtained, the CPU times obtained from the ILS-NSZM algorithm make it the more attractive option.The Zbus matrix was used to model the MV/LV systems as a single network because it allows for the establishment of an electrical relationship between primary and secondary nodes. A change in one network is therefore reflected in the electrical characteristics of the other one. The results obtained show that the NSZM method proposed in this paper allows the SA and ILS algorithms to discover lower operating costs in the MV and LV networks, as the power flow circulation between both networks can reach a global, optimal equilibrium.The integrated DSP problem was studied both with and without considering the installation of DGs. It has been demonstrated that when DGs are considered in the DSP problem, a reduction in investment and operating costs is achieved. Furthermore, the results obtained show that the sensitivity analysis, based on the impedance matrix employed in this paper, discovers lower operating and investment costs than those reported in previous literature. This demonstrates, therefore, that the use of the impedance matrix is desirable in order to solve the allocation and sizing problems of DGs, in the DSP problem.Building on the present study, the methodology presented here should be extended further, to consider a three-phase model for both networks. Furthermore, different DG technologies, and the allocation and sizing of energy storage devices, should be considered in the DSP problem. Moreover, the variability and uncertainty of renewable DGs, demand, and energy prices should be considered in the proposed methodology. In addition, the use of a different metaheuristic technique, combined with the NSZM method, will be developed further to improve the quality of solutions and their processing time.



---

*Source: 2617125-2022-02-15.xml*