Heuristic Algorithm Based Optimal Power Flow Model Incorporating Stochastic Renewable Energy Sources

Today’s electricity grid is rapidly evolving, with increased penetration of renewable energy sources (RES). Conventional Optimal Power Flow (OPF) has non-linear constraints that make it a highly non-linear, non-convex optimisation problem. This complex problem escalates further with the integration of RES, which are generally intermittent in nature. In this article, an optimal power flow model combines three types of energy resources, including conventional thermal power generators, solar photovoltaic generators (SPGs) and wind power generators (WPGs). Uncertain power outputs from SPGs and WPGs are forecasted with the help of lognormal and Weibull probability distribution functions, respectively. The over and underestimation output power of RES are considered in the objective function i.e. as a reserve and penalty cost, respectively. Furthermore, to reduce carbon emissions, a carbon tax is imposed while formulating the objective function. A grey wolf optimisation technique (GWO) is employed to achieve optimisation in modified IEEE-30 and IEEE-57 bus test systems to demonstrate its feasibility. Hence, novel contributions of this work include the new objective functions and associated framework for optimising generation cost while considering RES; and, secondly, computational efficiency is improved by the use of GWO to address the non-convex OPF problem. To investigate the effectiveness of the proposed GWO-based approach, it is compared in simulation to five other nature-inspired global optimisation algorithms and two well-established hybrid algorithms. For the simulation scenarios considered in this article, the GWO outperforms the other algorithms in terms of total cost minimisation and convergence time reduction.


I. INTRODUCTION A. MOTIVATION
Optimal power flow (OPF) in power systems research was first introduced in 1962 by Ebeed et al. [1] and since then multiple extensions and solutions have been proposed to solve OPF. It has particular significance when the distribution system operators (DSOs) aim to maintain reliable and The associate editor coordinating the review of this manuscript and approving it for publication was Jahangir Hossain . economical system operation in an electric power system. The main objectives of OPF are optimising generation cost, power loss minimisation, maintaining voltage stability and reducing greenhouse gas emissions, all while maintaining optimal settings of various system constraints. Special care must be taken to ensure that the constraints on power generator capabilities, the current carrying capacity of the line, the generator bus voltage and the power flow balance are all satisfied. During the process of optimisation, the optimal performance of the system is achieved when scheduled generator power, complex power flow in the lines and the voltage vector of buses are in accordance with the required operating state of the system. Traditional OPF involves just conventional fossil-fuel-fired base generation sources, and this already yields a highly non-linear, non-convex and mixed integer optimisation problem [2]- [4]. However, it is evident from recent literature in this area that the increased penetration of renewable energy sources (RES) has created new challenges at the planning and operational stage. The OPF problem escalates when the uncertainties associated with solar photovoltaic generators (SPGs) and wind power generators (WPGs) are considered, along with conventional power generation sources, when optimising generation cost.

B. LITERATURE REVIEW
In a typical power network, a large number of power plants are interconnected with the grid to supply energy. These generation sources have multiple constraints and multiple cost functions. Since the inception of OPF about half a century ago, numerous traditional optimisation techniques have been proposed for application in this field. These include non-linear programming, interior-point methods, quadratic programming and mixed-integer linear programming [5]- [7]. Because of their fast convergence and robustness in finding an optimal solution, some of these techniques have been successfully adopted by industry. However, a key issue with such optimisation methods is the requirement to first linearise the optimisation function. For this reason, the non-convex, non-smooth and non-differentiable properties of the optimisation function are often approximated.
To address this problem, heuristic optimisation algorithms have also been proposed. These aim to find an optimal solution for the power system without modifying the original cost function [8]. Generally, heuristic algorithms fall into two main categories, namely single solution based and population-based. Simulated annealing and tabu search best represent single solution based heuristic algorithms in this area [9], [10]. Population-based heuristic algorithms have been proposed, including genetic algorithms (GA), particle swarm optimisation (PSO), crow search algorithms (CSA), artificial bee colony (ABC), cuckoo search optimisation (CSO), differential evolution (DE) and success history-based adaptive differential evolution (SHADE) algorithms [11]- [17].
Various bespoke heuristic approaches have also been proposed in the context of OPF to enhance the efficiency of the search methods. The authors of [11] proposed a modified GA to solve the OPF problem by introducing an enhanced genetic operator for an improved problem-specific optimisation. When tested on the well-known IEEE RTS 96 and IEEE-30 bus system, the modified GA yields improved elitism and fitness scaling features, compared to a basic GA. Similarly, a quadratic cost function for OPF that considered multiple valve-point loading effects is more efficiently solved by an improved PSO algorithm (compared to the standard PSO approach) [18]. Karaboga and Akay proposed a population-based heuristic algorithm, ABC, that demonstrates competitiveness with other proposed methods for OPF (in part because of the robustness of the algorithm, but also possibly because fewer parameters were controlled).
With regard to modern heuristic techniques, it is essential to achieve a suitable balance between exploration and exploitation. When efficiently driven, the former emphasises the investigation capabilities of the algorithm in the search domain of unknown regions. The latter, by contrast, enhances the ability of the algorithm to find the global solution on the basis of the information provided by the exploration strategy. These two aspects are contrary to each other and thus remain a major challenge for the research community [20]. Note that standard ABC performs well in terms of the exploration process because of its randomness; however, the relatively poor exploitation phase can result in poor convergence [21].
Differential evolution (DE) has been suggested as a way to enhance ABC [17], [22]. This includes, for example, use of 'onlooker bees' with a predefined probability and knowledge of the current best solution. Gao et al. [22] exploited a chaotic system to, not only improve the initialisation phase, but also to efficiently modify the search mechanism to find an optimal solution. This approach depends on knowledge of the current best solution to improve the exploitation aspect of standard ABC. Tanabe and Fukunaga [17] proposed an advanced variant of DE, which yields an algorithm they call SHADE. Here, the settings of successful control parameters are used to guide the selection of future control parameters. This aims to ensure a suitable balance between the exploration and exploitation processes. Furthermore, a comparatively fast convergence rate is achieved for non-linear, multimodal and constrained optimisation problems.
The efficiency of SHADE is further increased when combined with an effective constraint handling technique, namely the superiority of feasible solution (SF) approach [23]. However, the resulting SHADE-SF [24] sometimes attains premature convergence (i.e. becomes trapped in a local solution) and the convergence rate can be very slow. Furthermore, the proposed techniques in [24], [25] are verified only on the IEEE-30 bus system which does not guarantee good performance over medium and higher bus systems (IEEE-57 and IEEE-118). For example, the authors of [24] ran 24000 iterations to determine their optimal solutions, requiring several minutes to converge to a local solution because of the high computational load. Similarly, the authors of [25], [26] ran their proposed techniques for fewer iteration numbers (i.e. 200), and did not get the constant convergence curves for other algorithms. For higher iteration numbers (i.e. 1000), there is a higher probability that other algorithms may outperform the proposed algorithms to find better solutions with less computational time. This implies that the algorithm's exploration and exploitation capabilities were not fully explored. This could be impractical for industrial power plant applications, which require fast and robust algorithms to handle uncertain demand. When an algorithm converges to a local solution in the search space, it might satisfy all the constraints but it could yield an inferior value for the objective function i.e. a much better solution may exist. In practice, this could mean spending more money to balance the same demand that could otherwise be achieved with algorithms that find the OPF global solution.
In the economic dispatch (ED) problem, system constraints (especially limitations on network parameters) may often have been ignored; however, complying with network constraints is essential for OPF. Reference [27] mentions system constraints, but does not explicitly address the question of how to satisfy these constraints. Furthermore, in the ED problem, emission aspects and voltage profiles are generally ignored; however, these are again all important in the case of OPF.
It is clear, therefore, that in a mixed network consisting of thermal power generators (TPGs), SPGs and WPGs, further research into OPF is required. A number of optimisation issues in this context are addressed in the present study, which has a particular focus on uncertainty modelling of SPGs and WPGs. The biggest challenge for incorporating SPGs and WPGs into the electricity grid is their intermittent nature. Normally, RES are owned by private operators, from which the grid DSOs sign an agreement for purchasing scheduled power. However, since electricity generation from these RES are uncertain, sometimes the power output may be more than the scheduled power, leading to underestimation of the available power level. The DSOs generally bear the penalty cost, since surplus power goes wasted if not utilised. By contrast, overestimation is the scenario in which the generated power is less than the scheduled power. To mitigate power demand, the DSOs must therefore keep spinning reserve power, adding to the ongoing operating cost of the system.

