A Bio-inspired Heuristic Algorithm for Solving Optimal Power Flow Problem in Hybrid Power System

In recent studies, emphasis has been placed on optimal power flow (OPF) problems in traditional thermal, wind, and solar energy sources-based hybrid power systems. Various metaheuristic algorithms have been proposed to find optimal solutions to the OPF problems in the hybrid power system. The OPF, due to the quadratic nature of its primary objective function, is a nonlinear, nonconvex, and quadratic optimization problem. In this study, we have proposed a bio-inspired bird swarm algorithm (BSA) to find an optimal solution to the OPF problem in the hybrid power system because it performs well in the case of optimizing the well-known Rastrigin quadratic benchmark function. In this study, uncertainty of utility load demand and stochastic electricity output from renewable energy resources (RESs) including wind and solar are incorporated into the hybrid power system for achieving accuracy in operations and planning of the system. We have used a modified IEEE-30 bus test system to verify and measure the performance of BSA and a comparison is made with well-known evolutionary metaheuristic algorithms. The proposed BSA consistently achieves more accurate and stable results than other metaheuristic algorithms. Simulation-based optimization results have shown the superiority of BSA approach to solve the OPF problems by satisfying all constraints and minimum power generation cost 863.121 $/h is achieved in case study 1. Simulation-based experiment results have indicated that by imposing the carbon tax (ton/h) the power generation from RESs was increased. In case study 2, the proposed BSA approach has also outperformed and minimum electricity cost 890.728 $/h is achieved as compared to other algorithms.


A. BACKGROUND AND MOTIVATION
Generally, an electric power grid or power system consists of power generation plants, electricity transmission and distribution systems. The power generation plants generate electrical power through electricity generation units or power plants. The transmission system carries the electricity through transmission lines from electricity generation plants to utilities or load centers. The electricity distribution system feeds the electricity through distribution lines to nearby homes, agricultural units, industries, and commercial buildings. The traditional electric power grids are responsible for producing electricity and carrying it to residential, industrial, and commercial consumers through electricity transmission and distribution lines. Two authorities, 1) independent system operator (ISO) and 2) electric utility, are responsible to control operations and planning of the power system in a country. The ISO is an independent authority established by the government to ensure reliability of electricity generation and the transmission system in the electrical power grid. An electric utility is an authority that engages in feeding the electricity through distribution lines to consumers by balancing the demand and supply of the electrical load.
In the electrical power grid operations and planning, the VOLUME 4, 2021 control always resides on the generation side and the power generation plants adjust their electricity generation according to the changes in electricity demand from consumers. Sometimes power generation plants produce surplus electricity, which is transmitted to the nearby area by transmission lines or stored [1]. Therefore, it is of practical importance to balance load demand and electricity supply in the power system. For this purpose, many techniques have been applied in the research literature. On the generation side, to address optimal power flow (OPF) problems in the power system is considered as a technique for finding stable and secure operating points of electricity generation plants and their optimal scheduling on an hourly basis [2][3][4][5].
In 1962, Carpentier first introduced economic dispatch problem extension as the OPF problem in traditional thermal energy sources-based power systems [6]. The OPF is one of the well-known and well-studied research areas in the power system. It can be defined as: "To find out the stable and secure operating points (levels) for electricity generation plants in order to meet load demand of utility in power system, generally with attention to minimize electricity generation cost" [6]. In traditional thermal energy sources-based power systems, the OPF is a nonconvex, nonlinear, and quadratic problem due to the quadratic nature of its primary objective function to reduce electricity generation cost. The primary objective function of OPF problem has been modeled as quadratic curve and its various forms such as valve-point loading effect quadratic curve, piecewise quadratic curve, and prohibited operating zones quadratic curve for the traditional thermal energy source [6,7]. Researchers have also proposed various techniques for solving the OPF problems considering other objectives, in addition to the primary electricity generation cost minimization objective. These objectives include minimizing voltage deviation, power loss in transmission lines, and emission pollution and enhancing voltage stability index [6][7][8][9].
In the last decade, integration of environment friendly and clean electricity output from renewable energy sources (RESs) including wind and solar into thermal power systems have become necessary due to the rising demand for electricity and global warming issues. Therefore, the power systems are striving towards a sustainable system future due to rapidly growing integration of RESs in power systems. On the electricity generation side, the RESs such as solar photovoltaic (PV) units and windfarms are being owned by private parties in a power system. The ISO purchases scheduled renewable electricity from private parties in order to cater the growing consumers' load demand. Wind power generation depends upon stochastic wind speed at different times of day. Similarly solar PV power generation depends upon uncertain solar irradiance during the day time. Due to the fluctuant and intermittent solar and wind power output, the available power from solar PV units and windfarms may be more or less than wind-scheduled power at different times of day. In an overestimation scenario, the ISO is required to have a spinning reserve based on utility load demand, when power supplied by solar PV units and windfarms operators is less than wind-scheduled power. The ISO has to increase the reserve cost associated with reserve electricity generation units to balance the supply and demand in this scenario. An underestimation scenario may arise when actual renewable energy received from RESs is greater than scheduled power. In that case, the surplus power output from RESs is wasted and ISO bears a penalty cost if it is not stored or transmitted to a nearby area [8,9]. Incorporating stochastic power generation from wind and solar into the system raises the complexity of power system operations and planning. The utility load demand is also uncertain in nature due to variation in consumers' load demand that directly affects spinning reserve cost in the power system. Moreover, considering the uncertainty of utility load demand has significant importance to achieve accuracy in the operations and planning of the system. Therefore, an effective technique is required to reduce the overall electricity generation cost.

