An efficient load flow solution for distribution system with addition of distributed generation using improved harmony search algorithms

The distribution system analysis is very important to know the system condition at any given point of time. This distribution power flow plays an important role. Many power flow techniques have been proposed in the past for distribution power flow (DPF). The implicit Z-bus Gauss–Seidel (GS) method is good when dispersed generations are modelled as PQ nodes in the distribution system but, when dispersed generations are modelled as PV nodes, it leads to divergence problems. In a Newton-like method, a slight change in the execution of the bus-type switching logic leads to an unlike convergence characteristic. Though backward/forward (BW/FW) method is effective and efficient, it slows down under heavy-load condition. Now, with the addition of distributed generation (wind/solar, etc.) at the distribution level, the classical configuration of the distribution system is also changed. This poses a need for more robust and efficient technique. This paper investigates the application of a new algorithm known as improved harmony search (IHS) for DPF. The IHS can determine multiple power flow solutions which fix the limitations of other methods while determining a solution for the highly stressed condition. The results for a standard test system are taken with IHS, Newton–Raphson (NR) and BW/FW method without and with the impact of the wind energy system, to show the potential of the proposed algorithm.

solution technique. Moreover, power conversion components can be modelled as negative PQ load, ideal voltage source, PV node or its Norton/Thevenin equivalents. The steady-state behaviour is the most fundamental calculation of any system. In power systems, this calculation is the steady-state load flow problem, also known as power flow. Several Newton-like methods are developed to solve the problem in [1,2]. The first model of a backward/forward method for distribution load flow assumes all nodes as PQ since in the distribution system it is assumed that there is no generating bus. Sweeping need not to be calculated, and inversion of Jacobeans as Newton-like methods, resulting in fast convergence, is developed [3]. Firstly, probabilistic load flow is presented in [4] and further reformed and applied to power system load flow problems.
A modified Newton-Raphson (NR) load flow scheme is proposed by [5], which includes generator reactive power limits automatically. This method has slow convergence characteristic but converges monotonically. A slight change in the execution of the bus-type switching logic leads to an unlike convergence characteristic. The ladder network theory [6] and backward/forward method [7] are mostly adopted for load flow study of distribution system due to their computational efficiency and accuracies. An improved backward/forward method is proposed in [8], where the mismatch of a specified voltage and calculated voltage is tested after each backward sweep. But this method arises some limitation if the substation point (specified voltage node) lies in the middle or in between the network rather than on the one side of the network; then the algorithm becomes complex. Hence, the numbering of the nodes should be sequential. The reconstruction of bus injection to branch current (BIBC) and branch current to bus voltage (BCBV) matrix takes place spontaneously by just changing the two elements of the node-by-node sparse matrix (S), reflecting the change in network configuration. But the whole algorithm is applicable to balanced load and balanced system. Also, the convergence speed of the algorithm slows down with higher R/X ratio which is proposed in [9]. [10] presented the behaviour of the transformer in BW/FW-type load flow. The nodal admittance matrix is calculated by primitive admittance matrix using incidence matrix, but in a certain condition as why-delta, zero sequence component cannot update.
The harmony search (HS) algorithm has benefits like easy to implement and fast convergence but still has a disadvantage of precipitate convergence [11], which has been overcome with the help of improved harmony search (IHS) algorithm [12], solving a complex problem of distribution system reconfiguration using HS algorithm. The results are obtained on standard 33-bus and 119-bus system to validate the proposed algorithm. The results are compared with available standard references and found that the proposed algorithm is performing well. The HS algorithm is applied for the optimal allocation of distributed generator, and the total harmonic distortion (THD) analysis has been done to analyse the power quality. The proposed technique is applied on 12-bus unbalance distribution system and is showing that the placement of DG is fast as compared to other techniques [13].
In [14], a genetic algorithm is proposed to solve the power flow problem. An abnormal solution is obtained, which is not possible with the NR method. A genetic algorithm is proposed for distribution system load flow by [15] and also used for integrated AC/DC system load flow in [16]. A real-coded genetic algorithm for load flow solution is proposed [17]. For the placement of various wind distributed generators, a power system optimization (PSO)-based probabilistic load flow model is proposed in [18]. A probabilistic distribution load flow to study the effect of wind turbine generator connected to a distribution system is proposed in [19]. There are three different classes of wind turbine models considered in the distribution system as shown in Fig. 1.
In constant power factor wind turbine (WT), by controlling the WT trigger angle, reactive power can be adjusted and reactive power can be calculated by knowing the required power factor. In variable reactive power WT, reactive power can be formulated by known active power. Constant voltage WT is handled by considering the connected node as PV node.
In this paper, a constant power factor model is used for modelling the WT. The purpose of the paper is to show the potential of the new proposed algorithm. However, the algorithm is flexible enough, and any other WT model or distribution component model can be conveniently utilized.
The rest of the paper is organized as follows: "Distribution system modelling" section elaborates a brief modelling of the distribution system. "Wind turbine performance model" section describes the wind turbine model which is further considered as DG. The improved harmony search algorithm for power flow analysis is discussed in "Improved harmony search algorithm" section. The results of analysis without and with wind generation as DG are presented in "Results and discussion" section and the conclusions are given in "Conclusions" section.

