Electric Power Grids Distribution Generation System for Optimal Location and Sizing—A Case Study Investigation by Various Optimization Algorithms

Abstract

In this paper, the approach focused on the variables involved in assessing the quality of a distributed generation system are reviewed in detail, for its investigation and research contribution. The aim to minimize the electric power losses (unused power consumption) and optimize the voltage profile for the power system under investigation. To provide this assessment, several experiments have been made to the IEEE 34-bus test case and various actual test cases with the respect of multiple Distribution Generation DG units. The possibility and effectiveness of the proposed algorithm for optimal placement and sizing of DG in distribution systems have been verified. Finally, four algorithms were trailed: simulated annealing (SA), hybrid genetic algorithm (HGA), genetic algorithm (GA), and variable neighbourhood search. The HGA algorithm was found to produce the best solution at a cost of a longer processing time.

1. Introduction

In the past decade, increased distributed generation (DG) has led to profound changes in electricity distribution networks. Several factors have driven DG (defined as production of electricity close to consumption centers) including new technological advances in the production of electricity on a small scale, a preference for the use of renewable resources, difficulties in network expansion, and a growing interest in incorporating demand and active agents in the electricity markets [1]. DG can contribute to reducing losses, improving voltage profile, improving reliability, and postponing investments [2]. However, as noted in [3], harnessing the benefits of DG depends largely on its location, sizing, and network features. That is why the past decade's alternative methodologies for proper location and sizing have been explored [4].

In [5] a literature review of techniques used for the location and optimal sizing of DG in distribution networks is presented. The authors classify the techniques according to analytical methods, metaheuristics, and mathematical programming [6]. It should be noted that for the problem under study metaheuristics techniques have significant advantages over classical mathematical programming because of the nonlinear and non-convex relationships in the location and sizing of DG. On the other hand, the main weakness of metaheuristics is that they do not guarantee obtaining a global optimum. However, metaheuristics may provide a solution, or set of high quality solutions. Another advantage that lies with metaheuristics techniques is that they allow the use of more detailed models of network operation than the analytical model. To apply mathematical programming techniques to the problem of optimal location and sizing of DG, it is necessary to use linearization or approximations to the equations in the balance of power [7].

Metaheuristic optimization techniques have been widely used in the location and design of the DG. These techniques include genetic algorithms, tabu search, and colonies of particles. In [8], a model of multi-objective optimization was presented to determine the location and optimal sizing of DG using the technique of simulated annealing. The elements were modelled power losses, the number of generators, the voltage profile and power injected by the DG. The method sought a solution that would improve the voltage profile and reduce losses to the minimum DG units in the system. In [9], the authors presented a population-based metaheuristic based on the parasitic reproductive behavior of some species of cuckoos (cuckoo search algorithm) for the location and optimal sizing of DG. The advantage of this algorithm was that it required a few parameters to calibrate. The aim of the study was the reduction of active power losses.

Hybrid methods combine two or more search techniques in order to exploit their potential and compensate for their deficiencies. The most common methods combine population hybrid techniques with methods that enhance some kind of local search or alternatively, heuristic methods with classical mathematical programming. In [10], a combination of simulated annealing and genetic algorithms for optimal location of DG comes with network distribution. The objective was to minimize losses. This shows that the combination of a genetic algorithm (GA) with simulated annealing was more effective than using only a GA. In [11], a method to maximize the benefit to network operators and owners of distributed generation in a deregulated electricity market hybrid algorithm was presented. As well as simultaneously optimizing the benefits to the distribution company, and to the owner of the DG, their method also considered the uncertainty of demand and energy prices [12].

In [13], a hybrid algorithm was presented that improved stress profiles and reduced emissions using a particles’ gravitational colony search to determine the proper location and sizing of the DG that minimized loss. The aim of this article is to contribute to the discussion on the effectiveness of heuristic and metaheuristic methods for optimal dimensioning and location of DG. Four different techniques were implemented and compared, namely (i) simulated annealing; (ii) variable search environment; (iii) genetic algorithm; and (iv) a hybrid method that combines variable search environment with a genetic algorithm. To test the efficiency of these methods, they were applied to various tests in a distribution system (34 bars) that is widely used in the technical literature.

