Cross‐entropy method for distribution power systems reconfiguration

Cross‐entropy (CE) is a powerful simulation method for the solution of continuous and combinatory optimization problems. The work presented here utilizes the CE method for the optimal topology of distribution power systems (DPSs). The optimal network switches are determined for the reduction of active power loss. The adapted CE method is tested on three case studies, namely, the 33‐node, 83‐node, and 880‐node DPSs. The results are compared with other reconfiguration algorithms to demonstrate the superiority of the proposed algorithm. The impact of the distributed generation is also investigated. The effective integration of the photovoltaic panels at midday, when their production is highest and meets the peak demand, is showed. Finally, the real‐time reconfiguration strategy based on the switching effort reduction is proposed and enhanced via an adequate selection of the initial switch states.


INTRODUCTION
Distribution power systems (DPSs) have gained increased importance in recent times due to the huge penetration of green energies. Because these renewable sources are characterized by their intermittent and stochastic behavior, the operating points of DPS change continuously, and to meet the optimal operation, the topology needs to follow these changes. Network reconfiguration achieved by changing the states of some maneuverable switches is widely in use by electrical engineers. It would be more practical in real time to manipulate the candidates' switches than to opt for other costly choices. 1,2 In some segments of the DPS, switches either are open or closed. 3 This option in the radial networks is used for minimization of losses as well as voltage deviations. 4 The mathematical representation of the aforementioned option includes binary variables with a nonlinear objective function and constraint. The integer nonlinear programming (INLP) methods are well suited for such cases as these tools can deal with large-scale problems. However, they must not have exponential complexity.
Usually, four families of methods are used for reconfiguration of DPS. The first concerns the method designed for a particular objective function structure, for example, integer programming, the cutting plane, 5 or branch-and-bound. 6 The second contains metaheuristic methods, which have the ability to decrease the search space and do not need derivation of the objective function. In addition, these approaches are more appropriate for INLP problems and assure the exploration of the search space. Furthermore, they have the ability to avoid the convergence to local optima of the feasible region and go to a global solution, but they are very time consuming and their solutions are not stable. In recent FIGURE 1 The topology of a distribution power system (HDE), 20 and the artificial immune systems. 21,22 The drawback of these approaches is that they are time consuming and their operation depends on careful selection of parameters. The third family is concerned with the heuristic based search family such as the method proposed by Baran and Wu 23 and the minimum cost maximum flow (MCMF) proposed by Ababei and Kavasseri. 24 These methods provide a fast, robust, and stable solution but do not guarantee the global optima.
The last family contains the determinist approach, which is the recent trend. For example, Taylor and Hover 25 have introduced three new formulations of the reconfiguration problem: the quadratic programming (QP), the quadratic constrained programming (QCP), and the second-order cone programming (SOCP). Recently, Jabr et al 26 have proposed the mixed integer convex programming (MICP) formulation and Ahmadi and Marti 27 proposed the linearized power flow with the mixed-integer quadratic constrained programming (MIQCP). Although these methods are very fast, they cannot always guarantee global optimum and are computationally time consuming. Most recently, Rubinstein 28 proposed the cross-entropy (CE) method for solving a rare-event (RE) simulation problem, and thereafter, this method is extended for the solution of optimization problems. 29 It is reported by Sebaa et al 30 and Ernst et al 31 that, in several applications, the CE method for engineering optimization problems achieves better results than other metaheuristic methods.
The CE method belongs to Monte Carlo (MC) methods. It uses the importance sampling technique (IS) in which this propriety facilitates its ability to deal with large-scale problems with poor simulation resources. MC methods are suitable for problems when other techniques are too difficult to solve or poses a heavy burden computationally. Another motivation behind the choice of the CE method is that it converges to the global optima when there are many local optima. This paper introduces the CE method to the problem of reconfiguration of the distribution networks and is organized as follow. Section 2 describes the problem statement, and Section 3 focuses on the presentation of the CE approach, especially when discrete variables are involved in the optimization problem. Finally, Section 4 presents case studies where the three known distribution networks are tested.