Distribution system modelling
Distribution system components mainly include two parts: (1) shunt elements such as spot loads, distributed loads, and capacitor banks and (2) series elements such as transformers and line segments. Figure 2 shows a generalized series element model.  An equivalent representation of three-phase conductor with mutual coupling between ground and phases wires is shown in Fig. 3 [6,20]. And the line segment equations are given below In the admittance form Rearrange (2) and calculate I a, I b and I c . Same procedure can be applied to node V So, the three-phase π-model of a line can be written in the nodal frame as where Y abc = Z −1,abc . The similarity between single-phase and three-phase admittance matrix is that each element being exchanged by a 3 × 3 matrix. The current injection method can also represent the line charging (shunt capacitance). The charging currents are Figure 5 shows the capacitance in a three-phase line and its equivalent in current injection form.

Wind turbine performance model
The power output from a wind generator is a function of the speed of the wind. The power output changes with the change in wind speed. To compute the power output given from the wind turbine (WT), it is necessary to know about the wind speed and its (4) profile at that particular location in order to quantitatively determine the value of the profile that can be used in the modelling and simulation. A probability distribution function (PDF) is assumed, mathematically which can be represented by (5) [21].
where k is the shape parameter, c is the scale parameter and v is the wind speed. The cumulative distribution of Weibull distribution is mathematically represented by (6) [21,22] Once the uncertain and variable nature of the wind is determined and modelled as a random variable with the help of the probability functions, the output power of the WT may also be written as a random variable through a conversion from wind speed to output power. Ignoring some nonlinearity, the output power from WT can be determined from given wind speed by (7), where w =output power from wind turbine (typically in kW or MW); w r = rated power from WECS; v i =cut-in wind speed (m/s or miles/h); v r = rated wind speed (m/s or miles/h); v o = cut-out wind speed (m/s or miles/h). Figure 6 shows the wind speed Weibull probability distribution function for k = 2 and different values of c.

Power flow problem
The rectangular form of power flow problem can be represented as a set of nonlinear equations given below where E i and F i are the real and imaginary parts of ith node voltage. G ij and B ij are the (i, j)th element of bus admittance matrix and N = total number of buses.
For a PQ node, power mismatch can be given as where P i(spe) and P i(cal) are the specified and calculated power at ith node. For a PV node, active power and voltage mismatch will be where V i(spe) is the specified voltage at ith PV node and V i(cal) is the calculated voltage and is given by The power flow problem can be formulated as a minimization-type optimization problem with an objective function Fit as follows where N is the total number of buses. M is the total number of PV buses, considering first bus (i.e. i = 1) as a slack bus.
It is clear that the optimal value of the objective function will be zero (or on the brink of zero). Under a stress or heavily loaded condition when Fit doesn't go to zero, the method will try to minimize it as much as possible. The unknown variables are: (1) voltage magnitude ( |V | ) at all buses except slack bus and voltage PV buses with in specified limits V min , V max ; and (2) voltage angle (δ) at all buses within specified limits δ min , δ max .