2. Background Mathematical Formulation

The objective function of the proposed problem is to improve the voltage profile and reduce system losses. To this end, the indexes defined in [14] described below are taken. The rate of voltage profile, denoted as IPT, is defined by Equation (1). This index takes into account the voltage, bars, and load expressed as Power of the system as a load factor. This makes it more important to maintain proper high-voltage in bars under higher demand.

$$IPT={\displaystyle \sum }_{\left\{i=1\right\}}^N=V_{i{L}_i}$$

(1)

Variables: $V_i$: High-voltage in the bar i (p.u); $L_i$: Load in bar i (p.u); and N: Total number of bars. The rate of improvement of voltage profile, denoted as IMPT is given by the Equation (2).

$$IMPT=\frac{IP{T}_{\left\{wDG\right\}}-IP{T}_{\left\{woDG\right\}}}{IP{T}_{\left\{woDG\right\}}} \times 100$$

(2)

Variables: $IP{T}_{\left\{wDG\right\}}$ is the index of the system voltage profile with DG (p.u.) and $IP{T}_{\left\{woDG\right\}}$ is the profile index stress the system without DG (p.u.). Note that the IMPT denotes the percentage improvement in IPT with DG. The second objective is to reduce active losses and compare losses with the system with and without DG, given by Equation (3).

$$IP{L}_{\left\{wDG\right\}}={\displaystyle \sum }_{\left\{K=1\right\}}^N{I}_{\left\{K,wDG\right\} }^2{R}_k\times D_k+ L{L}_{\left\{wDG\right\} }$$

(3)

A similar expression can reflect woDG conditions. The equation variables are described as follows: $IP{T}_{\left\{wDG\right\}}:$Index line losses with DG; $R_k$: Resistance in line K (p.u/km); $D_k$: Length in line K (km); $I_{\left\{\mathrm{K},wDG\right\}}^2:$ Current line K with DG (p.u); $L{L}_{\left\{wDG\right\}}$: Index line losses with DG; $L{L}_{\left\{woDG\right\}}$: Index line losses without DG; $I_{\left\{K,woDG\right\}}^2:$ Online K stream without DG (p.u).

While the percentage reduction in losses in the line is income tax given by Equation (4):

$$IRPL=\frac{IP{L}_{\left\{wDG\right\}}-IP{L}_{\left\{woDG\right\}}}{IP{L}_{\left\{woDG\right\}}} \times 100$$

(4)

The objective function is to minimize the reduction rates of losses and to improve the voltage profile. In this case, you should be using weighting factors $W_1$ and $W_2$ and dimensions for each rate in order to assess their importance in the optimization process. To be noted, maximizing the voltage profile can be achieved by Equation (2) and set by the objective function from an Equation (1). Minimizing the voltage profile can be readily available from Equation (4) and it can be by the objective function from an Equation (3). Both the voltage profile cases the objective function can be determined from Equation (1) to Equation (3).

The optimization problem to solve is described by expressions form Equations (5) to (17).

$$f\left(x\right)=W_{1IRPL}+W_{2IMPT}$$

(5)

where the subject of the previous equation is:

$$0\le W_m\le 1; m=1,2$$

(6)

$${\displaystyle \sum }_{\left\{m=1\right\}}^2{W}_m =1$$

(7)

$$u_i{P}_{Gi}-P_{Di}- V_i {\displaystyle \sum }_{k=1}^{nb}\left[V_k\left(g_{ik}\cos {\theta }_{ik}+b_{ik}\sin \left({\theta }_{ik}\right)\right)\right]=0$$

(8)

$$u_i{Q}_{Gi}-Q_{Di}- V_i {\displaystyle \sum }_{k=1}^{nb}\left[V_k\left(g_{ik}\sin {\theta }_{ik}+b_{ik}\cos \left({\theta }_{ik}\right)\right)\right]=0$$