LOSS REDUCTION
It is beneficial to recall some concepts of the graph theory field before proceeding to the concept of DPS reconfiguration. The topology of DPS depends on the arrangement of its elements. In a graphic frame, this topology is determined by a finite set Y = {y 1 , y 2 , y 3 , … , y N } of vertices or buses and a finite set B = {b 1 , b 2 , b 3 , … , b M } of links or branches.
The Y and B sets correspond to the definition of a graph that is denoted as G(Y, B). 24 The orientation of each branch can be a priori or indefinite as shown in Figure 1.
The reconfiguration problem can be seen as the selection of the branch x i from the branches in fundamental path i containing m i branches that should be switched off. Initially, the fundamentals paths have to be evaluated.
These paths are chosen in such a way that, when the initially switched-off branch is switched on with any branch of the same path that is initially switched off, the network keeps its connectivity and the maneuver does not affect any other path. This paper assumes that the set Ω i is a priori that is known.
Noting that ind 1,i is originally the switched-off branch, then the state space is Ω = Ω 1 × Ω 2 · · · × Ω Nl and the candidate solutions have the form X = (x 1 , x 2 , … , x Nl ).

Performance function
Losses cannot be avoided and are one of the technical and economic problems, which are presented by the Joule's formula where R k and I k are the resistance and the current of the branch k, respectively. Defining the network constraints as follows: (iv) connectivity and radiality.
In this paper, the back forward sweep technique that includes two steps is used to ensure that constraints (i) and (ii) are satisfied by the load flow (LF) routine. A faster LF routine is achieved with the branch exchange approach. For such an approach, first, the LF needs to be performed on the original network represented by the vector X orig = (ind 1,1 , ind 1,2 , … , ind 1,Nl ). For example, the LF for the vector X = (ind k1,1 , ind k2,2 , … , ind kNl, Nl ) is carried out by exchanging the branches (ind k 1 ,1 and ind 1,1 ), (ind k 2 ,2 and ind 1,2 ), and so on. Furthermore, the constant load model is adopted throughout this work and as Equation (1) is composed of discrete variables, a discrete optimization solution, such as the CE method is utilized.

CE METHOD FOR DISCRETE OPTIMIZATION
Originally, the CE method is developed to estimate occurrence of an RE in stochastic networks, and then, its suitability is confirmed for several optimization problems. A sequence of probability density functions (pdf) that estimate this RE problem are generated until they converge to a degenerated pdf and the global optimum is reached.
Two steps could summarize the CE method.
-Generation of random variables, vectors, or paths according to specific probability density functions (pdfs).
-Updating the pdfs based on the previous step best samples.

CE for optimization
Let G(·) be the objective function to be minimized in the domain Ω. Suppose that the task is to find out the global minimum of G(·) over Ω and the global optima is X * and G(X * ) ≤ * G (X * ) = * = min X∈Ω G (X) .
As the CE method for optimization is an extension of the CE method for RE estimation, the event associated with (1) is the probability that is represented by {G(X) ≤ }. Considering that X belongs to some pdf f (x, u) on Ω domain, hence X depends on u and the level . Thereby, the randomness aspect is introduced to the original problem to provide a stochastic model. When is close to * , the probability λ = Pr(G(X) ≤ ) become very small and it is considered as an RE probability. The IS technique is one of the methods that could estimate this probabilitŷ where X k is an independent and identically distributed (iid) samples generated from a known pdf g. The pdf that gets a zero variance estimator is g * (x) = f (x)I {G(x) ≤ } / . However, it depends on the unknown probability . The role of the CE method is to determine an IS pdf f (·, v) that is close to g * in terms of the CE distance (or Kullback-Leibler distance) The parameter v should maximize the following expectation (u is the parameter of the previous pdf): One more time, using other IS with pdf f (·, w), the next expectation should be maximized based on v is the likelihood ratio between the two pdfs. The parameter v can be computed based on the called stochastic counterpart (SC) problem.
Here, X k are iid samples generated from f (·, w). It is easy to solve the SC problem if the pdf f (·,v) belongs to a special distribution family like exponential, normal, or Bernoulli lows. This will be done by the computation of gradient with The typical CE algorithm for optimization is as follows.
1. Give the input data a. the initial parameter vector v 0 , the pdf family f (·, v), and the domain Ω b. the sample size L and the quantile ratio and the degeneration error Set t:=1 the iteration counter and Compute the objective functions G (X k ) and sort them in increasing order iii. Set t ≔ G(X elite ) iv. Solve the SC problem and get v t