B. LITERATURE REVIEW
In the research literature, various studies have been documented using two types of optimization algorithms for solving the OPF problems in the power system. These optimization algorithm types are traditional mathematical algorithms or methods and metaheuristic algorithms. Numerous mathematical optimization methods including linear programming [10], linear/quadratic programming [11], sequential linear programming [12], newton method [13], generalized benders decomposition (GBD) [14], nonlinear programming [15][16][17], mixed integer nonlinear programming (MINLP), [18], interior point method [19,20], and simplified gradient method [21] have been applied to solve the OPF problems. In these traditional methods, nonlinear objective function and constraints are converted into linear form before solving the OPF problem because the mathematical method cannot handle the nonlinear properties of the problem [22]. This convergence in constraints and objective functions may affect the accuracy of operations and planning of the power system.
The OPF problem primary objective -electricity generation cost minimization, is considered in all of the aforementioned studies [6,[23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. Moreover, other objectives such as reducing power loss in transmission lines, emission pollution, etc. are considered in some studies. In these studies, the performance of proposed metaheuristic algorithms has been measured on one or more IEEE-30, IEEE-57 and IEEE-118, bus test systems. The power systems are rapidly growing with RESs integration due to increasing electricity demand and global warming issue due to traditional thermal energy sources. In the literature, some studies have been documented [38][39][40][41][42] to find optimal solutions to the OPF problems in traditional thermal and wind energy sources-based hybrid power systems with a focus on minimizing the overall cost of electricity generation. In study [38], gbest guided ABC (GABC) has been utilized to find optimal solutions to the OPF problems in the thermal and wind energy sources-based hybrid power systems. In which objectives were to minimize electricity generation cost and emission pollution. In study [39], modified bacteria foraging algorithm (MBFA) is employed for solving the OPF problem in the traditional thermal and wind energy sources-based hybrid power system. A doublyfed induction generator model is utilized to justify the OPF problem inequality constraints and problem is formulated with various objective functions in above study. In study [40], authors applied ant colony optimization (ACO) and MBFA for solving the OPF problems in the traditional thermal and wind energy sources-based hybrid power systems. In study [41], authors utilized multi-objective glowworm swarm optimization (GWSO) to solve the OPF problems in the thermal and wind energy sources-based hybrid power systems. In all of the aforementioned studies, [38][39][40][41], the authors have utilized a modified IEEE-30 bus test system to verify and measure performance of the applied approaches. In study [42], the authors have adopted self-adaptive evolutionary programming (EP) for solving the OPF problems in the traditional thermal and wind energy sources-based hybrid power systems. In all of the aforementioned studies [38][39][40][41][42], the authors have applied well-known Weibull probability density function (PDF) for modeling uncertainty of stochastic wind speed to incorporate wind power output into hybrid power systems.

C. PROBLEM STATEMENT
In study [8], authors have proposed success history-based adaptation differential evolution (SHADE) to find optimal solution to the OPF problems in hybrid power system. The OPF problem objectives -to reduce electricity generation cost and emission pollution are considered. In study [9], the author has proposed a fuzzy logic technique based on PSO finding optimal solution to the multi-objective OPF problems in hybrid power systems, by considering objectives to reduce active power output cost and power loss in transmission lines.
In both studies [8,9], uncertainty of stochastic solar irradiance and wind speed are incorporated into the power system to solve the OPF problems. However, the utility load demand uncertainty has been ignored in these studies [8,9]. On the other hand, SHADE and PSO may be inefficient to find optimal solutions to the nonlinear and quadratic OPF problems because DE has a premature convergence property [54] and PSO is incapable of searching neighborhood existing solutions in nonlinear quadratic optimization problems [55].

D. CONTRIBUTION AND PAPER ORGANIZATION
In study [56], authors have proposed a new bio-inspired bird swarm algorithm (BSA). In which, it has been observed that the BSA has good diversity and can flexibly regulate its four different search strategies such as foraging, vigilance, producer and scrounger to explore the search space. Moreover, BSA can improve its convergence speed without affecting the stability and accuracy of optimal solutions by making better balancing among exploration and exploitation of search space. In fact, under suitable interpretations, DE and PSO mutation operators are distinct forms of the proposed BSA approach. In which, the bird's social behaviour such as the scrounger formula is similar to the DE mutation operator and the foraging formula is similar to the PSO. Moreover, the BSA has prominent distinguishing features, in addition to the merits of the DE and PSO.
In study [56], optimization results have proved the superiority of the BSA as compared to DE and PSO to optimize the Rastrigin function F 9 , which is a well-known quadratic benchmark function for performance evaluation of optimization algorithms. The primary objective function of the OPF problem is to reduce power generation cost, which follows a quadratic nature function. Inspired by the study [56], we have proposed BSA to find an optimal solution to the OPF problem in the hybrid power system. This research study is an extension of [57] our conference paper already published. It have following knowledge contribution in academic research; • A method based on BSA is proposed for finding an optimal solution to the OPF problem in the hybrid power system by incorporating uncertainty of the utility load demand and stochastic power output from RESs. We have used a modified IEEE-30 bus test system to measure and evaluate the performance BSA for finding an optimal solution tothe OPF problem in the hybrid power system. Simulation based optimization results have shown the superiority of BSA approach to solve the OPF problems (see Table 10 and Figure 8).
The organization of paper is as follows: The uncertainty modeling of utility load demand, stochastic solar irradiance, and wind speed are described in section II. The section III consists of OPF problem formulation, objective function, and different constraints. In section IV, we have explained the proposed metaheuristic method. The simulation-based experiment results for the proposed BSA approach and other algorithms are specified in V and concluding remarks are written in section VI.

II. UNCERTAINTY MODELLING OF WIND SPEED, SOLAR IRRADIANCE, AND UTILITY LOAD DEMAND
The most significant aspect of uncertainty modeling is to use an appropriate PDF for predicting the values of uncertain or random variables. We have used the Lognormal PDF and Weibull PDF for modeling uncertainty of solar irradiance and stochastic wind speed by adopting same strategy proposed in study [8]. We have utilized the same formula and approach of Gaussian (normal) PDF presented in studies [59,60]for modeling uncertainty of utility load demand. The Hong's point estimate method (PEM) proposed in study [61] has been used for calculating the utility load demand on load buses.
In this section, first, we have described the model for handling the uncertainty of wind speed to incorporate wind power output. Secondly, we have explained the model for incorporating stochastic solar power in the power system. Lastly, we have described Gaussian PDF for handling uncertainty of utility load demand.

A. UNCERTAINTY MODELING OF STOCHASTIC WIND SPEED
In the research literature, the well-known Weibull PDF has been mostly applied for modeling the wind speed v (m/s) [38][39][40]42] to incorporate the stochastic nature wind electricity generation into the power system. The uncertainty modeling stochastic wind speed using Weibull PDF can be defined as: where, scale factor c and shape factor k are parameters and its mean (M wbl ) can be calculated as: where, Γ() represents the gamma function and it can be formulated as: For the performance evaluation of BSA, we have modified standard IEEE-30 bus test system, in which two traditional thermal energy sources at 5 and 11 generator buses are replaced with wind energy sources (windfarms). The total number of wind turbines (WTs) in each windfarm and values of c and k of Weibull PDF parameters are given in Table3.  In a windfarm, the actual power output of WT depends upon wind speed v (m/s) it meets. In both windfarms, each WT has a 3 MW rated power output. The cumulative rated power generations of windfarms connected at generator buses 5 and 11 are 75 MW from 25 WTs and 60 MW from 20 WTs, respectively. The WT electricity generation can be VOLUME 4, 2021 formulated as follows [8]: where, v in is cut-in, v out is cut-out and v r is cut-in rated wind speed met to WT and its rated power output is represented as P wr . We have considered v r = 16 m/s, v out = 25 m/s, and v in = 3 m/s for wind power output calculation purpose.
The histograms in Figure 2 indicate wind power output based on wind speed Weibull distribution plotted in Figure 1, from windfarms connected at bus 5 and bus 11. It is observed from WT power output Eq. 4 that WT provides rated power output P wr when wind speed is in range [v r v out ]. The wind power output is calculated as [8]: The WT provides wind power in continuous form when wind speed is in range [v in v r ] and wind power output for this region is measured as [8]: The probabilistic model for solar irradiance I (W/m 2 ) follows lognormal PDF and it is generally utilized for handling the stochastic solar power generation in hybrid power systems [58]. The uncertainty of solar irradiance I (W/m 2 ) mathematically can be written as [58]: where µ in mean of solar irradiance and σ is standard deviation of solar irradiance. Lognormal PDF mean (M lgn ) is defined as: Assumed values for mean (µ) of solar irradiance and standard deviation (σ) have been specified in Table 4. In a modified  Figure 3.
The solar PV unit available or actual solar power generation subject to solar irradiance I (W/m 2 ) and solar PV unit electricity generation modeled as [8]: where I std and I c indicate the solar irradiance and certain solar irradiance point in a standard environment, respectively. The solar PV unit rated power output is represented as P pvr .
In this study, we assumed I std = 800 W/m 2 and I c = 120 W/m 2 . The rated power outputs related to solar PV units that are connected at generator bus 8 and 13 are P pvr = 35 MW and P pvr = 50 MW, respectively. The histograms in Figure  4 represent the stochastic solar power output based on solar irradiance I (W/m 2 ) distribution shown in Figure 3, from solar PV units.

C. UTILITY LOAD DEMAND UNCERTAINTY MODELING
The utility load demand is also stochastic due to variation in consumers' load demand that directly affects spinning reserve cost in a power system. Therefore, modeling the utility load demand uncertainty has a significant impact on solving the OPF problem and achieving accuracy in planning and operations of the power system. In this study, utility active (real) load is considered as a random variable and power factor is considered as constant. According to the constant power factor, the change in reactive power of each load bus or PQ bus depends upon its active load in the power system. We have used a Gaussian PDF [59,60] to model the uncertainty of utility load demand on each load bus in transmission system.
The prediction of load demand on a load bus follows Gaussian PDF as: where, P d,i is active load demand, σ i is standard deviation and µ i is mean value of active load at i th load bus. We have utilized modified IEEE-30 bus system for performance evaluation of our proposed method, in which four thermal energy sources at 5, 8, 11, and 13 generator buses are replaced with RESs due to emission pollution and global warming issues. The down-arrow (↓) symbol represents load or PQ bus shown in Figure 5. The active (real) load demand mean value µ has been considered equal to the base active load of each load bus, and 10% of the mean value µ has been considered as standard deviations σ of active load on each load bus for modeling uncertainty of utility active load demand across the modified IEEE-30 bus test system. Plotted histograms in Figure 6 based on Gaussian PDF and represent the distribution of active load demands on load or PQ buses, which are achieved after execution of 8000 Monte-Carlo simulation scenarios.
In the research literature, many approximation methods based on the analytical approach have been documented [62] for handling the uncertainty in power systems. The common uncertain source method, the discretization method, the truncated Taylor series expansion method, and the point VOLUME 4, 2021 developed an efficient PEM to measure the moments of Z = h(x), where Z represents an uncertain or random quantity and it is a function of n uncertain variables. The PEM is a simple to use method for measuring the moments Z and does not require derivatives of h(x) or any iteration as compared to other approximation methods such as the discretization method and Taylor series expansion method. The PEM can be utilized directly with a deterministic computer program. Based on the above facts, we used Hong's PEM for calculating the approximate load on load buses.
Hong's PEM concentrates upon statistical information obtained through initial few moments of a random variable on m concentrations for every variable. In which m × n points concentration matching is considered to achieve first noncrossed moments (m×n) and crossed second-order moments of each uncertain variable. The simplest case of Hong's PEM is the 2m scheme, in which skewness of an uncertain variable's PDF and correlation between the uncertain variables are considered for predicting the value of random quantity. Another particular case of Hong's PEM is the 2m +1 scheme, in addition to the skewness of the uncertain variable's PDF and correlation between the uncertain variables, kurtosis of the uncertain variable's PDF is also considered in this scheme to calculate the random quantity. More detail of Hong's PEM is available in [61][62][63].
In this research work, we have utilized Hong's PEM [61] by taking Gaussian PDFs of active load demands on modified IEEE-30 bus test system load or PQ buses as input from plotted histograms in Figure 6, to calculate active load on each bus. The deterministic load on each load bus obtained from using particular schemes of Hong's PEM such as 2m and 2m+1 is specified in Table 5 and plotted in Figure 7. The deterministic active load demand on each load bus based on both 2m and 2m+1 schemes is similar. Therefore, we have used the simplest 2m scheme of Hong's PEM for calculating load demand on each load bus to find an optimal solution to the OPF problem in the hybrid power system.

III. THE OPF PROBLEM FORMULATION
The OPF in the traditional thermal energy sources-based power system is a quadratic nature nonconvex and nonlinear problem, in which stable and secure settings of operating points in electricity generation plants are obtained for minimizing certain objectives. The OPF problem objective written as [45]: s. t.: g(x, u) = 0 h(x, u) ≤ 0, where, u is a set of control variables and x is a set of state variables. The function f (x, u) represents the objective of the OPF problem. The function g(x, u) represents equality constraints and h(x, u) represents inequality constraints.
The power flow in system is controlled by control or independent variables, while state variables described power system state. The control variables consist of all bus generators active power excluding slack (swing) bus active power, all generators or energy sources voltage magnitudes, shunt compensator at selected buses, and transformer tap in power system network. The state variables consist of generators' reactive power, swing bus active power output, line loading of transmission lines, and voltages magnitude at load buses.

A. CONSTRAINTS
Balancing both active (real) power and reactive power follow equality constraints, while security constraints and equipment's operating limits of transmission lines and load buses output from all energy sources also must be equal to demand and loss of reactive power. In the power system, equality constraints can be written as according to study [45]: where, the term P Gi indicates the active power and term Q Gi represents reactive power generation from energy source attached at i th bus. The term P Di represents the load demand of active power and term Q Di is demand of reactive power at i th load bus. The term V i represents i th bus voltage magnitude, while term V j represents j th bus voltage magnitude. The term θ ij = θ i − θ j shows difference of voltage angles θ i and θ j at i th and j th buses. The term G ij is transfer conductance and B ij is transfer susceptance in i th and j th buses.

2) Inequality constraints
In power systems, load buses and transmission lines security constraints, secure and stable equipment's physical settings, and equipment's operating lower and upper limits are considered as inequality constraints. These are mathematically described as follows [45]: P min ws,j ≤ P ws,j ≤ P max ws,j j = 1, ..., N W F , P min pvs,k ≤ P pvs,k ≤ P max pvs,k k = 1, ..., N P V , Q min ws,j ≤ Q ws,j ≤ Q max ws,j j = 1, ..., N W F , Q min pvs,k ≤ Q pvs,k ≤ Q max pvs,k k = 1, ..., N P V , The active (real) power output boundary limits of traditional thermal energy source and RESs such as windfarms and solar PV units are represented in constraints (15)- (17) and N T G , N W F and N P V are total number of traditional thermal energy sources, windfarms, solar PV units, respectively. The reactive power output boundary limits of all energy sources or generators including traditional thermal energy sources, windfarms, and solar PV units are defined by constraints (18)- (20), respectively. On load buses and generator buses voltage magnitude boundary limits are defined by constraint (21) and (22), respectively. N G is number of generator buses . In each inequality constraints max and min superscripts show boundary limits of the corresponding parameter. N L is number of load buses.