(9)

$$P_{ik}=V_{i{g}_{ik} }^2- V_i{V}_{k{g}_{ik}}\cos \left({\theta }_{ik}\right)- V_i{V}_{k{B}_{ik}}\sin \left({\theta }_{ik}\right)$$

(10)

$$Q_{ik}=V_{i{b}_{ik} }^2- V_i{V}_{k{g}_{ik}}\cos \left({\theta }_{ik}\right)- V_i{V}_{k{B}_{ik}}\sin \left({\theta }_{ik}\right)$$

(11)

$$S_{ik}^2=P_{ik}^2+Q_{ik}^{2 }$$

(12)

$$P_{Gi}^{\min }\le V_i\le P_{Gj}^{\max }$$

(13)

$$V_i^{\min }\le V_i\le V_i^{\max }$$

(14)

$$S_{ik}\le S_{ik}^{\max }$$

(15)

$$N_{DG}\le N_{DG}^{\max }$$

(16)

$$u_i\in \left\{0,1\right\}$$

(17)

where $P_{\left\{Gi\right\}}$ and $Q_{\left\{Gi\right\}}$ are the active and reactive power respectively, delivered by a unit of DG if it is located in the bar i. Note that not all bars have DG. DG for each unit must be assigned to a binary variable (called $u_i$). For simplicity, it is not considered that the DG inject or take reactive power from the network, so $Q_{\left\{gi\right\}}= 0$. $P_{\left\{DI\right\}}$ and $Q_{\left\{DI\right\}}$ correspond to power demands active and reactive bar i, respectively. Also nb is the number of bars, ${\theta }_{\left\{\mathrm{ik}\right\}}$ is the angular opening between the bars and k; $g_{\left\{\mathrm{ik}\right\}},$ and $B_{\left\{\mathrm{ik}\right\}}$ are the real and imaginary, respectively, of the nodal admittance matrix parts. Constraints (8) and (9) represent the balance equations and reactive power, respectively. Restrictions (10)–(12) represent the equations of active power flow, reactive and apparent power, respectively. The restrictions (13)–(15) consider the power limits injected by the DG, limits voltage network, and load flow limits, respectively. The constraint (16) indicates the maximum number of units DG needs to consider and restriction (17) indicates the binary nature of the variables $u_i$ (1 if with DG and 0 if without DG). The model described by Equations (5)–(17) corresponds to a problem in mixed integer nonlinear programming which is highly dimensional and non-convex having multiple local optima, which justified its solution using the search methods illustrated in this article.

3. Methodology on Hybrid Optimization Algorithm

To address the problem of optimal location and the sizing of DG described in the previous section four techniques were used as combinatorial optimization: Simulated Annealing, Variable Descending Search Environment, Genetic Algorithm and Hybrid Genetic Algorithm. A brief description of each technical solution as adopted in this study is presented [15].

3.1. Simulated Annealing

Simulated Annealing (SA) emulates the annealing process in steel and ceramics, which involves heating and then slowly cooling the material to vary its physical properties. This procedure was introduced in [16]. In each iteration of the SA some neighbours of the current status are evaluated and probabilistic decision made between making the transition to a new state or remaining in the current state. If the neighbour solution enhances the value of the objective function is accepted with probability 1, otherwise the probability of accepting by the Metropolis criterion given by Equation (18) where the parameter c corresponds to the temperature.