CE for the reconfiguration problem
The reconfiguration of DPS for loss reduction is a discrete optimization problem, 29 where the state space is After solution of the SC problem (7), it is easy to prove that Considering ∑m k=1 p k, = 1. The updated probabilities are where P 0 and P t are the vectors of initial and previous probabilities. As this optimization involves discrete random variables, the likelihood ratio is set to W = 1. Therefore, (9) becomes where the Elite is the set of the variables X k that verifies G (X k ) ≤ and N Elite = #Elite. The CE algorithm for the reconfiguration of DPS is as follows.
The pdf family f (·, P t ) (8) e. The objective function for the original configuration G (X orig ) f. The sample size L, the quantile ratio and the degeneration error 2. Set t:=1 the iteration counter and N Elite = L 3. While the ||P t −P t−1 || ≥ i. Generate idd random vectors X k (k = 1:L) belonging to f (·,v t−1 ) ii. Compute the objective functions G (X k ) and sort them in increasing order iii.

End of the while loop
At each iteration, the probability p t+1 k, will be updated. At the iteration t, the solutions are sorted according to their objective function, then only the best N Elite samples are kept. The probability p t+1 k, is simply the occurrence of the branches ind k,j in the loop j and of course among elements of Elite .

33-node DPS
This test case shown in Figure 2 is known as Baran and Wu's network. 23  The original or the total circuit loss is 202.69 kW. For the application of the CE method, the program is written and run with MATLAB ® that operates under Core i5, 2.5 GHz and 8 Gb RAM.
Usually, the quantile ratio is selected between 0.01 and 0.1. It has been suggested that the size of elite set should be sufficiently large to provide a good estimation of pdfs (9). In this paper, ε = 0.1 and L = 50 are selected for the desired FIGURE 2 A 33-node distribution power system with five fundamental loops (in color)

FIGURE 3
Evolution of the probabilities using the proposed cross-entropy method 32 Figure 3, the CE method converges to the global solution within 3 iterations, the optimal solution is to switch off the lines {32, 37, 7, 9, 14}, which leads to the reduction of network losses to 139.55 kW.
For the purpose of comparison, Table 1 presents the results obtained by approaches proposed in literature and the CE method. It shows the optimal configuration, the power loss, and the convergence time for the Baran's network. The convergence time is achieved with 5 iterations, which can be multiplied by 3/5 as the optimal performance is reached by 3 iterations. This, in turn, reduces the proposed method's convergence time to 1.1713 seconds. All the presented approaches   have met the optimal configuration, except the heuristics approaches of Baran and Wu 23 and Ababei and Kavasseri 24 that differ from the global optima by 3.57%, but their convergence time is short.

94-node DPS
It is also known as the Taiwan DPS and has been introduced in the work of Su et al. 34 It contains 94 nodes, 11 feeders, 83 normally closed switches, and 13 open lines. Therefore, the number of fundamental loops is 13. Its nominal voltage is 11.4 kV as depicted in Figure 4. The quantile ratio and the sample size are 0.1 and 180, respectively. The CE approach converges within 13.476 seconds to the global solution, which is clearly fast than the convergence time of other methods. However, the Baran and Wu's method is the fastest and achieves the global optima. Results obtained by approaches proposed in literature and the CE methods are shown in Table 2.

