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

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 intensiﬁcation and diversiﬁcation 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 eﬃciently and eﬀectively 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 eﬃciency 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.


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.
ese 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][5][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].e number of papers regarding DSP for a medium voltage [10][11][12][13][14][15] is higher than for low voltage [16][17][18][19].
Moreover, only a few papers consider both levels [20][21][22][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.
e 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.e 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.e 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.
e 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.e 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].
is 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.e 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. is 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].is 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).
e 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.e 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]. is test is a pairwise comparison, used in hypothesis-testing situations.e Wilcoxon signed ranks test is chosen because it is less influenced by sample size (number of case problems) than other nonparametric tests [33].
e 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.e 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 2 Journal of Electrical and Computer Engineering 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.is 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.
e 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.e difference from other papers within the literature on the subject is that all the installing decision variables are represented in a single vector.is 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.e 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.e 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. is 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.
e proposed methodologies are tested on a distribution system, with the results being presented in Section 4. e results are analyzed and discussed in that section.Lastly, conclusions and future work are discussed in the final section.

Mathematical Formulation of the Problem
e impedance matrix is used in this mathematical model to represent both networks as a single system [34].e 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. ( 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.

Journal of Electrical and Computer Engineering
Equation (11) determines the injected power in each DT.Equation (12) determines the power generated from the installed DG.
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.

Methodology Proposed
e MINLP model presented here is solved using two wellknown, single solution-based metaheuristics (simulated annealing and iterative local search), combined with a new neighborhood search method based on the impedance matrix (NSZM).e 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.
us, the NSZM method enhances the search process of both metaheuristics.
e 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].e 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.e transition between solutions is controlled by each metaheuristic.
e main aspects of the new NSZM method and the SA and ILS algorithms used in this paper will be discussed in the following.

Codification.
e integrated DSP problem considering DG is represented using the new codification vector of Figure 1. is 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.
e vector used can be partitioned into five areas.e first area contains the locations and sizes of the existing and new primary feeders (size n 1 + n 2 ).e second area includes the locations and sizes of the existing and new secondary circuits (size n 3 + n 4 ).e third area relates to the locations and capacities of the DTs (size n 5 ).e fourth area relates to the locations and sizes of the existing and new substations (size n 6 + n 7 ).
e fifth area relates to the locations and capacities of the DGs (size n 8 ).

Neighborhood Search Method Based on the Zbus Matrix (NSZM).
is method proposes a defined neighborhood structure and uses the Zbus matrix to obtain electrical and physical information about the solutions network.e information obtained from the Zbus matrix allows one to explore the solution space more efficiently and effectively.
e 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.e symbols °and ⊘ are the Hadamard product and division, respectively.
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 Z BUS n×n � R BUS n×n + jX BUS n×n and I BUS n×1 � I BUS n×1 ∠Θ 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): Journal of Electrical and Computer Engineering 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 →BUS * ,k 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.
e 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. is information is obtained from equations ( 26)- (29).
Equation ( 26) is used to find out if a specific node is a terminal node.e node k is a terminal node, if equation (26) is true for all i � 1, 2, . . ., n.
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.
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.
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.
e physical information obtained from equations ( 26)-( 29) is used in the neighborhood structure to make the changes proposed by the criteria.ese 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.

Neighborhood Structure.
e 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.e neighborhood structure considered is explained in the following: e branch with max ΔU ij   is selected to connect the link.Next, the branch with min I ij   of the loop formed is selected to be removed.e 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.e loads are transferred from the nearest existing DT to the new one, by removing the circuit with min I ij   of the loop that is formed.ese currents are calculated using equation (23).Finally, the operative limits of the DTs involved in the process are checked and corrected.(iii) DT Removal.
e DT to be removed is chosen randomly from the three DTs with the smallest installed capacities.e loads of this DT need to be transferred to another existing DT. e 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

Primary feeders
Secondary circuits DTs Substations DGs 3 6 Journal of Electrical and Computer Engineering than one DT that fulfills these conditions, the DT is chosen randomly from among them.Loads are transferred by removing the circuit with min I ij   of the loop formed.ese 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.
e neighborhood structure uses equations ( 26)-( 29) through the process to determine the physical network information of the solutions.