B. OBJECTIVE FUNCTION
To reduce electricity generation cost from traditional thermal energy sources and RESs and including emission cost (e.g., carbon tax) is our objective for finding an optimal solution to the OPF problem in the hybrid power system. The minimization of electricity output cost objective is stated as follows: The first term ξ T (P T G ) represents power generation cost of thermal energy sources, which based on valve-point loading effects quadratic curve. The second term N W F j=1 ξ w,j (P w,j ) represents the cost of power output from windfarms. In objective function, the third term N P V k=1 ξ pv,k (P pv,k ) is cost of solar power output from solar PV units. The last term ξ E is the cost of emission (e.g., carbon tax). The detailed modeling of traditional thermal and RESs power generation cost functions and emission cost are provided herein.

1) Thermal power cost curve
The traditional thermal energy source required fossil fuel to produce electricity. The power generation cost of fossil fuel based energy sources can be calculated by regular quadratic fuel curve, valve-point effects quadratic fuel curve, and piecewise quadratic fuel curve [38]. In some practical cases, thermal power is generated from traditional thermal energy sources using different fossil fuels like natural gases,coal and oil. The power output cost from these types of energy sources is calculated using the piecewise quadratic fuel cost curve.
In this study, we assumed that the same fossil fuel based traditional thermal energy sources are used for power generation purposes. Therefore, to calculate the power cost related to thermal energy source, we used two forms of fuel cost curve; 1) quadratic fuel curve and 2) valve-point effects quadratic fuel curve. The generated power (MW) from thermal energy source followed a quadratic relationship with fossil fuel cost ($/h) as [45]: where, P T Gi is i th thermal energy source power output. a i , b i , and c i are i th thermal energy source cost coefficients. In a traditional thermal energy source-based power system, the objective function is modeled using valve-point loading effects quadratic fuel cost curve for more precise and realistic measuring of thermal power cost. Because traditional thermal energy source's steam turbines have multi-valve, in such case a variation befall in fuel cost curve. In such case, the fossil fuel cost curve of thermal energy source is measured using valve-point loading effects quadratic fuel cost curve as follows [45]: where, coefficients d i and e i represent the valve-point loading effect of i th TG. P min T Gi represents the minimum power output of i th TG during operation. For this study, cost coefficient values of fuel cost curve and valve-point loading effects related to T G 1 and T G 2 are specified in Table 6.

