Active power ouptut optimization for wind farms and thermal units by minimizing the operating cost and emissions

Received Jul 18, 2019 Revised Jan 30, 2020 Accepted Feb 10, 2020 In recent years, many works have been done in order to discuss economic dispatch in which wind farms are installed in electrical grids in addition to conventional power plants. Nevertheless, the emissions caused by fossil fuels have not been considered in most of the studies done before. In fact, thermal power plants produce important quantities of emissions for instance, carbon dioxide (CO2) and sulphur dioxide (SO2) that are harmful to the environment. This paper presents an optimization algorithm with the objective to minimize the emission levels and the production cost. A comparison of the results obtained with different optimization methods leads us to opt for the grey wolf optimizer technique (GWO) to use for solving the proposed objective function. First, the method used to estimate the wind power of a plant is presented. Second, the economic dispatch models for wind and thermal generators are presented followed by the emission dispatch model for the thermal units.Then, the proposed objective function is formulated. Finally, the simulation results obtained by applying the GWO and other known optimization techniques are analysed and compared.


INTRODUCTION
Recently, the integration of the wind energy into electrical grids has increased significantly because it is clean and cheap in comparison to conventional energy sources. Therefore, in case of integration of wind farms (WF) with an existing grid with conventional sources, it is necessary to take also into account the wind power plants (WPP) in the economic dispatch (ED). In literature, many works have been done in order to discuss diverse methods of economic dispatch in which WFs are integrated into power systems.
In fact, in [1], an economic dispatch for combined wind thermal systems is realized using flower pollination algorithm (FPA). Besides, the authors in [2], present a new approach for ED problems in which WFs are installed in the network using particle swarm optimization (PSO) technique. Also, the authors in [3], discuss an economic load dispatch (ELD) in which solar and wind units are included in addition to thermal units. To resolve the ELD, the firefly algorithm optimization technique is used.
Nevertheless, although the fact that the wind energy is clean, the thermal power sources are considered as pollution sources. As a result, even the emissions have to be taken into account during the planning stage. In this paper, we propose an objective function that allows getting the optimal real power generation of the units which consist of thermal generators and WPPs. The proposed objective function aims to reduce both the production cost and emissions simultaneously. In this work, the potential power of a wind farm is considered as a constraint of the proposed objective function so that to not exceed the available wind power. The grey wolf optimizer (GWO) method is utilized in order to solve the proposed optimization algorithm.

WIND POWER ESTIMATION
The probability density function (pdf) of the weibull distribution is used to estimate the wind speed in any place. The probability density function f v (v) of the wind speed is expressed as follows [4]: where: c : is the scale parameter k : is the shape parameter v : is the wind speed The method cited in [5] is used for modelling the wind generator curve, as it is shown in Figure 1. This method is chosen because it is considered to be the most simplified method to simulate the power output of a wind generator [6]. The power provided by a wind turbine can be represented as [5]: where: v c : is the cut-in wind speed of the wind turbine v f : is the cut-off wind speed of the wind turbine v r : is the rated wind speed of the wind turbine Figure 1. Power curve of a wind turbine [7] The transformation of the probability distribution function pdf of the wind speed to the wind power can be expressed as follows [4]:

WIND POWER POTENTIAL OF A SITE
The wind power potential is given by [8]: where: A : Swept area of the rotor blades in m 2 ρ : Air density (kg/m 3 ) and is calculated as follows: with: ρ 0 = 1,225 kg/m 3 is the air density value at sea level and H m is the site elevation in m.
In this work, it is considered that H m is equal to 343m. Γ : gamma function The standard gamma function is expressed by [9]:

ECONOMIC DISPATCH
The goal of the economic dispatch is to minimize the operating cost of generators contributing to provide the load demand. In this study, the valve point effect is taken into account in the economic dispatch. Hence, the operating cost of each conventional generator can be given by the following [10]: where: p i : is the generation output of generator i a i , b i , c i : are the cost coefficients of generator i e i and f i : are the cost coefficients of generator i reflecting valve point effects The operating cost of the WPP is the sum of three components which are the direct cost, the penalty cost and the reserve cost. These three costs are explained in detail below. The system operator must pay a wind power generation cost C w i to the wind producer which may not exist if the wind power plants are owned by the power operator [4].
where: w i : is the scheduled wind power from the wind power generator i d i : is the direct cost coefficient of the wind power generator i If the scheduling of the wind power is less than it would be due to an underestimation of the available wind power, a penalty cost C p will appear that can be expressed as [4]: where: k p i : is the penalty cost coefficient for the wind power generator i f w (w) : is the weibull distribution function for wind power w r i : is the rated wind power from the wind power generator i The penalty cost may be equal to zero if the wind farm is owned by the power system. If the scheduling of the wind power is more than it would be due to an overestimation of the available wind power, a reserve cost C r i will appear that can be given by [4]: where: k r i : is the reserve cost coefficient of the wind power generator i The total operating cost is expressed as follows: (11) where: X : the position of each particle having the dimension (n=n1+n2) n1 : the number of the conventional generators n2 : the number of wind power generators

EMISSION DISPATCH
The emission dispatch aims to allocate the optimal power outputs for various generator units with the objective to minimize the emissions. In this study, the emission function is represented by the quadratic function [11]: (12) where: α i , β i , γ i : The emission coefficients of generator i

PROBLEM FORMULATION
In this work, the objective function aims to find the optimal allocation of power output from a combination of conventional generators and wind power plants by minimizing both the operating cost and the emissions. The proposed objective function is expressed as follows: where: w 1 , w 2 : are weight's factors, where w 1 and w 2 ≥ 0 and w 1 +w 2 =1 [12] h T : is the total price in penalty factor PPF in $/kg The total PPF is expressed as the summation of PPF of thermal units [13]: where: h i : is the price penalty factor of the generator i In this paper, we use the equation below to calculate the PPF of the i-th generator because it was proved that the min/max PPF is better than other PPF types [13]: The objective function proposed in (13) is subjected to the following constraints:  Operating limits of the conventional and wind power generators: where: p i max : is the maximum limit of real power output of the conventional generator i p i min : is the minimum limit of real power output of the conventional generator i w r i : is the rated wind power of the wind power generator i  The estimated available wind power of the generator i which is obtained using (4):  Power balance The total power generated from the conventional and wind power generators should be equal to the total demand P D plus the network losses P L .
where the power losses are calculated using B matrix technique and is expressed by the [14]: where B ij , B 0i and B 00 are coefficient of transmission loss. In this work, we use the B matrix cited in [15].