Modifications of the Zbus Matrix to Reflect Changes Made by the Neighborhood Structure.
e Zbus matrix can be modified to reflect changes that are made to a network without the need to rebuild it.ese 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 Figure 2).ese changes, produced by the neighborhood structure, are reflected in the Zbus matrix by the addition of links in the network (see Figure 2).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.e 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.is procedure is based on the elimination of a node in a matrix where the independent variable of this node is equal to zero.e 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 nodes p − 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 nodes B − 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)).
e procedure for upgrading or degrading the capacity of an element is explained next.Between k − 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)).e 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.

Evaluation of the Configurations.
e 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].
e 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.
e combination of objective function plus penalty costs is called the "fitness function" (F fit ) and has been widely used in [9,10,21,23,30].Penalty costs are expressed in monetary units.In this document, the SA and ILS algorithms are used and explained in the following.

Simulated Annealing (SA) Algorithm.
e SA algorithm is a stochastic optimization procedure that models the process of annealing in solids [35]. is technique simulates the heating and controlled cooling of a material to increase the size of its crystals and reduce their defects.e process ends when thermal equilibrium is achieved.e 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 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. exp ermal equilibrium can be achieved by defining an adequate cooling scheme.e computational time, and the quality of the final solution, will depend on the selection of the cooling scheme and the tuning of its parameters.e cooling scheme employed is presented in equations ( 32) and (33).Equation (32) calculates the initial temperature T 0 , equation (33) calculates the new value for temperature T k+1 at level k, and equation (34) calculates the new value for the number of transitions N k+1 at level k.For more detailed information, the cooling scheme used is explained in [36].(34) e algorithm stops when the incumbent solution is not improved for a defined number of temperature levels k, or when a fixed number of temperature levels k for the cooling process is achieved [32,36].
e search process for the SA algorithm is repeated until a certain number of transitions N k is reached.Next, the temperature T k 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.e cooling process ends when the stop criterion is met.

Iterated Local Search (ILS) Algorithm
. ILS is a metaheuristic, based on a simple idea [37].e 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.is 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.
e 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.e 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.is process is repeated until a certain number of global iterations is completed.

General Methodology of the SA-NSZM.
e 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 Subsection 3.6 e NSZM method is then used to obtain a new solution for the SA algorithm, as explained in Subsection 3.2.e 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.
e search process is repeated until N k transitions are completed.For the next temperature level, T k and N k 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.e cooling process ends when the stop criterion is achieved.Figure 3 shows the proposed SA-NSZM methodology.

General Methodology of the ILS-NSZM.
e 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 (n d , i local , and j global ) are initialized.Following this, the NSZM method is used to obtain a new solution for the ILS algorithm, as explained in Subsection.e new configuration represented by the Zbus matrix is evaluated, using the load flow and the fitness function. If , the new solution is accepted.Otherwise, the new solution is rejected.
is local search process is repeated until i local 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 n d times.ese perturbations are repeated until j global global iterations are finished.
e ILS-NSZM algorithm ends when the global search process has been repeated to the point at which j global global iterations are finished.Figure 4 shows the proposed ILS-NSZM methodology.

Description of the System.
e SA-NSZM and ILS-NSZM methodologies were tested on the real, medium-low voltage distribution system presented in [23].e 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.
e full system database is presented in [38].