$$Prob\left(accept x^{\prime } \right)=\{\begin{array}{cc}1 & , f\left(x^{\prime } \right)<f\left(x\right)\\ \exp \left(-\frac{f\left(x^{\prime } \right)-f\left(x\right)}{c}\right)& , f\left(x^{\prime } \right)\ge f\left(x\right)\end{array}$$

(18)

SA assesses unattractive solutions in the early stages, then as the temperature parameter is reduced, the search becomes more selective, lessening the declines in the objective function.

3.2. Variable Search Environment Descendin

Environment Variable Search (EVS) is a metaheuristic based on a local search in a changing neighbourhood (also known as environment structure) [17]. Variations in EVS are given the names down, reduced, basic or general EVS. This paper considered an extension of EVS known as Variable Search Environment Descending (VSED) in which the current solution obtained from the change in a local search is implemented; as long as this one has found a better solution. VSED is illustrated below [4] as the following:

• Initialization: Select the set of environments, structures $N_k, k = 1,\dots , k_{\max }$ to be used in the descent. Find an initial solution x;
• Iterations: Repeat until no improvement is obtained (until there is no more optimization that we can get).
  In the following sequence:

  (1)

  Make $k\to 1$.

  (2)

  Repeat until $k = k_{max}$ the following:

  (a)

  Exploration of the environment: Find the best solution $x^{\prime }$ of the kth neighborhood of $x(x^{\prime } \in N_k\left(x\right))$

  (b)

  Move or not: If the obtained solution $x^{\prime }$ is better than x, do $x\to x^{\prime },k\to 1$; otherwise do $\to k+1$.

In the study presented in this article, environments or neighbourhoods were defined as the size (increase or decrease the capacity of the DG) and location (DG move to a neighbouring node).

3.3. Genetic Algorithm

Genetic Algorithms (GA) solve optimization problems by simulating natural selection or “survival of the fittest” [18]. The general routine of a GA is to generate an initial population of random or pseudo-random solutions. Each individual in the population is defined by a string of bits. In this case, the objective function is to evaluate the power flow as a consequence of the location and sizing of the DG [19]. For a given tournament a number of individuals are selected. The number of tournaments is equal to the size of the population. Recombination is made at one randomly selected point. The mutation is created by changing a bit (zero to one) randomly with a given probability of occurrence. Individuals generated in the process of recombination and mutation replace existing individuals if they are better than their predecessors. Two stopping criteria are considered; the maximum number of iterations or the maximum number of iterations without improvement of the objective functions.

3.4. Hybrid Genetic Algorithm

Hybrid methods (HGA) seek to combine the advantages of two or more metaheuristics for high quality solutions [20]. The most common hybrid methods combine population methods (e.g., Genetic Algorithms) with local search methods (Simulated Annealing and Variable Search Environment) or exact methods (linear programming and nonlinear) [21]. In this paper, a method population (GA) combined with a local search method (VSED) was implemented. The flowchart of the implemented algorithm is illustrated in Figure 1 below. The HGA structure retains essentially the GA structure described in the previous section; however, after mutation and before replacing the individuals of the next generation it performs a local search in order to find better quality individuals in the current generation [22,23]. As already described in the GA, only individuals exhibiting improvement in the objective function are included in the new generation [24].

4. Numerical Simulation Test and Investigation Results

In Figure 2, a distribution system of IEEE standard 34 bars is illustrated. The maximum network demand is 15.8 MW with a power factor of 0.95 lag. It may include up to 4 DG units, each one represented with a maximum power of 2.0 MW. Candidate solutions with DG units are penalized in the objective function (to minimize the loss in power) in order to make them less attractive. Note that according to the formulation (Equation (17)) the algorithm may select other number DG units up to the maximum.

The distribution of demand in each of the network nodes is shown in Figure 3. It can be seen that a significant part of the demand is in the first bars.

4.1. Coding Solutions

A chain of 4 binary bits was used to represent each used generator. The first bit represented the state of the generator (generator 1, the generator does not exist 0); the remaining three bits represented the power level of the generator. Thus, 1000 denoted minimum capacity (0.25 MW), while the chain 1111 represented a generator with maximum capacity (2.0 MW).

Table 1. Codification of Distribution Generation (DG) size.

Code	Size (MW)	Code	Size (MW)
1000	0.25	1100	1.25
1001	0.50	1101	1.50
1010	0.75	1110	1.75
1011	1.00	1111	2.00

illustrates the equivalence.

Table 2. Solution candidates encoding.

Bus 1	Bus 2	Bus 3	Bus 4	…	Bus nb-1	Bus nb
0100	0011	1011	0111	…	0011	1111

illustrates the code of a particular candidate solution comprising 4*nb bits, where nb was the number of bars in the system. Bars 1, 2, 4 and nb-1 have no generator; while bars 3 and nb have generating capacities of 1.0 and 2.0 MW, respectively.

4.2. Weighted Factors Calibration

Before starting the optimization process weighting factors must be assigned to the two objectives under study Improved Routine Protoco (IRPL) and intensity modulated proton therapy (IMPT) in order to assess the relative fitness of particular solutions. This requires a diagnosis without DG base case. This diagnosis was made by a load flow calculation with and without DG. The load flow analysis was performed by software Matpower [25] with DG modelled as bars where active power could be injected. In

Table 3. Values of IRPL and IMPT for candidates with randomly generated solutions.

Simulation Trial	IRPL	IMPT	IRPL/IMPT
1	41,884	1513	27,682
2	50,127	2806	17,864
3	43,297	2124	20,385
4	57,487	2478	23,199
5	34,329	1878	18,280
6	30,287	1079	26,070
7	44,238	2567	17,233
8	30,234	1969	15,355

, Figure 4 and Figure 5, the voltage profile and line losses, respectively, with and without DG are illustrated.

4.3. Results Using Simulated Annealing

To start the SA, a base solution that meets the criterion of maximum number of generators is generated. From there neighbouring solutions are explored according to the criteria stated in Equation (18). This means that early in the process the probability of accepting poorer quality solutions is high, but as the process evolves that probability is restricted. Thus, the principle of seeking diversity is privileged and at the end of the simulation it intensifies the search for better solutions. After performing several trial runs with different SA algorithm parameters, initial and final temperature was calibrated at 3.0 and 0.5, respectively, and at each iteration the temperature was reduced by 0.01. In

Table 4. Best result using simulated annealing (SA) .

Bar	Size (MW)	W1 IRPL	W2 IMPL	FO
8	1.5	29,578	43,747	73,325
17	1.75
22	2
29	1.5

and Figure 6 the best solution found (after 190 iterations) and the convergence of the algorithm illustrated by the SA. In this case the time for calculation was 1.2 min. A computer Intel Core i3 2.4 GHz with 4 GB of RAM was used in all simulations.

4.4. Results Using Variable Search Environment Descending

The size and location of the DG defined the VSED environment. A base solution of individuals was explored to better the environment by following the instructions as described in the methodology section. The search continues until it either passes a certain number of iterations or until no improvement is achieved over previous iterations. In

Table 5. Best result using Variable Search Environment Descending (VSED).

Bar	Size (MW)	W1 IRPL	W2 IMPL	FO
5	2.0	39,274	38,313	77,857
10	2.0
21	1.5
31	1.5

, the best found VSED is illustrated. Here 22 iterations were performed by evaluating the structure of eight individuals with a computation time of 2.1 min.

4.5. Results Using the Genetic Algorithm

The implemented GA solution used a coding chromosome structure as illustrated in Table 2. An initial population of possible solutions was generated pseudo-randomly which led to new solutions of better quality achieved by implementing a sequence of selection, recombination, and mutation as described in the methodology section. The GA parameters of initial population size, mutation rates, and recombination were calibrated through repeated runs. It was noted that with small initial populations, poor quality responses were obtained. By increasing the number of individuals in the initial population the quality of solutions improved, but the computing time increased. The best solution was found with a population of 100 individuals and mutation and recombination rates of 10% each. The computation time to find the best solution was 2.9 min after 90 iterations. The best result is shown in

Table 6. Best result using Genetic Algorithms (GA).

Bar	Size (MW)	W1 IRPL	W2 IMPL	FO
10	2.5	38,887	37,836	76,723
15	1.5
26	1.5
29	1.25

, while Figure 7 shows the process of convergence for different tests. It can be seen that even though the initial populations were of different quality, the tests converged to solutions of similar quality.

4.6. Results Using the Hybrid Genetic Algorithm

As mentioned above the HGA implemented combined the GA described in the previous section with VSED. That is, after applying the traditional operators of GA and prior to the population replacement in each iteration, the fitness of individuals was increased using VSED. The same parameters were calibrated as with traditional GA.

Table 7. Best result using Hybrid Genetic Algorithm (HGA).

Bar	Size (MW)	W1 IRPL	W2 IMPL	FO
12	1.7	40,573	39,201	79,774
13	2.0
22	1.0
30	1.25

presents the best solution found. It was observed for all tests that the computational time was considerably greater than that required by the other algorithms implemented; however, the quality of the response was better. The best answer was found after only six iterations in a computation time of 12.6 min.

5. Comparative Performance Index of Algorithms

The best result was obtained using HGA; however, this is the method that took longer to find a high quality solution. The fastest converging test algorithm was the simulated annealing, but solutions obtained with the other methodologies were of better quality. The VSED required less iteration to converge compared to the SA, but it’s computing time per iteration was greater and to be evaluated in terms of iteration, two structures of the neighborhood were needed. HGA generated the highest quality solutions when using large initial populations of 100 individuals. Figure 8 compares the computation time for the methods used, where the HGA stands out as relatively slow. The best solution found (see Table 7) showed an 80.2% reduction in losses compared to the base case without DG.

6. Conclusions

In this study, four algorithms were compared to obtain the optimum location and sizing of DG in distribution systems: SA, VSED, GA and HGA. The implemented methodologies based on these optimization algorithms were all successful in finding high quality solutions. It was observed that the optimum percentage of DG penetration of the test system varied between 6 and 7 MW installed at the start bar and the end bar for substantial improvement using DG units of 1–2 MW. In addition, the optimum location and sizing for DG allowed a substantial improvement of the voltage profile of the network and reduced losses by 80.2 per cent in the test system. In conclusive agreement with the presented research the results are dumped, which are emphasized by the three loss values:

(1)

For the optimized solution 134.7321 kW is the power generated and target achieved by the proposed hybrid genetic algorithm.

(2)

Losses for setting opening closer to the optimal solution, which is 134.9930 kW on fluctuation.

(3)

The value of a configuration that is not near the optimum and possibly a local optimum, 137.5293 kW. The ordinate represents the number of times it is achieved that value losses in five runs and the abscissa represent the fraction of crossbreeding.

It can be noted that the values determined by the range between 0.8 and 0.95 for the optimal solution search time is found to acceptable for the expected optimization margin. Furthermore, stress profiles were analyzed and found to have the lowest value of 0.9424 pu (31 bar), which implies a voltage drop across the bar of 6%, which is less than the maximum value allowed by national law: 8% in medium voltage. Contradictory to other countries, this value decreases to 5%. A new run was made taking into account this restriction, finding that there is no configuration having a bar with a drop level less than or equal to this voltage value. The configuration that approaches the performance of this restriction corresponds to the opening of the switches 7, 9, 14, 28 and 32 with a value of 134.9930 kW losses since the bar with less tension accumulated a voltage drop of 5.76% (31 bar).

Furthermore, the fastest algorithm in converging was the Simulated Annealing. Although, this methodology is faster, the solutions obtained with the other methodologies were of better quality. The VSED required fewer iterations to converge more than SA, but the computing time per iteration is greater, and to be evaluated by iteration, two structures of the neighborhood regarding GA it was found that the HGA is to find high quality solutions, large initial populations (100 individuals) must be used. Comparisons of the computation time for each of the methods was done, where the HGA had the best solution found, it was observed to have a reduction in losses of 80.2% compared to the base case without DG. Finally, the fastest method was SA; however VSED and the HGA showed better results as concluded in

Table 8. Comprehensive obtained performance between optimization algorithm.

Optimization Algorithm Type	Iterations Unit	Time in Minutes
SA	190	1.2
VSED	22	2.1
GA	90	2.9
HGA	6	12.6

.