2) Wind power cost function
We have assumed that private parties hold RESs such as solar PV units and windfarms and ISO purchases scheduled power from private parties according to the signed agreement. The wind power cost follows directly proportional relationship to wind-scheduled power and can be formulated as: where, λ ws,j is a function to calculate power cost of windscheduled power P ws,j received from j th windfarm and g j is cost coefficient of wind-scheduled power related to j th windfarm. The distributions of windfarms power output are shown in Figure 2. The available power from windfarms can be less or more than wind-scheduled power because of fluctuant and stochastic wind power output. In an overestimation scenario, when wind power supplied by windfarms operators is less than the wind-scheduled power, the ISO is required to have a spinning reserve based on utility load demand. The wind power reserve cost λ wr,j for j th windfarm can be calculated as [39]: λ wr,j (P ws,j − P wav,j ) = K wr,j Pws,j 0 P ws,j − P w,j × f w P w,j dP w,j , (27) where, P wav,j , K wr,j , and f w (P w,j ) represent available wind power, coefficient of windfarm reserve cost, and wind power for j th windfarm, respectively.
In a scenario, when wind power supplied by windfarms operator is greater than the wind-scheduled power, if it not possible to reduce power generation from thermal energy sources, the windfarms surplus electricity is dumped and ISO bears penalty cost. The penalty cost λ wp,j related to j th windfarm can be formulated as [39]: where, K wp,j is penalty cost coefficient of j th windfarm.
The cost related to any windfarm power is calculated by adding penalty cost, reserve cost, and its wind-scheduled power cost. The cost coefficients and wind-scheduled power are specified in Table 7. The total wind power cost ξ w,j (P w,j ) for j th windfarm can be calculated by adding windscheduled power cost and both penalty and reserve costs as [39]: ξ w,j (P w,j ) = λ ws,j (P ws,j ) + λ wr,j (P ws,j − P wav,j ) + λ wp,j (P wav,j − P ws,j ).