Improved harmony search algorithm
In general, optimization techniques use gradient information to find the appropriate direction towards an optimal solution. However, a discrete variable cannot have the derivatives. To overcome this, Geem et al. [23] introduced a new meta-heuristic algorithm named harmony search (HS) that is a phenomenon mimicking conceptualized using the musical process of searching for a perfect state of harmony. This coordination in music is analogous to find the optimal solution to an optimization process. An instrumentalist always aims to produce a bit of music with perfect harmony. On the other hand, an optimal solution to an objective function of the optimization problem must be the best solution obtainable to the problem under the given objectives and limited by constraints. Both processes aim to produce the finest or optimum. The harmony search algorithm is used for network optimization [24]. The harmony search algorithm is used for optimal adjustment of PID controller parameter for an interconnected hydrothermal power system, and it shows that controller based on HS has comparatively less time of performance in comparison with genetic algorithm (GA) based controller [25]. The HS method can be successfully applied to various kinds of problems as structural analysis, mechanical component design, water distribution network, medical imaging, games and many others. But, its application in distribution load flow has not yet been investigated. The various steps of the improved harmony search algorithm are: Step 1: Problem formulation: In order to apply harmony search, an objective function should be formulated as where x i is the set of variables. X i = x L i , x U i , x L i and x U i are the lower and upper range of variables. N is the number of variables.
Step 2: Algorithm parameter setting: Once objective function is ready, the parameters of the algorithm should be set with convinced values as harmony memory size (HMS), harmony memory considering rate (HMCR), pitch adjusting rate (PAR), number of improvisation (NI), and fret width (FW). Random length only for a continuous variable is FW, which is generally called as bandwidth (BW).
Originally fixed values of PAR and BW were used. This fix value leads to higher number of iteration, resulting in slow convergence of the algorithm. There is a varying value for PAR and BW with iterations [26] as: where PAR min and PAR max are the minimum and maximum pitch adjusting rate. BW max and BW min are maximum and minimum bandwidth. NI is the number of solution vector generation, and I is the generation number. This modification of PAR and BW called the improved harmony search. This is utilized in this paper.
Step 3: Harmony memory initialization: Musician's harmony memory (HM) is an improvisation of many random harmonies. The minimum number of random harmonies should be HMS Step 4: Harmony improvisation: A new harmony is generated based on three basic operations expressed as follows: (1) Random selection: The new harmony is generating by randomly picking any value from the total value range x L i ≤ x i ≤ x U i with probability of (1-HMCR).
( Step 5: Update memory: In terms of objective, the new better harmony is replaced by worst harmony in HM.
Step 6: Termination criteria: Process stops when the maximum number of improvisation reached. Figure 7 shows the stepwise flowchart of the improved harmony search algorithm.

Results and discussion
The studies are performed on a standard 10-bus radial distribution system data [29]. The 10-bus radial distribution system consists of 9 load points. The single-line diagram of 10-bus section of the distribution system with wind generation at bus 10 is shown in Fig. 8.

Wind power based on wind distribution
The probability of wind speed is calculated as a Weibull curve defined by average wind speed and shape factor k. Generally, k = 2 for inland sites, 3 and 4 for coastal and island sites, respectively [30]. The wind is categorized into "bins" of 1 m/s. The instantaneous wind turbine power for each bin is multiplied by Weibull wind probability giving the net average turbine output power in that bin. The turbine average output power on a continuous 24-h basis is the sum of these powers. (20) : . . . : :  Assuming 11.55 m/s, rated speed of wind turbine and cut-in and cut-out speed is 3 m/s and 20 m/s. The rated power output power of a single turbine is 10 kW, and 75 similar wind turbines are assumed in a farm. Figure 8 shows the wind turbine power (kW) for each bin, Weibull wind probability in percentage and net output power of the turbine. The continuous average power comes 3.07 kW for each wind turbine. So, the total capacity of the wind farm is about 230 kW as shown in Fig. 9.

DLF without wind generation
For the system explained above, the distribution load/power flow studies are conducted by the three methods: (1) NR in rectangular form; (2) BW/FW sweep method [8] and (3) IHS algorithm.
In the study, the slack bus voltage and angle are fixed to 1 per unit and 0°. All buses are considered as PQ bus except one slack bus (bus 1). Table 1 shows the comparative study of voltage (per unit) and angles (°) for base case by all the three methods. The voltage magnitude and angles at various buses remain same with NR and BW/FW. Hence, they are represented by a single column.
Uniform load increase scenario has been utilized to increase the load as where i = 1,2,….n; n = number of load buses. P i Dbase and Q i Dbase are the base case active and reactive power demand at ith bus. λ is the load parameter. The uniform load increase scenario has been utilized for simplicity. However, the user can take any practical load increase scenario.
By gradually increasing the loads, Table 2 shows the voltage (per unit) and angles (°) of 10-bus radial distribution system by all the three methods. Again, voltage and angle (22)  remain very nearly same for NR and BW/FW method for loading parameter λ = 0.5. Again increasing the load, the Newton method fails to converge at λ = 1. Table 2 shows the voltage magnitude and angle for λ = 0.5 and 1. Again, increasing the load, both the methods (NR and BW/FW) fail to converge at loading parameter 1.2. Table 3 shows the voltage (per unit) and angles (°) profile with loading factor 1.2 by IHS algorithm.

DLF with wind power generation incorporated
By the loss-power sensitive, the appropriate location is found at bus 10. So, a total 230 kW wind generation is considered at bus 10. In this work, wind generation is modelled as negative PQ load. The load/power flow studies are again conducted by the three methods: (1) Newton method (NR) in rectangular form; (2) backward/forward (BW/FW) sweep method; and (3) improved harmony search (IHS) algorithm.  In the study, again the slack bus voltage and the angle are fixed to 1 per unit and 0°. Table 4 shows the comparative study of voltage (per unit) and angles (°) for base case by all the three methods. By gradually increasing the loads, Table 5 shows the voltage magnitude (per unit) and angles (°) for loading parameter, λ = 0.5 and 1.2. The voltage and angle are remained very nearly same for NR and BW/FW method at loading 0.5. At loading 1.2, Newton's method fails to converge.
Again, increasing the load, both the methods (NR and BW/FW) fail to converge at loading parameter 1.5, which has earlier failed at loading 1.2 in "DLF without wind generation" section, it shows the improved load handling capability due to the presence of wind generation. So, due to wind generation, the convergence capability shifts to higher loading side. Table 6 shows the voltage and angle profile with loading parameter 1.5 by IHS algorithm.

Conclusions
This work represents the application of improved harmony search algorithm for solving distribution system power flow problem. The algorithm is applied to standard 10-bus distribution system, and the results are compared with conventional NR and BW/FW method. The advantage of the proposed IHS method over conventional methods is its robustness. In many cases, during very heavy loaded condition the NR and BW/FW method fails to provide DLF solution. This method can be used to find the DLF solution in such situations where the NR and BW/FW method fails. The small-scale distributed generation is an attractive architecture for new distribution networks. Hence, in this paper wind generation system is also included in modelling DLF. The model concludes the adequacy analysis of distribution system power flow with significant wind generation. The results for standard 10-bus system demonstrate the potential of the proposed algorithm.