C. CONTRIBUTION AND PAPER ORGANIZATION
In the present article, a new objective function is formulated that considers direct, penalty and reserve costs of renewable sources, in addition to the generation costs of the TPG units [24]. Wind distribution is modelled using the Weibull probability density function (PDF) and solar irradiance is modelled with a lognormal PDF (see later section III for details and references). Generation cost is optimised and the effect on optimal scheduling of changes to reserve and penalty costs are investigated. Finally, with regard to emissions, fossil fuel driven TPGs emit harmful gases into the environment, while renewable sources do not. Hence, a carbon tax is imposed in some countries, usually in proportion to the emitted greenhouse gases [28]. For the relevant case study in the present article, a carbon tax is embedded into the objective function to investigate its effect on generator scheduling.
To summarise, the cost functions and associated optimisation algorithms that have been developed for OPF to date, have either relied on linear approximations or, when using modern heuristic techniques, have seen unsolved challenges in relation to the exploration and exploitation phases. The present research aims to address these limitations by means of grey wolf optimisation (GWO). First introduced by Mirjalili et al. [29], GWO has been proven to be flexible, easy to apply, scalable and, most importantly in the present context, has an inherent capability to strike a practically useful balance between exploration and exploitation [30] (see section IV).
To the present authors' knowledge, the application of GWO to OPF in the presence of uncertain power outputs from RES has not yet been documented in the literature. In this article, we propose this method to handle the OPF problem, which is to be our main research contribution. The performance of the new approach is evaluated and benchmarked against other well recognised evolutionary algorithms such as GA, PSO, CSA, SHADE-SF, ABC and two well-established hybrid algorithms, namely GA-PSO and ABC-CSO, for the modified IEEE-30 and IEEE-57 bus test systems. We investigate the potential for GWO to reduce both the total fuel cost and optimisation convergence rates. We evaluate ten scenarios, for example involving carbon tax and for different types of renewables. Hence, novel contributions are made in three main areas: the new objective functions for OPF; the use of the GWO approach to optimise objective functions both in small and medium-scale systems; and a simulation based investigation for selected case study examples to demonstrate the benefit of the proposed approach in terms of operation cost, computational time and scalability.
The rest of the article is organised as follows. Section II describes the mathematical model and associated constraints for OPF. Section III presents the uncertain SPG and WPG output models. Section IV develops the new approach for applying GWO to OPF with the presence of uncertain RES. The six algorithms under consideration are compared for several realistic case study problems in simulation in sections V (IEEE-30 bus), VI (IEEE-57 bus) and VII (hybrid algorithms). Finally, conclusions are presented in section VIII. Tables 1 and 2 give details of the main nomenclature and acronyms used throughout this article, whilst important parameters for the IEEE-30 and IEEE-57 bus system networks are summarised in Appendix Tables 3 and 8, respectively. Focusing initially on the IEEE-30 bus system for brevity and illustrative purposes, Fig. 1 shows that the network comprises three different power generation sources i.e. TPGs, WPGs and one SPG. The main characteristics of this model are almost similar to the one proposed in [24]. Outputs of the SPG and WPGs contain variations which need to be balanced with the help of reserve and other generator outputs collectively. Thus, the total generation cost includes the total operational cost of the TPGs, together with the penalty and reserve costs due to intermittency in the SPG and WPG outputs. Penalty and reserve cost details are provided in subsequent sections.

II. MATHEMATICAL MODEL
Although the present article focuses on OPF, power systems are in general, highly nonlinear and, therefore, many control and optimisation issues become hard-to-solve  problems unless linearised. As a result, GWO offers the potential to improve other power system problems, such as controlling flexible alternating current transmission system (FACTS) devices and the optimisation of the placement of distributed generators.

A. COST MODEL FOR THERMAL POWER GENERATORS
TPGs are operated on fossil fuel. The relationship between generator output power (MW) and fuel cost ($/hr) is straightforwardly expressed with the following quadratic equation, where a i , b i and c i are cost coefficients associated with the i-th TPG, while NT represents the total number of TPGs. However, cost function modelling with a valve point loading effect has a precise and more realistic impact on the quadratic cost function. In practice, TPGs operation is based on controlling the steam valves for turbine operation through distinct nozzles. When the individual nozzle of the multi-valve system operates at its full output, the highest efficiency of a TPG is achieved [31]. Fig. 2 illustrates a multi-valve loading effect on the quadratic cost function. To model these valve point loading effects, the basic cost function in Eq. 1 is altered with the addition of the absolute value of the sinusoidal function for a multi-valve steam turbine. Consequently, the TPG units total cost ($/h) becomes, Valve point loading effect on a quadratic cost function. VOLUME 8, 2020 In Eq. 2, g i and h i are cost coefficients representing the valve-point loading effect, while P min Tg,i denotes the minimum power which the i-th TPG produces when it is in operation. The coefficients used in this work are provided in Appendix Table 4.