880-node DPS
This test case is created using data from 135-node DPS and data from 201-node DPS. It contains 880 nodes, 7 feeders, 873 normally switched-off lines, and 27 tie lines. 24 Figure 5 shows the simplified representation of this system. Abbreviations: ACSA, ant colony search algorithm; CE, cross-entropy; DPS, distribution power system; HDE, hybrid differential evolution; MCMF, minimum cost maximum flow; MICP, mixed integer convex programming.  Table 3 summarizes the results obtained by approaches proposed in the literature and the CE method. The best result obtained by the proposed CE algorithm in 10 runs when L = 600 and = 0.1 are selected. However, application of the metaheuristic algorithms for this case, which is a large-scale system has not yet been reported. Inspection of Table 3 suggests that methods of Baran and Wu 23 and Ababei and Kavasseri 24 converge to a near optimal solution with a difference of 1.19% and 1.33%, respectively, from the global solution, whereas the MILP 26 method gives a solution with a difference of around 0.19% from the optima solution with a converge time of 398 seconds.
The QP 25 method with relaxation (in the relaxation the binary variables are considered as continuous by adding constraints, like y 2 −y = 0) converges to a near optimal solution with a difference of 0.003% from reaching to the global optima with a converge time of 10 hours and 22 minutes. The proposed approach reaches the global optima in 214.269 seconds, which, as yet, has not been achieved by other researchers.
Based on the analysis of results obtained for various networks, one could confidently state that the CE method is suitable for network reconfiguration and it converges to the global optima within reasonable computation time. In addition, the proposed method by nature is suitable for parallel computation techniques, which, in turn, aids the reduction of the computation time in the case of large-scale networks.

Convergence analysis
The CE results presented in Tables 1 to 3 are not stable as the sample size is relatively small. To obtain a robust and stable solution, the sample size should be at least equal to five times the problem size #P. In Table 4, for every problem size, the  algorithm has run 100 times. The probability to reach the global optima G * is equal to ∑ 100 i=1 I{G i −G * }/100, where G i is the solution found at the ith run and the suboptimality of the problem is ∑ 100 i=1 [(G i −G * )/G * ]/100. Over 100 runs, the proposed method has not been able to meet the global optima twice for the Baran and Taiwan test cases and once for the Ababei test case. The CE algorithm converges every time to the degenerated pdfs within a few iterations. Moreover, the convergence time increases quadratically with the problem size as depicted in Figure 6. However, as the sample size increases linearly with the problem size, the complexity of the algorithm is defined by O(n 2 ).
The outcome of the proposed approach is similar to that of those presented in the reference papers (shown in Tables 1 to  4). In addition, the proposed method does not require enormous simulation materials, as by its nature, the CE approach uses the elitism concept and at each step the new pdf is adjusted by the previous elite set. Hence, the reconfiguration, based on the CE method, is relatively quicker than other methods (the reference papers).

Limits on the maneuvers number
In order to include limits in the maximum number of switches, the following constraint is considered: where n max , n sw , and x i are the maximum number of switches permitted, the current number of switches, and the index of the switch in the line of loop i, respectively. This constraint is considered when generating samples by the acceptance-rejection method or the Gibbs sampling via the truncation of the samples do not follow Equation (11). However, the easy way to handle this constraint is to use the penalty function d. The new performance function becomes where the penalty function is and P orig loss is the total active power loss of the original configuration (network before reconfiguration process). Figure 7 shows the influence of the maximum number of switches on the total power loss of the 94-node DPS. This DPS contains 13 possible switches. Hence, the number of n max takes values from 0, 1, up to 13. For n max ≥ 9, the total active loss is constant because the optimal solution corresponds to the switches of only nine lines. A good performance can be achieved by the switches of only five lines, which corresponds to the reduction of 11.53% of the total power loss. Note that the optimal solution achieves a reduction of 11.68%.