3) Solar power cost function
Solar power cost is also directly proportional to solarscheduled power. The solar-scheduled power cost λ pvs,k is a function of scheduled power provided from k th solar PV unit, as follows [42]: where, scheduled power cost coefficient h k is related to k th solar PV unit and P pvs,k is delivered solar-scheduled power from k th solar PV unit. In Figure 4, the distributions of power generation from solar PV units are plotted. Similar to the windfarms power output behaviour, in the underestimation scenario, the available power from solar PV units can be more than solar-scheduled power and in an overestimation scenario, the available solar power can be less than solar-scheduled power. In such a case, the ISO requires a spinning reserve energy source. According VOLUME 4, 2021 to the concept presented in the study [42], we have modeled solar reserve cost λ pvr,k of k th solar PV unit power output, for an overestimation scenario. It can be formulated as [42]: λ pvr,k (P pvs,k − P pva,k ) = K pvr,k × f pv (P pva,k < P pvs,k ) × P pvs,k − E(P pva,k , < P pvs,k ) , (31) where, K pvr,k represents coefficient of reserve cost of k th solar PV unit, P pva,k is actual available power received by k th solar PV unit, and f pv (P pva,k < P pvs,k ) represents the calculation of overestimation scenario related to k th solar PV unit. E(P pva,k < P pvs,k ) represents the prediction (expectation) of power output from solar PV unit k th below solar-scheduled power P pvs,k .
For underestimation scenario, the penalty cost λ pvp,k related to k th solar PV unit power output can be formulated as [42]: where, K pvp,k represents the penalty cost coefficient for k th solar PV unit. f pv (P pva,k > P pvs,k ) represents the calculation of underestimation scenarios related to k th solar PV unit. E(P pva,k > P pvs,k ) represents the prediction (expectation) of power output from solar PV unit k th above solar-scheduled power P pvs,k .
Similar to windfarm power cost, the cost related to a solar PV unit power is calculated by adding penalty cost, reserve cost, and solar-scheduled power cost. The cost coefficients and solar-scheduled power related to solar PV units are specified in Table 8. The total solar power cost ξ pv,k (P pv,k ) for k th solar PV unit by adding solar-scheduled power cost and both penalty and reserve cost can be measured as [42]: ξ pv,k (P pv,k ) = λ pvs,k (P pvs,k ) + λ pvr,k (P pvs,k − P pva,k ) + λ pvp,k (P pva,k − P pvs,k ).