B. DIRECT COST OF WIND AND SOLAR PHOTOVOLTAIC POWER
TPGs are fossil fuel-fired. When compared to such conventional generators, SPGs and WPGs need no fuel for operation. In this case, when RES belong to DSOs, only the initial outlay or maintenance cost of the RES are assigned [32], [33]. However, when the ownership of RES belongs to private parties, scheduled power obtained from RES is charged in accordance with the mutually agreed contract.
The direct cost of the j-th wind power plant in terms of scheduled power is modelled as follows, where d w,j and P Ws,j represent the direct cost coefficient and scheduled wind power associated with the j-th WPG, respectively. Similarly, the direct cost of the k-th solar power plant is determined using, where d s,k and P Ss,k are the direct cost coefficient and scheduled solar power from the k-th SPG, respectively. Although there is only one SPG in the present case study example ( Fig. 1), the mathematical formulation utilises a k subscript here so as to develop and solve the generalised problem.

C. COST EVALUATION OF UNCERTAINTIES IN WIND POWER
Due to the intermittent nature of RES, two situations may be encountered with respect to the energy generation profile. Situation one arises when the generated power from RES is less than the expected value. This is referred to as overestimated output power. To compensate for overestimated power and to provide uninterrupted power supply to end consumers, the spinning reserve needs to be maintained by system operators on the generation side. The cost associated with reserve generating units, as required to address the overestimation problem, is termed the reserve generation cost [34]. The cost for the j-th wind power plant is determined using, where r w,j is referred to as the reserve cost coefficient pertaining to the j-th wind power plant, P Wa,j is the available power from the same wind power plant and f w (P W ,j ) is the probability density function for the wind power of the j-th power plant. The output power probability calculation from various WPGs at different wind speeds is discussed later (section III-B). The second situation arises when the the generated power from RES is greater than the estimated power. The surplus power is potentially wasted. In this case, the DSO aims to reduce output power from traditional TPGs. This situation is referred as underestimated output power from wind energy resources, and the associated cost is called the penalty cost. This penalty cost, paid by the DSO in proportion to the additional power generated from WPG, is determined as follows, where p w,j and P Wr,j represent the penalty cost coefficient and rated output power from the j-th WPG.

D. COST EVALUATION OF UNCERTAINTIES IN SOLAR PHOTOVOLTAIC POWER
The output power from SPGs in the network is also intermittent and uncertain in nature. The method to solve over and underestimation of SPG output is similar to that used for the WPGs. However, one distinct difference is that solar radiation follows the lognormal PDF as compared to the Weibull PDF for the wind distribution [35]. In the present work, the reserve and penalty cost models are built following similar concepts to those proposed by [36].
In the case when the generated output power is less than expected, the reserve cost is calculated as follows, where r s,k and P Sa,k represent the reserve cost coefficient and available power respectively, associated with the k-th SPG. The solar power shortage probability is represented by f s (P Sa,k < P Ss,k ), while E(P Sa,k < P Ss,k ) defines the expected power of the SPG below P Ss,k . The penalty cost for the k-th solar power plant is, where p s,k represents the penalty cost and f s (P Sa,k > P Ss,k ) shows the probability of surplus power generated by the k-th solar power plant as compared to P Ss,k . Finally, E(P Sa,k > P Ss,k ) is the expected surplus output power.

E. CARBON TAX BASED EMISSION MODEL
Traditional TPGs release greenhouse gasses into the environment. The emission of harmful gasses such as NO x and SO x into the environment increases when the generation from thermal power generators increases (in p.u. MW). This direct relationship is represented by Eq. 9 i.e. harmful emissions in tons per hour (ton/hr), where α i , β i , γ i , ω i and µ i are the emission coefficients of the i-th TPG. Appendix Table 4 shows the emission coefficients used in this research. These values are similar to those introduced by [37], except for a small adjustment to coefficient µ for the generator connected with bus-1.
In recent years, to produce clean energy and to protect the environment from the effects of harmful gases, notably to address global warming, many countries are imposing a carbon tax on greenhouse gasses emissions [28], [38], [39]. Due to the associated additional cost, the energy production sector is under enormous pressure to reduce such emissions or to produce a cleaner form of energy from RES. In the present work, a carbon tax is optionally imposed on the level of greenhouse gas emissions. The carbon emission cost ($/hr) is determined as follows, where C em and C t represent the carbon emission cost and carbon tax per unit amount of greenhouse gasses, respectively.

F. OBJECTIVE FUNCTIONS
The objective functions for OPF are formulated from the various model components discussed above. In this article, Objective 1 is based on the sum of the costs from Eqs. 2-8, whilst Objective 2 also includes the emissions from Eq. 9. Hence: Objective 1: Minimise, where NW and NS are the number of WPGs and SPGs in the network, respectively. To study the impact of carbon tax on generation scheduling, the second objective function is constructed by adding the emission cost to Eq. 11. Objective 2: Minimise, Both OPF objective functions, Eqs. 11 and 12, are subject to system equality and inequality constraints, as discussed below.

1) EQUALITY CONSTRAINTS
Equality constraints represent typical load flow equations in a power system. These constraints are used for power balancing of both real and reactive powers generated to the total demand and loss in a system. The equality constraints are stated below [24], (a) Active power constraints: (b) Reactive power constraints: In Eqs. 13-14, σ ij = (σ i − σ j ) represents the voltage angles difference between bus-i and bus-j and NB represents the total number of buses. The active and reactive power demand at bus-i is represented by P Di and Q Di whilst the active and reactive power generation is represented by P Gi and Q Gi , respectively. Power generation can either be from conventional power generators or through RES. The transfer conductance and susceptance between bus-i and bus-j are represented by A ij and B ij , respectively.

2) INEQUALITY CONSTRAINTS
The inequality constraints define operating limits for the equipment and components in the power system. These also relate to the security constraints on load buses and lines.
(a) Generator constraints: V min (b) Security constraints: where S max lq in Eq. 23 and the similar terms in Eqs. 15-22 represent the constraint limits. In particular, Eqs. 15-17 define active power limits on the TPGs, SPGs and WPGs while for the same generators, Eqs. 18-20 define reactive power capabilities. Furthermore, Eqs. 21-22 apply voltage limits on the generator buses and load buses (PQ). NG and NL represent the number of generator buses and load buses, respectively. Finally, line flow limits on apparent power oscillations VOLUME 8, 2020 are defined by Eq. 23 for the total number of transmission lines (Nl).
It is pertinent to mention here that, after achieving an optimised solution for power flow, the equality constraints are satisfied automatically via the power balance equations. By contrast, the inequality constraints are control variables. These include the generator active power and generator bus voltages, and are intended to be self-limiting. When the optimisation algorithm is applied to choose a feasible solution, the selected value of these control variables lie in the bounded range. However, inequality constraints associated with the slack bus generator, the reactive power of other generators, line capacities and voltage limits on load buses, all require special attention. Hence, section IV of the article describes the handling of inequality constraints with control variables in more detail.