Parameters of the SA-NSZM and ILS-NSZM Algorithms.
e parameters of both algorithms were adjusted by a trial and error process [35,37].e factors used for the cooling scheme of the SA-NSZM algorithm were η � 1%, λ � 5%, α � 1.05, and β � 0.95.e cost of the initial solution F(x 0 ) for Cases A and B was 1.85 million USD, and for Cases C and D it was 2.03 million USD. e initial number of transitions N 0 for Cases A and C was 293, and for Cases B and D it was 308.N 0 was defined as the number of variables of each problem.e initial temperature T 0 , according to equation (32), becomes 6175.6 for Cases A and B and 6798.3 for Cases C and D. e stopping criteria used were 15 temperature levels if the incumbent solution is not improved and 60 temperature levels for the cooling process.According to Journal of Electrical and Computer Engineering equation (33), the final temperature T f for 60 levels becomes 284.50 for Cases A and B and 313.20 for Cases C and D.
e factors used for the ILS-NSZM algorithm were n d � 3, i local � 100, and j global � 100.

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. e 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.e   1 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 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]. is test is a pairwise comparison, employed in hypothesistesting situations.e 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].e Wilcoxon signed ranks test uses the average of the best objective function values of Table 1.e Wilcoxon signed ranks test is chosen as it is less influenced by sample size (number of case problems) than other nonparametric tests [33].
e Wilcoxon test uses the difference between the performance scores of the two algorithms on i th out of n problems.
e 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.e 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.
is 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].
Table 2 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).e results obtained by way of the Wilcoxon test show that the null hypothesis is not rejected.e differences observed between the algorithms are therefore explicable by chance alone.is means that the NSZM method enhances both algorithms equally and can be applied to enhance any single, solutionbased metaheuristic in order to solve problems related to the expansion planning of power systems, considering the installation of DGs.
A more detailed comparison of the best results obtained is shown in Tables 3 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] (2) (3) (2) ( 12 Journal of Electrical and Computer Engineering 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. ( (2) (2)

Journal of Electrical and Computer Engineering 13
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.us, 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 shows that the SA-NSZM and ILS-NSZM algorithms discovered lower operating costs, using less (3) ( (5) (2) (2) ( (2) ( 14 Journal of Electrical and Computer Engineering 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.e development of efficient and effective methodologies that can solve the DSP problem, considering DG integration, is therefore highly desirable and necessary.
In Table 3, it can be seen that the use of the NSZM method proposed in this paper allows the SA and ILS (2) (2) ( Journal of Electrical and Computer Engineering 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, Table 3 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.e 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.
Tables 3 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].
e LV buses with DGs installed for Cases B and D are presented in Table 5.One can see that the locations and sizes obtained for DGs differ between cases.

Conclusions
is 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.
e 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.e proposed methodologies have been validated and tested using a real MV/LV distribution system, taken from the specialized literature.e 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.
e results obtained with the test show that the differences observed between the two algorithms are explicable by chance alone.is 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.
e 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.e 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.e 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.
is 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 threephase 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.

Figure 1 :
Figure 1: Codification scheme for the integrated DSP problem considering DG.
Metaheuristics.A single solutionbased 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 Journal of Electrical and Computer Engineering neighborhood search, variable neighborhood descent method, guided local search, iterated local search, and their variants.
e stochastic mechanism accepts a new solution depending on the values of the temperature level T k and the objective function for the current and new solutions (F current fit and F new fit ).If F new fit ≤ F current fit the new solution is accepted, and it replaces the old one.If F new fit ≥ F current fit

Table 1 :
Objective functions in millions of USD.

Table 2 :
Wilcoxon signed ranks test results.

Table 3 :
[23]arison with results reported in[23]in millions of USD, Cases A and B.

Table 4 :
Comparison of best solutions in millions of USD-cases C and D.

Table 5 :
Locations and types of DG installed for Cases B and D. : Installation cost of a new DT at bus i, type d ($) C N DG i,g : Installation cost of a new DG at bus i, type g ($) Percentage of worsening for the objective function of the initial solution (% /100) F(x 0 ): Objective function of the initial solution β: Rate of change for the temperature α: Rate of change for the number of transitions Z P ij,p : Impedance of a MV wire at buses i − j, type p [Ω] Z S ij,c : Impedance of a LV wire at buses i − j, type c [Ω] : Decision variable to install a new DT at bus i, type d Power generated by a substation at bus i at load level l (MVA) U BUS i,l : Bus voltage at bus i at load level l (kV) I BUS i,l : Bus current at bus i at load level l (A) ZIP P i,l : ZIP load model at MV bus i at load level l ZIP S i,l : ZIP load model at LV bus i at load level l Z BUS ij : Zbus matrix [Ω] Ω EP : Set of existing MV branches Ω ES : Set of existing LV branches Ω ESS : Set of existing substations Ω N DT : Set of new DT Ω N DG : Set of new DG Ω NL : Set of the number of levels of the load demand curve Ω NP : Set of new MV branches Ω NS : Set of new LV branches Ω NSS : Set of new substations Ω PF : Set of new and existing MV branches Ω PN : Set of buses of MV network Ω SC : Set of new and existing LV branches Ω SN : Set of buses of LV network Ω SS : Set of new and existing substations Ω T DT : Set of types of DT Ω T DG : Set of types of DG Ω TP : Set of types of MV wires Ω TS : Set of types of LV wires Ω TSS : Set of types of substations.
i : Lower voltage limit at bus i (kV).s : Decision variable to install a new substation at bus i, type s.I P ij,l : Current of MV branch i − j at load level l (A) I S ij,l : Current of LV branch i − j at load level l (A) S DT i,l : Power injected by a DT at bus i at load level l (kVA)