4) Emission cost
The combustion of fossil fuels in traditional thermal sources of energy is the core cause of greenhouse/harmful gases including SO x , CO x , and N O x emission into the environment.
With growing global environmental concerns, to regulate the power system for accounting the minimum emissions is necessary. The total emission E (ton/h) into the environment by a traditional thermal energy source can be formulated as follows [6]: where, emission coefficient for i th TG are α i , β i , γ i , ω i , and µ i and values of these coefficient related to T G 1 and T G 2 are listed in Table 6.
In the recent decade, due to global environmental issues, many countries are imposing a carbon tax to minimize carbon emission into the environment. Therefore, carbon tax widely has been applied to curb greenhouse gases and encourage investment in clean forms of energy, [59]. Carbon tax (C tax ) imposed on emitted greenhouse gases and emission cost ($/ton) can be calculated as:

IV. PROPOSED APPROACH TO SOLVE THE OPF PROBLEM
Various nature-inspired metaheuristic algorithms have been developed as a substitute to the mathematical methods for solving optimization problems, in research literature.
Population-based BSA [56] is a new stochastic swarm intelligence algorithm. To address the optimization problems, intelligence of bird swarms extracted from bird's social behaviours has been utilized. The birds in the swarm improve their fitness through social behaviours and interactions with other birds in the swarm. The working model of BSA is based on three types of bird behaviours such as foraging, vigilance, and flight behaviour. The BSA can flexibly regulate its four different search strategies such as foraging, vigilance, producer, and scrounger to explore the search space. Based on these facts, the BSA can improve its convergence speed through better balancing between exploitation and exploration of search space without affecting the stability and accuracy of the optimal solution. Therefore, the proposed method based on BSA may provide a more stable and accurate solution for the OPF problems in the hybrid power system. A bird's social behaviours and interactions with other birds in the swarm can be made understand based on welldefined rules as follows: Rule 1: Individual birds in a swarm may switch into two types of behaviours; 1) foraging behaviour and 2) vigilance behaviour, on the basis of the random or stochastic decision. Rule 2: In a swarm, the individual bird may update or improve fitness through social behaviour and by promptly recording self and swarm's best memory or previous experience to explore the food patches in a specific area during foraging behaviour. The bestrecorded experience or memory about searching food items can be utilized to explore food patches, and social behaviour and information are shared immediately. food reserves and a bird having greater food reserves than other birds will be at the swarm's center. Rule 4: The birds in a swarm have flight behaviour due to foraging behaviour or any other reason. During the flight behaviour the birds in the swarm can be often switched again into two types of birds; 1) producer birds and scrounger birds on the basis of their food reserves. Birds that have food reserved between lowest and highest are randomly switched into scrounger and producer. Rule 5: After arrival at a new place, birds divide into producers and scroungers. The producers search food items or patches and randomly followed by scroungers to search food patches. Precise pseudo code of proposed method based on above defined rules of BSA is given in Algorithm 1. Three types of bird's behaviours such as vigilance behaviour, foraging behaviour, and flight behaviour are briefly described here.