G. LOAD BUS MODELLING
In OPF studies, generator reactive power capability has an important role. With regard to TPGs, narrower implementation ranges are defined in this study, compared to the ranges defined by e.g. [32], [39], [40]. This is because, in recent years, the reactive power capabilities have evolved. Wind turbine reactive power profiles and other relevant features are now commercially available [41]. With the help of the Enercon FACTS wind turbine capability curve, it is clear that a wind turbine can deliver reactive power from 0.4 p.u. to 0.5 p.u. during its active power output range. The reactive power absorbing capability of the generator can be enhanced with the help of negative reactive power delivery.
A rooftop SPG can be modelled as load bus (PQ) with Q = 0. However, large photovoltaic generation facilities are equipped with converters. Considering the dynamic behaviour of converters, it is desirable to conduct full-scale generator modelling for P-Q capability [42]. Some articles in the literature consider controller and converter models when conducting a detailed analysis of SPG reactive power capabilities [43]. In [44], the authors extended their study to analyse the impact of variation in radiation and ambient temperature on photovoltaic capability. In the present study, the generator active (P) and reactive (Q) power parameters are set according to Appendix Table 6, whereas the reactive power capability of SPG is set between 0.4-0.5 p.u..
In the OPF problem, system parameters like real power loss in the network and voltage deviation are also important. Some of the power losses in the transmission system are unavoidable because of the inherent resistance in transmission lines. Hence, the network losses are determined as follows, where A ij is the transfer conductance and σ ij = (σ i − σ j ) is the voltage angles difference between bus i and bus j.
In a power network, voltage deviation indicates the relative voltage quality in the system. To formulate a voltage deviation indicator in the network, a nominal value (i.e. 1 p.u.) is taken as a reference value for cumulative voltage deviation for all load buses. This is expressed as follows,

III. STOCHASTIC SOLAR POWER, WIND POWER AND UNCERTAINTY MODELS
It is well-known that PDFs can be used for mean power calculations of wind turbines [34], [36], [45], [46]. The wind speed (v) m/s follows a Weibull PDF and is determined using a scale parameter (c) and shape parameter (k) as follows, The Weibull distribution mean is defined, To compute the gamma (x) function, In the modified IEEE-30 bus system model, conventional TPGs at bus-5 and bus-13 are replaced with WPGs. In the proposed case studies, we use these PDF parameters to compute wind speed. Figs. 3 and 4 show the wind frequency distribution based on Weibull fitting [24]. The output curve is achieved after running 8000 Monte-Carlo scenarios. Wind turbine design requirements are specified in [47], i.e. the highest turbulent class IA of turbine and maximum average speed of 10 m/s at hub height. For the simulations reported below, k and c are carefully chosen to ensure both diversity and realistic geographic locations for wind farm sites, with their values given by Appendix Table 5. For both wind farms, the value of shape parameter is 2 which corresponds to the moderately gusty winds. In Northern Europe and most other locations around the world, this value for the shape parameter is often assumed [48]. Also, to gain the capacity factor as for a real wind farm (30-45%), the values of scale parameter for both wind farms are assumed to be c = 9 and c = 10 [49]. Similarly, the IEEE-30 bus system bus-11 conventional generator is replaced with the SPG unit. In [35], the author investigated frequency distributions of global radiations at important metrological stations in Taiwan. According to the study, the lognormal function describes the frequency distribution quite better where weather conditions are more dispersive. In this study, the parameters for lognormal distribution are determined using the corresponding mean and standard deviation of the global irradiation in Taiwan. The output of the SPG has a direct relation with solar irradiance (I ) which follows lognormal PDF. The solar irradiance probability is dependent on the standard deviation (λ) and mean (ψ) when it follows a lognormal PDF, as follows,   The lognormal distribution mean is defined, simulations with a sample size of 8000. Values for selected parameters of the lognormal PDF are assumed using the corresponding mean and standard deviation of the global irradiation in [35] and summarised in Appendix Table 5. We use these values in the simulation study, with the exception of section V-F, in which they are modified in order to observe the impact of parameter variation on the total cost. Wind turbine output is expressed [24], where v in , v o and v r represent the cut-in, cut-out and rated output wind speed of the wind turbine, respectively, while P Wr defines the rated output power of the wind turbine. According to the product data sheet of Enercon E82-E4 where I sr represents solar irradiance in a rated environment i.e. 800 W /m 2 , I c represents a specific irradiance point, here 20 W /m 2 and P Sr is rated output from the SPG.

B. WIND POWER PROBABILITY MODEL
With reference to Eq. 31, it may be observed that WPG output is categorised into three distinct features. This is due to the fact that wind speed is not constant in all regions. The wind turbine output power is zero when it encounters a wind speed (v) which is below cut-in speed (v in ) or above cut-out speed (v o ). The wind turbine produces rated output P Wr when it encounters the rated wind speed (v r ) or below cut-out speed (v o ). Hence, the output of the wind turbine for the first and second eventuality in Eq. 31 for being 0 or P Wr is determined as follows [51], The output power of a wind generator is continuous between cut-in and rated speed of wind. Hence, the probability for the continuous region is determined as follows [24], C. SOLAR POWER OVER/UNDERESTIMATION COST As observed from the histogram in Fig. 6, the SPG unit has stochastic output power because of the variance in solar irradiance. The dotted line shows the scheduled output power needed to supply to the grid. It is important to note that scheduled power has no fixed value, rather there is a mutually agreed power level between the DSO and the private party which sells solar power. For the calculation of under and overestimation costs of the SPG, the following equations are used in the model [24], where P Ss+ and P Ss− represent the surplus power and shortage power, as lying on the left and right half plane of schedule power (P Ss ) in the histogram of Fig. 6. Similarly, f ps+ and f ps− are relative frequencies for the occurrence of P Ss+ and P Ss− . N + and N − represent number of discrete bins on the right and left planes of P Ss for PDF generation.