Distributed generators
In this section, a 94-node network is considered to study the impact of the photovoltaic panels (PV) that are installed at the generation node, on the network global loss. These PVs are assumed to have the same capacity as the traditional generators connected to the same node. Figure 8 shows the profile of PVs power. Furthermore, the loads are varying in time and are supposed to have a household's profile. Other profiles such as the industrial profile or a mixture of these two profiles could also be investigated. Only the load profile at the 15th node is presented in Figure 9.
The objective function must include various criterions of the classical objectives such as voltage deviation improvement, load balancing, and the global active loss. In addition, attempts should be made to reduce the number of switching maneuvers during the day. To satisfy all of the above objectives, a compound of these criterions would guarantee the aim min ⟨ P loss P orig loss where x m and x old m respectively are the current and the previous indices of the switched-off branch in the loop m (I), which is known as the indicator function.
The optimization problem described by (14) should be repeated every 15 minutes (96 samples every 24 hours). Based on how the initial states of switches are selected, two cases could be analyzed. The first case considers the optimal configuration of the network as discussed earlier and the second case considers the original configuration of the network. From the results presented in Figures 10 to 12, one can clearly observe that the first choice leads to a minimum number of maneuvering of switches in 1 day, whereas the second case shows the importance of switching. Furthermore, Figure 13 shows that the voltage drop for both cases is low (around 5%) and similarity exists between their voltage profiles. In addition, Figures 14 and 15 depict similarity in power loss and current balancing profiles.    Figure 16 confirms that the supplied non-green power reaches its minimum value at midday, ie, when the PV production is the highest. Moreover, the blue curves in Figures 13 to 16 represent the scenario without the green sources. These Figures show that, at midday, the impact of the PVs is important. However, the losses are slightly lower because, in the case "with DG," the importance is given to the consumption of green energy than to the loss minimization.
The same analysis is conducted for an 880-node power system assuming that this network is subject to the same amount of PVs penetration at nodes 1, 3, 4, 5 with a maximum power of 44.32 MW at 14 hours (02:00 pm). The results are presented in Figures 17 to 23 and confirm the advantage of the green energy. The average number of switches (the maximum number of switches is 457) that change state or initialize is 11.20. However, in general, the average number of switches without the initialization is 11.51.

CONCLUSION
The proposed CE algorithm is successfully applied to the problem of reconfiguration in DPSs. The results demonstrated that the optimal state of line switches can be obtained within a few iterations with reasonable convergence time and the computational complexity of this algorithm displays quadratic O(n 2 ) characteristics.
Furthermore, based on analysis of the results, the superiority of the proposed method is confirmed over its counterpart heuristic and metaheuristic algorithms. It can even compete with direct mathematical methods, such as QP and MILP. These findings are further reinforced when the proposed algorithm is applied to a DPS containing many distributed generators, such as PVs that provide power over a half day (daylight) or shorter period, and in addition, when the load and the generation are varied.

CONFLICT OF INTEREST
The authors declare that there is no conflict of interest regarding the publication of this article.

NOMENCLATURE Y
A finite set of N buses or nodes B A finite set of M branches G(Y, B) Axiomatic representation of a graph Nl Number of the fundamental paths Ω i The set of the branch indices belonging to the ith path m i The number of branches in the ith path ind j,i The jth branch belonging to the ith path Ω The state space X = (x 1 , x 2 , … , x Nl ) ∈Ω The candidate solution G(X) A performance/objective function R k and I k Resistance and the current of the branch k X* The global optima * The minimum of G(·) f (·, u), g(·, v) pdfs with a parameters u or v following a special family of pdf I {·} The indicator function is equal to 1 if the logic expression {·} is true otherwise it is 0 D(·, ·) The Kullback-Leibler (cross-entropy) distance E f (·) The expectation of the function f (.) W (·, u, w) Likelihood ratio between the two pdfs f (·, u) and f (·, w) L The sample size The quantile ratio Degeneration error N Elite The size of the elite set Elite t The iteration counter p t k, the probability of the occurrence of the ind k,j at the iteration t ORCID Karim Sebaa https://orcid.org/0000-0002-8896-2823