A. FORAGING BEHAVIOUR
Stochastic decision (Rule 1) is taken according to the probability P of bird foraging food. Individual birds in swarms switched into foraging behaviours if probability P is greater than randomly selected constant value from a uniform normal distribution (0,1), otherwise the bird has vigilance behaviour. The best recorded experience or memory about searching food items can be utilized to explore food patches, and social behaviour and information are shared immediately (Rule 2). It can be mathematically modeled as [56]: where, C represents cognitive and S represents social accelerated positive coefficients. p i,j represents the past best position (local optimal) of i th bird, while g j is the swarm shared best (global optimal) past position. The function rand(0,1) represents uniform distribution of numbers in (0, 1). The term x t i,j (i ∈ {1, 2, ..., N }) represents N virtual birds' position at time t, having vigilance or foraging behaviour and term (j ∈ {1, 2, ..., D}) represents dimensions of available search space in which birds take flight.

B. VIGILANCE BEHAVIOUR
According to Rule 3, birds would not travel directly towards the swarm's center. However, individual birds may struggle to travel towards the swarm's center and birds' movement may be affected by competition with each other. Individual bird's movement or vigilance behavior modeled as [56]: where, positive constants a1 and a2 are within range [0,2]. mean j is the average position of j th bird's swarm and k(k ̸ = i) represents a random nonnegative integer (k ∈ {1, 2, ..., N }). The sumF it is best fitness values sum of birds swarm and pF it i is i th bird best fitness value. The smallest positive constant ε is used for avoiding zero division error. VOLUME 4, 2021 Individual birds travel towards the swarm's center because of indirect and direct effects. The swarm average fitness value is measured in the form of indirect effect and induced by environments. The direct effect is made by specific interference and A2 is used to simulate it. If k th bird (k ̸ = i) best fitness value is better as compared to i th bird best fitness value, in that case A2 > a2. It indicates the k th bird may suffer a smaller interference as compared to i th bird. Therefore, k th bird would quickly move towards the swarm center as compared to i th bird. However, there are some random movements of birds towards the center of the swarm. In the case of minimizing optimization problems, the bird's smallest fitness value is considered a better fitness value.

C. FLIGHT BEHAVIOUR
The birds in a swarm have flight behaviour due to foraging behaviour or any other reason. During the flight behaviour the birds in the swarm can be often switched again into two types of birds; 1) producer birds and scrounger birds on the basis of their food reserves (Rule4). The producer behaviours and scrounger behaviour can be written as, respectively [56]: where, k(k ̸ = i), represents a nonnegative integer (k ∈ {1, 2, ..., N }). FL is the following coefficient for hunting food -bird's scrounger behaviour (Rule 5). The term randn(0, 1) is random number normal distribution.  performance evaluation and its further detailed is available in Table 9.
We have implemented the proposed method and other algorithms in MATLAB R2017a and used the MAT-POWER6.0 package for load flow calculation. The execution of simulation-based experiments have been performed on Microsoft Windows 10 64-bits with Intel Core(TM) i7-5500U CPU @2.40 GHz and RAM @8.00 GB. We have conducted two case studies for performance evaluation of the proposed BSA method. Initially, we have solved the OPF problems in a hybrid power system by considering the objective to reduce electricity generation cost. In second case, carbon tax (e.g., emission cost) is included in power generation cost minimization objective function to reduce emissions pollution.