IV. OPTIMISATION TECHNIQUE
GWO was proposed by Mirjalili et al. 2014 [29] and, in a relatively short period of time, has already attracted significant research interest. It is inspired by the leadership and hunting behaviour of grey wolves which live in the form of a pack. It has been widely used for different optimisation problems and can show improved characteristics over other swarm intelligence techniques: its initial search is based on relatively few parameters for which no initial derivation is required. Furthermore, the approach is flexible, straightforward to apply, scalable and most importantly for the present work, it helps to strike an accurate balance between exploration and exploitation.
In the real world, grey wolves adopt a social hierarchy that has been categorised into four different levels: alpha (α), beta (β), delta (δ) and omega (ω) wolves. From the top to bottom of the leadership hierarchy, α wolves are known to be the superior. Their role is decision making in the pack. Alpha wolves are followed by β wolves, whose role is to help α wolves in decision making and to carry out other important activities in the pack. At the bottom of this hierarchy come the δ and ω wolves. Omega wolves are also known to be the scapegoats. They are subordinates to all other wolves.
In regard to GWO, accurate determination of prey location is treated as the optimisation problem (fittest solution), while the position of the wolves relative to the prey determines the best solution. The position of the α wolves is said to be the best solution found so far in the search space, because they are expected to be closer to the prey than other wolves in the pack. Similarly, β and δ wolves determine the second and third best solutions in the search space because of the hierarchical classification and corresponding position towards the prey/optimal solution. To allocate their position in the search space, these wolves are represented as X α , X β and X δ . Fourth level ω wolves update their position X ω in accordance with the relative position of the α, β and δ wolves. Finally, hunting for prey is achieved by adopting four main steps, namely encircling, hunting, attacking and searching again. These steps are represented in the GWO as follows.

A. PREY ENCIRCLING
GWO starts with a step that is analogous to chasing and encircling the prey. To mathematically model the encircling behaviour of grey wolves corresponding to the prey location, the following equations have been proposed, where, and t indicates the current iteration, while − → X (t) and − → X p (t) are position vectors representing the current location of the grey wolf and prey in the search space, respectively. The coefficient vectors − → A and − → C are determined as follows, To control exploration and exploitation, the components of − → a are linearly decreased from 2 to 0 over the course of an iteration. Note that − → r 1 and − → r 2 are random vectors whose values are chosen between [0, 1]. To reach prey position (X p , Y p ), the current position of a grey wolf (X , Y ) is updated with Eqs. 38-41. The value of − → a is assumed the same for all the wolves in a population. A wolf can update its position according to the best agent in different places by setting the values of − → C and − → A .

B. HUNTING
After finding the prey location, the grey wolves encircle it. The α wolves guide the pack for prey hunting, while β and δ wolves also contribute. Initially, the α, β and δ wolves location are saved as the 'best' location, representing their better knowledge to recognise prey location. The remaining search agents, mainly ω wolves, update their location in accordance with the position of the best search agents. For α, β and δ wolves, position location is determined as follows, At iteration t, the distance between − → X (t) and the three best hunt agents ( − → X α ), ( − → X β ) are ( − → X δ ) are determined using Eqs. 42-47, in which A 1 , A 2 and A 3 are random vectors as defined in Eq. 40. Finally, wolves movement towards prey is updated by Eq. 48.

C. ATTACKING THE PREY (EXPLOITATION)
Hunting ends when grey wolves attack the prey. It is possible when the prey stops moving around, and grey wolves start exploiting its position. Mathematically, prey approaching behaviour of grey wolves is modelled when the value of the exploration rate − → a is decremented from 2 to 0 over the course of an iteration. The optimum location of a prey is represented as 0 in the search space. Note that, the fluctuation range − → A in Eq. 40 decreases with the decremented value of − → a . This is due to the fact that − → A chooses a random value between [−2a, 2a] and the value of − → a is decremented with every iteration to locate prey for attacking. Exploitation is emphasised when agents attack prey whilst the value of − → A lies between [−1, 1]. This shows that the agent is ready to carry out an attack, since they are one step behind, between the prey position and their current position in the search space.
The parameter − → a is linearly updated as follows, where T indicates the maximum iteration number, set to 1000 in this study.

D. SEARCHING AGAIN FOR PREY (EXPLORATION)
The α, β and δ wolves position in the search space guide the whole pack to search for prey. Initially, all wolves diverge from each other to first locate the prey, before subsequently converging to attack the identified prey. This behaviour of divergence and convergence is obtained when the value of − → A is randomly chosen between −1 > A > 1. The algorithm tries to search the global candidate solutions when the value of | − → A | > 1 forces grey wolves to diverge from the prey. Similarly, the value of | − → A | < 1 helps grey wolves to converge towards the fitter prey. These random values of − → A enhance the search space area for wolves and obligate them to diverge into a comparatively larger area to search for prey. This facilitates the GWO algorithm to search globally for an optimum solution. Due to these exploration and divergence characteristics, the GWO algorithm can find fitter prey than other approaches.
During a new search, if a wolf finds a better prey closer to it, that wolf becomes an α and, based on distance, other wolves are divided into β and δ wolves. Here, − → C is another important component which emphasises the exploration process in the GWO algorithm. Random weights containing values between [0, 2] are assigned to each prey in the search space via Eq. 41. If C > 1, the prey needs to be emphasised, while if C < 1, that prey is de-emphasised. The value of the − → C component helps to avoid a local optimal solution, and helps the GWO algorithm in general to adopt more random behaviour once the optimisation process has started. In contrast to the values of A, the value of C is not decremented linearly. This is part of the deliberately engineered behaviour of the GWO algorithm i.e. to provide random values not only during the initial iterations but right through to the final iteration, in order to maintain good exploration. With the help of this component, local optimum stagnation is avoided, not only during the initial iterations but also in the final iterations when local optimum stagnation is otherwise frequent.
In nature, many obstacles are faced by wolves before attacking prey. Due to this fact, a rapid approach to the prey is essential. This behaviour in the GWO algorithm is achieved with the help of the C vector. When the algorithm assigns random weight values C for the prey, it becomes harder and it is further to go for grey wolves to approach the prey and vice versa. Finally, the GWO algorithm terminates when the end criteria are met. VOLUME 8, 2020