OPTIMIZATION ALGORITHOMS 7.1. Grey wolf optimizer algorithm for solving the economic and emission dispatch
In this work, the GWO technique which is a new metaheuristic optimization method is proposed to solve the economic and emission dispatch problem due to its many advantages such as being simple in principle, having a fast seeking speed and being easy to realize [16]. Besides, in [17], twenty-nine test functions were used to benchmark the performance of the GWO in terms of exploitation, exploration, convergence and local optima avoidance. The results illustrateted that GWO could provide highly competitive results compared to well known heuristics methods. The GWO method is invented by Seyedali et. al. in 2014 [17]. Grey wolves mostly live in a group called a pack. They are known by having a very strict social dominant hierarchy. In fact, in the hierarchy, the wolves within the pack are divided into alpha, beta, omega and delta. The first level of hierarchy of grey wolves is alpha. The alpha is responsible for making decisions about hunting and its decisions are dictated to the pack. Also, the alpha is considered to be the dominant, so his orders should be followed by the pack. The second level of hierarchy is beta. The beta wolves help the alphas in decision making. The delta wolves are in the third level of hierarchy. The delta have to obey to alpha and beta; however, they dominate the omega. In the lowest level, there are omega wolves that have to obey to all the other dominant wolves in the pack. In hunting, there are three main steps followed by a grey wolves: a) Finding, chasing, and approaching the prey. b) Encircling and harassing the prey until it stops moving. c) Attaking the prey.
In the grey wolf algorithm, it is considered that the alpha α is the fittest solution, while beta β and delta δ are considered respectively as the second and the third best solutions. The omega wolves follow the other three wolves. The equation to update the position of the pack is [17]: where A ⃗ ⃗ and C ⃗ are coefficient vectors, t is the current iteration, X ⃗ ⃗ p is the position vector of the prey, and X ⃗ ⃗ is the position vector of the grey wolf. The vectors A ⃗ ⃗ and C ⃗ are calculated using the [17]: where r 1 and r 2 are randoms verctors in [0,1], and componenets of are linearly decreased from 2 to 0 during iterations. The position vector of the grey wolf is updated based on the positions vectors of the first three best solutions which are α, β and δ. In this regard, the following equations are used [17]: The steps for solving the economic and emission problem using the GWO algorithm are as follows: Step 1. Generate randomly the grey wolf position of each search agent x i k .
Step 2. Evaluate the objective function for each search agent (13).
Step 3. Initialize X ⃗ ⃗ α , X ⃗ ⃗ β and X ⃗ ⃗ δ . X ⃗ ⃗ α is the position of the first best solution and X ⃗ ⃗ β and X ⃗ ⃗ δ are respectively the positions of the second and third best solutions. Step 5. Check that the constraints are satisfied (16), (17) and (18). Step 6. Evaluate the objective function of each search agent (13).
Step 8. If the maximum iteration is not reached, return to step4. Otherwise, stop the algorithm.

Particle swarm optimization technique (PSO)
The particle swarm optimisation (PSO) technique is a metaheuristic algorithm invented by Kennedy and Eberhart in 1995 [18]. In this method, the initialization of a group of particles is done in a random manner in the d-dimensional search space, where d is the size of the decision variables in the optimizatin problem. To each i-th particle a position vector x i , a velocity vector v i and a position Pbest i are associated. The particles exchange effectively information during an iterative process so that to find the optimal solution. In each iteration, the PSO algorithm searches for the optimal solution by updating the velocity (28) and the position (29) of each i-th particle taking into consideration its previous best position Pbesti and the best position of the group gbest. At each iteration k, the equations allowing to update the velocity and the position of each particle are given by [19]. v i k+1 = w k v i k + c 1 r 1 (Pbest i k − x i k ) + c 2 r 2 (gbest k − x i k ) (28) where: r 1 and r 2 : Uniformly distributed random numbers in the range [0 1]. w : Inertia weight. c 1 and c 2 : Acceleration coefficients.
The following equation is used to calculate the inertia weight [19]: where: k : the current iteration. k max : the maximum number of iterations. w min and w wax : are the lower and the upper bounds of the inertia weighting factors, respectively.

Bat algorithm (BA)
Bat algorithm (BA) is a metaheuristic optimization algorithm invented in 2010 by Yang [20]. In this optimization method, the echolocation behaviour of bats is used [20]. In order to look for prey, bat fly where: u (0,1) : is a uniform random number ranging from 0 to 1 x Gbest : the best solution found by the swarm A random walk is used for a local search and is expressed by [21]: where: φ : is a scaling factor allowing to limit the step size of the random walk A i : is the loudness N (0, σ) : is a normal random number with a standard deviation σ and a mean equal to zero when the bats are near their target, they decrease the loudness A i and increase the pulse rate r i . This can be expressed by [21]: where α and γ are constants.