A. CASE STUDY 1: GENERATION COST MINIMIZATION
The simulation-based experiments have been conducted for performance evaluation of the BSA-based proposed method and other metaheuristic algorithms to solve the OPF problems in a hybrid power system. The OPF problem objective to reduce the electricity output cost is written in Eq. 23, excluding last term (emission cost or carbon tax). In Eq. 23, the first term ξ T (P T G ) , represents the power output cost related to thermal energy sources. The second term N W F j=1 ξ w,j (P w,j ) is power output cost related to windfarms and the third term N P V k=1 ξ pv,k (P pv,k ) is power output cost of solar PV units.
The electricity generation cost of an individual windfarm is calculated by adding three types of wind related cost such as 1) wind-scheduled power cost, 2) reserve cost, and 3) penalty cost. Wind-scheduled power cost of an individual windfarm follows a direct relationship to wind-scheduled power and a high spinning reserve is required when windscheduled power is kept high. In such a case, overall wind power cost increases while at a lower rate penalty cost decreases. The wind speed and wind electricity generation are highly dependent on the value of scale parameter c of Weibull PDF and the lowest wind power cost achieved at an intermediate value. Similarly, the electricity generation cost of an individual solar PV unit is also calculated by adding penalty cost, reserve cost, and solar-scheduled power cost related to solar PV unit. Solar-scheduled power cost related to solar PV units also follows a direct relationship to solarscheduled power. It is observed that solar power output cost did not monotonically increase with the values of lognormal PDF parameters such as standard deviation σ and mean µ for solar irradiance. Therefore, serious attention is required to choose the suitable value of solar-scheduled power associated with solar PV units. If the mean µ value is kept low, then it is suggested to choose a smaller value for solar-scheduled power.
In Table 10, optimal values of objective function, parame-  Figure 8a. All load or PQ buses voltage magnitude profiles obtained during performance evaluation of BSA and other algorithms are plotted in Figure 9a. For case study 1, the trade-off between optimal power generation cost and convergence time of the BSA-based proposed method and other algorithms are presented in Table 11.

B. CASE STUDY 2: GENERATION PLUS EMISSION COST MINIMIZATION
In case study 2, the minimum electricity generation cost expressed in Eq. 23 including emission cost (e.g., carbon tax) is an objective function for performance evaluation of the BSA-based proposed method. In this study, a carbon tax rate C tax ($20/ton) is imposed on released N O x , SO x , and CO x from fossil fuel-based energy sources. It is observed that the generation of the clean form of energy from solar and wind have increased and emission (ton/h) pollution decreased in the power system by imposing a carbon tax.  h are obtained from DE, ABC, PSO, SHADE, and HSA, respectively. For case study 2, the convergence of optimal power generation costs using the proposed BSA method and other algorithms is plotted in Figure 8b. All load or PQ buses voltage magnitude profiles are plotted in Figure 9b. The trade-off between optimal power generation cost and convergence time of proposed method and other algorithms is represented in Table 12. The simulation-based experiment results indicate that power generation from RESs increased when the carbon tax (ton/h) was imposed on carbon emission from fossil fuel-based energy sources.

VI. CONCLUSION
We have proposed a new bio-inspired bird swarm algorithm for finding optimal solutions to the OPF problems in the traditional thermal, wind, and solar energy sources-based hybrid power system, in this study. In which, we have incorporated utility load demand uncertainty and stochastic nature power generation from RESs. The power generation cost for thermal energy sources is measured using a valvepoint loading effects quadratic fuel curve. The Gaussian PDF, Lognormal PDF, and Weibull PDF have been used for modeling uncertainty of utility load demand, stochastic solar irradiance, and wind speed, respectively. The simulationbased optimization results have shown the superiority of the BSA to solve the OPF problems by satisfying all constraints and minimum power generation cost 863.121 $/h is achieved in case study 1. It has been observed from optimization results that the generation of the clean form of energy from RESs has increased and emission pollution has decreased in the hybrid power system by imposing a carbon tax. In case study 2, the proposed BSA approach has also outperformed and minimum electricity cost 890.728 $/h is achieved as compared to other algorithms. The comparative evaluation and simulation-based optimization results confirmed the superiority of BSA approach over other metaheuristic algorithms. The optimization performance of BSA in terms of accuracy, stability, and efficiency have made it attractive for application to real-time optimization problems. The simulation results have encouraged for further study. In future, the application of proposed BSA approach can be extended in a large-scale traditional thermal energy sourcesbased power system to solve the other optimization problems such as unit commitment, chance-constrained OPF, Transient stability constrained OPF.

VII. ACKNOWLEDGEMENT
This work was supported by King Saud University, Riyadh, Saudi Arabia, through Researchers Supporting Project number RSP-2021/184. The work of author Ayman Radwan was supported by FCT / MEC through Programa Operacional Regional do Centro and by the European Union through the European Social Fund (ESF) under Investigador FCT Grant