V. CASE STUDIES AND RESULTS FOR IEEE-30 BUS SYSTEM
In this section, we verify the effectiveness of our proposed optimisation framework and the chosen GWO algorithm, using the modified IEEE-30 bus system introduced earlier.
Numerical results using the GWO approach are compared with those obtained by GA [11], PSO [12], CSA [13], ABC [14] and SHADE-SF [24]. To perform the simulation work, a core i7 Mac book processor with 16 GB RAM is used. The MATPOWER packages proposed in [41] are used for the power flow calculations. We present six case studies. Case-1 is a benchmark simulation to optimise total generation cost. In Case-2, total generation cost is optimised when a carbon tax is imposed on emissions from conventional TPGs. Case-3 schedules power generation sources while considering stochastic WPG and SPG underestimation and overestimation costs. In Case-4 and Case-5, power generation costs are optimised while considering reserve and penalty costs. Finally, Case-6 describes how Weibull and lognormal variable variations affect the WPG and SPG capabilities. In a single run of the algorithm, a maximum of 1000 iterations are performed as the end criteria.

A. CASE-1: MINIMISING TOTAL GENERATION COST
By making use of Eq. 11, Case-1 performs optimisation scheduling of both TPGs and RES to minimise total generation cost. The direct cost coefficients of wind power are d w,1 = 1.6 and d w,2 = 1.75. The penalty cost coefficient for not fully utilising wind power is assumed as p w,1 = p w,2 = 1.5 and the reserve cost coefficient for overestimation is r w,1 = r w,2 = 3. These values are used for illustrative purposes. Finally, the PDF parameters for the WPGs and SPG are given in Appendix Table 5, with many of these parameters being almost same as in [24]. With these settings, Fig. 7 compares the convergence characteristics of different optimisation techniques. For this case study, Appendix Table 6 summarises the optimum results for all the control variables, such as total generation cost, reactive power (Q) and other important parameters. A voltage at the i-th bus is denoted by v i . Similarly, with the help of Eqs. 24-25, power loss (P l ) and voltage deviation (V d ) are determined. Note that P Wg,1 and P Wg,2 signify scheduled power from the two wind generation sources. From simulation results, it is found that that the GWO and SHADE-SF algorithms are more efficient with fast convergence and better solution quality when compared to the other well established algorithms for similar OPF frameworks. The minimum generation cost achieved by SHADE-SF and GWO are 782.30 and 781.40, respectively. Hence, for this scenario, GWO outperforms SHADE-SF and all other algorithms in terms of the total cost and elapsed time.

B. CASE-2: MINIMISING TOTAL GENERATION COST WHEN CARBON EMISSION TAX IS IMPOSED
A carbon tax (C t ) with a rate of $20/ton is imposed in this case study [52]. The objective is to minimise the cumulative cost by utilising Eq. 12. With the imposition of the carbon tax, penetration of RES is expected to increase, and this is evidenced by the simulation results. Appendix Table 7 provides the optimum power generation schedule of all relevant parameters, including total generation cost (with the carbon tax), reactive power of the generators and other important parameters required for OPF. In Case-2, a higher penetration of RES is achieved as compared to Case-1, when no penalty was imposed on carbon emissions. The extent of RES penetration in the optimum generation schedule depends solely on the emission volume and rate of carbon tax imposed. For this scenario, Fig. 8 compares the convergence characteristics of GWO and other techniques to reveal that GWO has the best performance in terms of total cost minimisation. An important factor that needs to be critically addressed in OPF problems is the load bus voltage. Operating voltages for all buses need to be within the range 0.95-1.05 p.u.. In this regard, Fig. 9 illustrates the voltage profiles for both Case-1 and Case-2. These results show that the requirements are satisfactorily met after optimisation. The remaining case  studies (Case-3 to Case-6) yield similar voltage profiles and so, for brevity, the equivalent plots are omitted.

C. CASE-3: SCHEDULED POWER VS THE COST OF WIND AND SOLAR POWER GENERATORS
Appendix Table 5 shows the Weibull PDF parameters used for the analysis in this case study, while section III-A discussed wind turbine parameter selection. The cost coefficients selected for this case study are similar to case-1. Note that the average cost of the TPGs is higher than the direct cost of RES. Similarly, the direct costs are higher when compared to the penalty cost for not fully utilising available wind power [52]. In these simulations, scheduled available power for the two wind farms is varied from zero (0) to the rated power, as plotted in Figs. 10 and 11. Total costs represent the sum of direct, reserve and penalty cost of the corresponding scheduled power. There exists a linear relationship between direct cost and scheduled power. When the scheduled power from RES increases, larger spinning reserves are required, which increases the reserve cost and consequently generation costs move upwards. Contrary to the reserve cost, penalty cost decreases at a lower rate with increased scheduled power from RES.
Similarly, when SPG output is over/underestimated from scheduled power, cost variations occur because of the associated penalty and reserve costs coefficients. Fig. 12 illustrates the change in solar power generation cost for scheduled power. To evaluate the total cost of SPG, operation and maintenance cost needs to be analysed. It is learnt from [53] that the cost ranges selected for this study are similar to those of onshore wind power plants. Therefore, in this study, the direct, reserve and penalty cost coefficients are assumed as d s = 1.6, r s = 3 and p s = 1.5, which are almost similar to [24]. It is important to note that the total cost of solar power generation does not increase uniformly with specified PDF parameters of solar irradiance. Cost plots for this case study show that when scheduled power from SPG is set to 20 MW, the minimum cost is achieved.

D. CASE-4: EFFECT OF PROBABILITY DENSITY FUNCTIONS ON WIND AND SOLAR POWER GENERATOR COST
The Weibull distribution scale parameter (c) has direct impact on WPG cost. This case study evaluates how, for a fixed arbitrarily scheduled power, WPG cost changes with the variations in scale parameter (c) of Weibull distribution. The value of the shape parameter for both wind farms is k = 2 because, the Rayleigh distribution is equivalent to a Weibull distribution with k = 2, corresponding to moderately gusty winds. The values for the cost coefficients are identical to those used in Case-1. The outputs of P Wg,1 and P Wg,2 are 25 MW and 20 MW respectively, which is one third of the total installed capacity. This assumption appears realistic, since real wind farms contribute a capacity factor between 30%-45% [49]. The relationships between P Wg,1 and P Wg,2 costs and the Weibull distribution scale parameter are shown in Figs. 13 and 14. Minimum costs are achieved when the value of the scale parameter is in the middle of the specified range. A higher valued scale parameter implies the prevalence of higher wind speeds with a certain probability. This is due to the fact that for a fixed interval, scheduled power remains the same, which increases penalty costs and hence the overall costs also rise. However, above a certain value of the scale parameter, the reserve cost reductions are insignificant.    earlier in Case-1. The minimum solar power cost is achieved when λ = 5.5. Also, when λ = 5.8, the reserve cost and penalty cost values are the same. For higher values of λ, penalty costs increase sharply, pushing the overall cost to a higher level. Finally, note that SPG output and solar irradiance have a direct relationship with λ. When λ is lower, then the output of the SPG is also low. To withstand this situation, high reserve powers are required, which increases the reserve cost. When λ is relatively high, high solar irradiance is expected, hence increasing output power from the SPG. In such a scenario, the penalty costs yield an increase in the overall cost. Keeping in mind these two scenarios, an appropriate value from the SPG always needs to be scheduled.