Gravitational search algorithm (GSA)
The GSA is a heuristic optimization algorithm invented in 2009 by Rashedi et al. [22]. GSA technique is based on the law gravity and motion. In this method, agents are considered as objects and their masses are used in order to measure their performance. The position of i-th agent is defined by [23]: where: x i d : is the position of ith mass in the dth dimension n : is the dimension of the search space At a specific time t, the force acting from mass j to mass i is given by the following [23]: where: M i and M j : are the masses of the objects i and j R ij : is the euclidiean distance between the objects i and j ε : is a small constant G : is the gravitational constant at time t First G is initialized with G0 and then will decrease according to the time using in [23]: where α is a descending coefficient, t is the current iteration and T is the maximum number of iterations.
The total force acting on the agent i in the d dimension is given by the following in [23]: where rand j is a random number in the interval [0 1]. The acceleration at time t of the agent i in the d dimension is expressed by [23]: The velocity and the position of each agent is updated according to the [23]: where:

SIMULATION AND RESULTS
In order to test the proposed optimisation algorithm, three different cases are simulated using the IEEE 30 bus system. The test system is constituted of five conventional plants and one wind farm. The five conventional thermal plants are at buses 1, 2, 22, 27 and 23. While, the WF is at bus 13. Table 1 illustrates the cost and emission coefficients and the maximum and minimum power output of each conventional generator in the test system [13]. The wind farm utilized in this paper has a capacity of 120 MW. It is considered that the wind turbine technology installed in the WPP is V90/3000 (Vestas) [24]. Therefore, the rated power and the swept area of the rotor blades of the wind turbines installed in the WPP are respectively 3 MW and 6361.7 m 2 . Besides, Table 2 shows the values of the wind farm production limits (WF min and WF max ), the direct cost coefficient, the reserve cost coefficient and the penalty cost coefficient of the WF in the test system. In the three cases, the load system is considered to be equal to 300 MW. Also, the value of the shape k and scale parameter c are considered to be respectively equal to 2.14 and 7.29 of tangier region cited in [25].  Case 1: the weight's factors are considered to be w 1 =1 and w 2 =0.  Case 2: the weight's factors are considered to be w 1 =0 and w 2 =1.  Case 3: the weight's factors are considered to be w 1 =0,5 and w 2 =0.5.
In each case, 10 runs are done using the GWO and three other metaheuristic optimization algorithms: the PSO, the BA and the GSA which are used in the purpose of comparison. The optimization parameters are shown in Table 3. In each case, the average values of the active power outputs, the losses, the total operating cost and the emissions obtained in the ten runs for each optimization method are  Table 4. As we can notice from Table 4, when the only objective to minimize is the total cost, the GWO method attained better results than PSO, BA and GSA. In fact, the GWO method permits achieving the lowest average value of the total operating cost.
In the second case, when the only objective to minimize is the emissions, the GWO performs better than the other optimization methods by producing the best results. Indeed, the average values of the emissions obtained using the GWO is 280.6020 kg/h, while the average values obtained with PSO, BA and GSA are respectively 280.7096 kg/h, 301.0234 kg/h and 299.9055kg/h. In the third case, when the objective is to reduce both the operating cost and the emissions, the GWO permits getting the lowest value of the operating cost in comparison to PSO, BA and GSA. Regarding the emissions, the PSO attained the best results. Neverthless, the GWO permits getting better results than the BA and the GSA. In fact, the average values of the emissions obtained with the GWO, BA, and GSA are respectively 293.9747 kg/h, 314.7205 kg/h and 301.7588 kg/h.
According to the simulation results, the maximal real power provided from the WF is 72.5881 MW which represents the wind power potential of the site. Therefore, considering the wind power potential of a site as a constraint of the multi objective function proposed in this paper, permits not exceeding the available wind power. In addition, the GWO method can be used to solve such optimization problems due to the fact that it allows getting better results by reducing both the operating cost and the emission levels simultaneously in comparison to some known optimization algorithms, as we can notice from the simulation results.

CONCLUSION
In this paper, we present an optimization algorithm based on grey wolf optimizer (GWO). The proposed algorithm allows obtaining the active power output of thermal and wind power plants installed in the grid with the objective to reduce the cost production and emission levels simultaneously. Three different cases using the GWO, PSO, BA and GSA are analysed and compared. The proposed optimization method allows minimizing significantly the emissions and the operating cost as illustrated in the simulation results.