E. CASE-5: OPTIMIZED COST VS RESERVE COST
Case-5 evaluates the relationship between optimized cost and reserve cost. Values for the solar and wind power generation reserve cost coefficients are increased from r w,1 = r w,2 = r s = r = 4 to r = 6 with an increment of 1. Penalty cost coefficients for RES are similar to those used in Case-1 and Case-2, with p = 1.5. For these parameters, the optimized schedules for all the generators are illustrated in Fig. 16. Reserve costs are varied and three different cases are considered i.e. r = 4 (Case-5a), r = 5 (Case-5b) and r = 6 (Case-5c). As shown in Fig. 17, increasing the reserve cost yields an inverse relationship with RES optimum power scheduling. This is because increased reserve cost coefficients imply higher costs for spinning reserve and a reduction in the RES contribution to optimum power scheduling. This gap is compensated for with increased output from the TPGs. In conclusion, with an increase in the reserve cost coefficient values, the contribution of both the SPG and WPG decreases, yielding an increase in the overall generation cost.
F. CASE-6: OPTIMIZED COST VS PENALTY COST For these simulations, most of the parameters are identical to those considered in Case-1, except for the penalty costs associated with the SPG and two WPGs. These are  varied from p = 1.5 to p = 5 in discrete steps of 1. Here, p w,1 = p w,2 =p s = p = 1.5 is increased to p = 3 (Case-6a), p = 4 (Case-6b) and p = 5 (Case-6c), similar to [24]. The reserve cost coefficient values for RES remains unchanged from those used in case-1 and case-2 i.e. r = 3. The optimised schedules for all the generators output are shown in Fig. 18.
The penalty cost is imposed when power generation from RES is higher than the expected power. In such a scenario, with relatively high penalty costs, there is a need to raise the scheduled power from RES if solar irradiance and wind speeds are high. The strategy to increase scheduled power from RES helps to reduce the penalty cost. In Case-5, when the reserve cost increases, the outputs from RES monotonically decrease. However, in Case-6, when the value of p is increased, the output from solar generation will occasionally be decreased. This is due to the highly non-linear relationship of wind and solar power reserve/penalty costs with the probability distribution of these sources. Fig. 19 illustrates the four different costs, i.e. TPG, SPG, WPG and total cost. Total cost is slightly increased due to the small upward fluctuations of solar and wind power generation. TPG costs, however, remain steady for different  values of p. The relationship between penalty cost, reserve cost and the voltage range is illustrated by Fig. 20. Here, the bus voltage range is specified as 0.95-1.10 p.u. and the reserve/penalty cost cases are combined to analyse the overall impact. For all these cases, the voltage profile of different buses is ideal because the voltage value lies within the specified limit. However, for bus-8, the generator voltage  shows a significant change for different values of the cost coefficients. This is because P Tg,3 is connected with bus-8, hence variations in the reactive power output when different cases are simulated, can yield abrupt changes in the output of bus-8.
To solve the OPF problem, the operating limits of the power system states or dependent variables should be satisfied. The state or dependent variables include the load bus voltages, the generator reactive powers and the line flows. The active power loss, voltage profile and voltage security in a power system strongly depend upon the flow of reactive power in the transmission lines [54]. Appendix Tables 6 and 7 specify limits on reactive power and scheduled reactive power profiles for all the generators, as also illustrated in Fig. 21. For all of these cases, the reactive power of all the generators successfully lies within the required limits. For the optimisation problems considered in this section, a minimum violation of reactive power constraints is desired. One advantage of the GWO approach, is that it allows for network operation components to lie close to the defined limits i.e. GWO provides an efficient method for handling the non-linear constraints aspects of the OPF problem.

VI. CASE STUDIES AND RESULTS FOR IEEE-57 BUS SYSTEM
To confirm the robustness and scalability of the proposed GWO algorithm, the modified IEEE-57 bus test system, representing a medium-scale power system, is investigated. The active and reactive power demands of this system are 1250.8 MW and 336.4 MVAR, respectively, at 100 MVA base. More details about the system are given in [38]. To execute the optimisation process for all algorithms, the population size is 50 and the maximum number of iterations is set to 1000. The main characteristics of the IEEE 57 bus system are summarised in Appendix Table 8 whilst the cost and emission coefficients are described in Appendix Table 9. For the following case studies in the IEEE-57 bus system, optimal solutions are obtained under similar equality and inequality constraints given in Eqs. 15-22.

A. CASE-7: MINIMISING TOTAL GENERATION COST
The aim of this case study is to minimise the basic quadratic fuel cost given in Eq. 11. The fuel cost obtained by GWO algorithm is 20440.32 $/hr and this value is the best solution compared with those obtained using the GA, PSO, CSA, SHADE-SF and ABC algorithms, where the fuel cost value is 20919.9 $/hr by GA, 21475.1 $/hr by PSO, 20905.4 $/hr by CSA, 20786.5 $/hr by SHADE-SF and 20462.4 $/hr by ABC as given in Appendix Table 10. The convergence characteristics of GWO and the other optimisation techniques are shown in Fig. 22. The second case study in the IEEE-57 bus system aims to minimise both the quadratic fuel cost and carbon gas emission. It is similar to Case-2 in the IEEE-30 bus system when an additional emission constraint is included in the objective function given in Eq. 12. With higher penetration of RES, the value of emission in GWO is reduced from 11.87 ton/hr to 6.28 ton/hr. The convergence characteristics of GWO and the other techniques are shown in Fig. 23.
Simulation results in Appendix Table 11 show that the GWO and ABC algorithms are more efficient to find global optimum when compared to the other algorithms when using the IEEE-57 bus system. This is, perhaps, because of the GWO's and ABC's search mechanism that prevents them from easily getting trapped in local optima. Furthermore, the GWO algorithm requires the least computation time, suggesting that GWO is a highly promising algorithm for solving many practical global optimisation problems with computationally expensive objective function and constraints.

VII. CASE STUDIES AND RESULTS USING HYBRID ALGORITHMS
Generally, in a constrained optimisation problem, heuristic algorithms adopt premature convergence and tend to be computationally expensive. The OPF problem involves a large scale system and hence could yield impractically long execution times. For this reason, two hybrid techniques, namely GA-PSO and ABC-CSO, are presented for comparison. Note that each optimisation technique, GA, PSO, ABC and CSO, has its weaknesses and strengths. GA-PSO combines the properties of GA and PSO, while ABC-CSO combines ABC and CSO, in order to balance the exploration and exploitation capabilities. During the first iteration, both hybrid algorithms determine the individual best result of the relevant algorithm. In the second iteration, the best positions selected are mutated separately by applying different steps for both algorithms and the cost is therefore calculated. More details about these hybrid methods are found in [55] and [56]. In the present work, these hybrid algorithms are evaluated using the IEEE-30 and 57 bus systems. As expected, the hybrid models prove more successful with better search quality than the basic methods (i.e. GA, PSO, ABC and CSO). The advantage of hybrid approaches over basic techniques is their robustness and flexibility. The results obtained from the hybrid algorithms are good in terms of generation cost and are better in terms of execution time than the basic methods. Indeed, the worst solution in iteration one obtained by the hybrid methods is still better than the best result obtained by the basic methods in the last iteration (i.e. after 500 iterations). Finally, the comparison with the proposed GWO approach, for the IEEE-30 and 57 bus systems, are discussed in Cases 9 and 10, respectively. This case elaborates on the standard OPF problem with a basic quadratic cost function for the IEEE-30 bus system. The first objective is to minimise the total generation fuel cost given by Eq. 11. The hybrid algorithm results are summarised in Appendix Table 12 (cf. Tables 6 and 7 for GWO). The obtained cost for GWO is 781.40 $/hr, compared to 781.994 $/hr and 783.415 $/hr for GA-PSO and ABC-CSO,  respectively. When the execution times are compared, GWO remains effective, requiring only 429 seconds to complete 500 iterations, compared to 496 and 587 seconds for the GA-PSO and ABC-CSO algorithms. Fig. 24 compares the convergence characteristics of GWO and the hybrid algorithms. For this example, the hybrid algorithms converge faster, but all these techniques demonstrate fast and stable convergence characteristics.
The performance of the hybrid models is also tested when an additional emission constraint is added to the objective function i.e. Eq. 12. Simulation results in the Appendix Tables 7 and 12 show that GWO achieves similar minimum generation costs (809.93 $/hr) compared to the hybrid approaches. It is important to note that multiple parameters increase the complexity of the hybrid algorithms. Appendix Table 12 provides details of average execution time to complete one iteration. Using Appendix Tables 6 and 7, it is calculated that GWO requires 0.96 seconds to complete one iteration, compared to the higher execution times required by GA-PSO and ABC-CSO for similar scenarios. However, the convergence characteristics in this scenario are similar to those shown in Fig. 24, except that the generation cost is increased because of the carbon tax imposition.    In this case study, the performance of the hybrid algorithms is tested on the IEEE-57 bus system, again with and without tax imposition, with the results given by Appendix Table 12  (cf. Tables 10 and 11 for GWO). Without carbon tax, the generation costs are 20440.083 $/hr, 20440.112 $/hr and 20440.32 $/hr for GA-PSO, ABC-CSO and GWO, respectively. The computational time taken by GA-PSO, ABC-CSO and GWO are 520, 1006 and 440 seconds, respectively. Fig. 25 illustrates the convergence on the IEEE-57 bus system for each algorithm. Table 12 also shows results when the emission constraint is added i.e. Eq. 12. By imposing a carbon tax at the rate of $20/ton, the carbon emissions have been significantly reduced by GWO and the two hybrid algorithms.
The optimum results given in Appendix Tables 6, 7, 10, 11 and 12 show that the hybrid algorithms can potentially achieve a better result than GWO but (for these examples)   at the cost of computational time. According to the ''no free lunch theorem of optimisation'' [57], there is no single optimisation technique which is best suited to solve all kinds of optimisation problem. In the present context, hybrid algorithms can be considered as a feasible solution for different OPF problems where generation cost saving is a priority. Whilst the present article has focused on the GWO approach, it can be pointed out that hybrid algorithms might also be made more computationally efficient, motivating further research into the utility of such approaches for OPF.

VIII. CONCLUSION
In this work, a recently developed evolutionary algorithm, GWO, was employed to optimise OPF problems whilst considering stochastic RES in the network. Different PDFs were used to model SPG and WPG uncertainty, and their integration methods were discussed. A number of case studies were investigated to evaluate the performance of the proposed algorithm and the results were compared with other well recognised evolutionary algorithms. Hence, novel contributions include the proposed objective functions that consider RES; the use of a GWO approach to address the non-convex OPF problem, and its application both in small and medium-scale systems with evaluation via simulation.
The safety of an electrical network is compromised if physical or security constraints on system components are compromised. Such a situation may lead to excessive losses, malfunctioning of the components and sometimes complete failure of the system. It is essential that the network runs within predefined limits. The new results in this article show VOLUME 8, 2020  the GWO proves to be very effective and reliable, with fast convergence rates to find global solutions for the considered objective functions. It outperforms other algorithms in terms of total cost and convergence time minimisation, whilst simultaneously addressing the necessary system constraints.
In this regard, the other algorithms sometimes adopt premature convergence, which can stop the algorithm from finding a global solution. By contrast, for the scenarios considered in this article, GWO maintains a satisfactory balance between exploration and exploitation, in order to find a global solution. Furthermore, when the elapsed time of GWO and the benchmark algorithms are compared, GWO remains very effective. Hence, the results in this article suggest that GWO could be applied to various non-linear, non-convex, multimodal and constrained optimisation problems in OPF.
More generally, the present work adds to the growing literature that suggests GWO offers a flexible and powerful tool for certain types of optimisation problem. In a different context, for example, GWO has been used for research into robot control systems. Indeed, the present authors are investigating future use of GWO in this area for nuclear decommissioning applications (see Bandala et al. [58]), with  Tables 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11 and Table 12. XIANDONG MA received the B.Eng. degree in electrical engineering, in 1986, the M.Sc. degree in power systems and automation, in 1989, and the Ph.D. degree in partial discharge-based highvoltage condition monitoring, in 2003. He is a Senior Lecturer with the Department of Engineering, Lancaster University, U.K. His research interests include intelligent condition monitoring and fault diagnosis of the renewable energy systems, distributed energy generation and smart grids, advanced signal processing, and instrumentation. He has published over 110 papers in the leading international journals, book chapters, and international conferences, and holds patents. He is a Chartered Engineer (C.Eng.) and a Fellow of the Institution of Engineering and Technology (FIET). He is also an Associate Editor of the IEEE SYSTEMS JOURNAL and the International Journal of Automation and Computing. VOLUME 8, 2020