*Article* **Hybrid Gravitational–Firefly Algorithm-Based Load Frequency Control for Hydrothermal Two-Area System**

**Deepak Kumar Gupta <sup>1</sup> , Ankit Kumar Soni <sup>1</sup> , Amitkumar V. Jha <sup>2</sup> , Sunil Kumar Mishra <sup>2</sup> , Bhargav Appasani <sup>2</sup> , Avireni Srinivasulu <sup>3</sup> , Nicu Bizon 4,5,6,\* and Phatiphat Thounthong 7,8**


**Abstract:** The load frequency control (LFC) and tie-line power are the key deciding factors to evaluate the performance of a multiarea power system. In this paper, the performance analysis of a two-area power system is presented. This analysis is based on two performance metrics: LFC and tie-line power. The power system consists of a thermal plant generation system and a hydro plant generation system. The performance is evaluated by designing a proportional plus integral (PI) controller. The hybrid gravitational search with firefly algorithm (hGFA) has been devised to achieve proper tuning of the controller parameter. The designed algorithm involves integral time absolute error (ITAE) as an objective function. For two-area hydrothermal power systems, the load frequency and tie-line power are correlated with the system generation capacity and the load. Any deviation in the generation and in the load capacity causes variations in the load frequencies, as well as in the tie-line power. Variations from the nominal value may hamper the operation of the power system with adverse consequences. Hence, performance of the hydrothermal power system is analyzed using the simulations based on the step load change. To elucidate the efficacy of the hGFA, the performance is compared with some of the well-known optimization techniques, namely, particle swarm optimization (PSO), genetic algorithm (GA), gravitational search algorithm (GSA) and the firefly algorithm (FA).

**Keywords:** load frequency control; automatic generation control; controllers; optimization techniques; multisource power system; interconnected power system; hybrid gravitational with fire fly algorithm; gravitational search algorithm; firefly algorithm

#### **1. Introduction**

A power network generally comprises several areas or power systems, interconnected through tie-lines. Distribution systems, transmission lines, and generation systems that may also include renewable energy sources are some of the prime constituents of the power network [1]. The real-time integration of these components and their operation in the dynamic environment cause differences in the active and reactive power demands. The variations in these quantities produce undesired oscillations in the system. These

**Citation:** Gupta, D.K.; Soni, A.K.; Jha, A.V.; Mishra, S.K.; Appasani, B.; Srinivasulu, A.; Bizon, N.; Thounthong, P. Hybrid Gravitational–Firefly Algorithm-Based Load Frequency Control for Hydrothermal Two-Area System. *Mathematics* **2021**, *9*, 712. https://doi.org/10.3390/math9070712

Academic Editor: Aleksandr Rakhmangulov

Received: 8 February 2021 Accepted: 22 March 2021 Published: 25 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

oscillations have to be damped, else they may adversely affect the operation of the power system, and may even lead to a power blackout. Many approaches have been adopted in the literature in the domain of power grid control and stability [2,3]. Y. Li et al. addressed the issues of scalability, privacy, and reliability in a multienergy system [4]. For a power system comprising multiple areas, the issue of optimal generation and distribution has been addressed in [5]. The power management of interconnected single-phase/three-phase microgrids for enhancing voltage quality is considered by J. Zhou et al. in [6]. Researchers have also modernized the power system using smart grid technologies with the objective to make the power system reliable, resilient, secure, and stable [7–10].

For an interconnected power system with multiple sources, the automatic generation control (AGC) methodology has been devised to limit the oscillations produced in the power system due to the mismatch in demand and supply. However, to ensure the control and stable operation of the power system, it is desirable that the oscillations lie within the acceptable range. Further, these oscillations must be controlled within a minimum time to stabilize the system. The main motivation of the AGC is to improve the performance of the interconnected power system by considering several performance metrics, such as load frequency and tie-line power. These performance metrics are highly correlated with generation capacity at the generation side, demand at the consumer side, and total losses at the transmission side. Any mismatch in these may result in deviations in load frequency, as well as in the tie-line power flow. This may lead the system to an unstable state, with severe consequences. Thus, system load frequencies and tie-line power (TLP) must be within the nominal range to realize a stable system. This is generally achieved using the load frequency control (LFC) method. Further, the power systems are stabilized by controlling the speed of the generators (for load frequency) and the TLP based on the area control system. The area control system has basically two objectives: to cater to the demands of its own customers and to respond to the demands of other control areas. In the context of area control systems, area control error (ACE) comprises load frequencies (LF) and tie-line power (TLP).

Automatic generation control has the following major responsibilities [11–13]:


Extensive research in the literature shows an attempt to bridge the gap in the modeling and analysis of the hybrid hydrothermal power system (HTPS) through a linearized approach. The linearized model of the power system is easy for performance evaluation. This has further motivated researchers to present the linearized model of the interconnected power system for AGC analysis. The power system with more than one interconnected area is referred to as a multiarea power system (MAPS). Many aspects of the AGC analysis in the case of interconnected power systems have been thoroughly discussed in [14,15] with several case studies. For the design of the control parameters of proportional-integral-derivative (PID) and proportional plus integral (PI) controllers, the maximum peak resonance method has been reported in [16]. Continuous and discrete mode analysis with a generation rate constraint for interconnected HTPS has been reported in [17]. Problems with load frequency control for hybrid hydrothermal have also been addressed in [18]. Further, Jha et al. have considered the PI controller for the load frequency control of hybrid hydrothermal systems [19]. The use of several artificial optimization techniques were illustrated with case studies in [20] for the LFC in various system operating conditions. A comparative analysis of several soft computing techniques for the LFC have been studied in [21] by Gupta et al. Further, for a power system with interconnected areas, load frequency control analysis has been carried out in [22] using a hybrid adaptive gravitational search and pattern search algorithm. Recently, Gupta et al. reported novel hybrid optimization techniques for addressing the issues of LFC in MAPS comprising multiple sources [23]. Koley et al. [24] presented the issue of LFC by considering a power system involving hybrid power plants such as thermal, wind, and photovoltaic

generation stations. The exemplary work of Khadange et al. introduces the hybrid guided gravitational search with pattern search (hGGSA-PS) optimization technique for MAPS in order to analyze the LFC [25]. AGC, using a coordinated design for two-area systems (TAS), was proposed in [26] by Khezri et al. The scope of some of the advanced controllers to achieve AGC for multiarea systems was considered by Gondaliya et al. in [27]. A new approach, referred to as secondary LFC, has been introduced in [28] for a power system with a multigrid configuration. A new fractional order PI controller was proposed by Celik et al. in [29]. Various optimization techniques for controlling the controller parameter for the LFC of multiarea power systems have been proposed in the literature [30]. For example, a wind-driven optimization algorithm has been proposed by Haes et al. in [31]. A social spider optimization technique was presented in [32]. Further, Nilkmanesh et al. have proposed a multiobjective uniform-diversity genetic algorithm (MUGA) for MAPS in [33]. The slap-swarm optimization technique was discussed for the LFC of MAPS by Sahu et al. in [34].

The overall objective of the design is to discuss the LFC of the interconnected MAPS. The novelty and contribution of this study can be highlighted as follows:


Thus, in summary, this paper presents a hybrid algorithm to tune the PI controller for the optimum LFC of an interconnected MAPS. This hybrid technique merges two well-known optimization techniques, the GSA and the FA. Effectiveness of the proposed controller tuned using this hybrid intelligent optimization technique is compared with controllers tuned using other well-known optimization techniques, such as the GA, the PSO, the GSA, and the FA. The proposed algorithm works well for different operating conditions (such as changes in load), which shows its robustness. The dynamic response of all the state variables has been improved in terms of settling time and overshoot. It is observed that the proposed controller outperformed the other techniques (GA, PSO, GSA, and FA) in terms of performance, stability, and robustness.

The remainder of this paper is presented in five sections. In Section 2, a brief introduction about the proportional integral controller is presented. The test system of hybrid HTPS is also presented in this section. Further, to analyze the test system, it is modeled on the basis of the state space approach. The preliminaries of the optimization techniques along with the design of hGFA for a hydrothermal power system under consideration is presented in Section 3. Section 4 covers the design methodology and the simulation of the HTPS. The analysis of the simulation results and discussion is presented in Section 5. Finally, Section 6 deals with the conclusion of the research work. Some of the system variables are summarized in the Appendix A.

#### **2. State Space Modeling of a Hydrothermal System**

This section deals with the PI controller and the modeling of the two-area hydrothermal system using a state space approach.

#### *2.1. The Proportional Plus Integral Controller*

We have considered the PI controller for a LFC analysis of the hydrothermal power system. If *e*(*t*) represents the error at the input of the PI controller, then its output *u*(*t*) can be represented by Equation (1) as given below:

$$u(t) = K\_{pr}e(t) + K\_i \int e(t)dt\tag{1}$$

where *Kpr*, and *K<sup>i</sup>* , are the controller parameters. In this paper, we consider the objective function as an integral time absolute error (*ITAE*) [23]. The *ITAE* cost function is as given by Equation (2) below.

$$ITAE = \int\_0^\infty t \times |e(t)| dt\tag{2}$$

#### *2.2. Test Model for the Two-Area Hydrothermal System*

In order to obtain a performance analysis of the hydrothermal power system based on hGFA optimization techniques, the test system as illustrated in Figure 1 is considered. The different blocks of the test system are modeled using the standard linearized method and are shown by their respective transfer functions. Here, the two-area power system is considered such that one area consists of a hydropower plant and the other area has a thermal power plant for power generation. Each plant has a PI controller, shown by the transfer function to achieve automatic gain control of the power system.

**Figure 1.** Test system model using transfer functions.

The model of the test system is obtained using the state space approach detailed in [11]. The state space modeling of the test system and its implementation are discussed in the next subsection.

#### *2.3. State Space Modeling*

As shown in Figure 1, the test system comprises two interconnected areas resulting in a hydrothermal power system. Each area consists of one PI controller to achieve AGC based on LF, as well as TLP for effective operation of the power system. With the given model and its transfer functions, the test system can be modeled through multidimensional state space analysis. Using a state space approach corresponding to Equations (3) and (4), the transfer function of the overall hydrothermal power system can be obtained.

$$
\dot{\mathfrak{x}} = A\mathfrak{x} + Bu \tag{3}
$$

$$y = \mathbb{C}\mathfrak{x} + Du \tag{4}$$

Here, *x*, *y*, and *u*, denotes the state vector, output vector, and input vector, respectively. The different matrices with the real constant model-dependent values of the state space model are represented by variables *A*, *B*, *C*, and *D*.

On analysis of the state space model, the following equations follow:

.

$$\propto = \left[\Delta f\_{1\prime} \Delta P\_{mech1\prime} \Delta P\_{v1\prime} \Delta f\_{2\prime} \Delta P\_{mech2\prime} \Delta P\_{1\prime} \Delta P\_{v2\prime} \Delta P\_{12\prime} \Delta P\_{ref1\prime} \Delta P\_{ref2}\right]^T \tag{5}$$

$$\boldsymbol{\mu} = \begin{bmatrix} \Delta \boldsymbol{P}\_{L1} \Delta \boldsymbol{P}\_{L2} \end{bmatrix}^{\mathrm{T}} \tag{6}$$

$$y = \begin{bmatrix} \Delta f\_{1\prime} \Delta f\_{2\prime} \Delta P\_{12} \end{bmatrix}^T \tag{7}$$

These variables are part of the test system's state space model, which characterizes the entire hydrothermal power system as described in [11]. The state space matrices to construct the transfer function is determined by analyzing the differential equations. These are illustrated below.

$$\dot{\mathbf{x}}\_1 = -\frac{1}{T\_{p1}}\mathbf{x}\_1 + \frac{\mathbf{K}\_{p1}}{T\_{p1}}\mathbf{x}\_2 - \frac{\mathbf{K}\_{p1}}{T\_{p1}}\mathbf{x}\_8 - \frac{\mathbf{K}\_{p1}}{T\_{p1}}\mathbf{u}\_1 \tag{8}$$

$$
\dot{\mathbf{x}}\_2 = -\frac{1}{T\_{t1}}\mathbf{x}\_2 + \frac{1}{T\_{t1}}\mathbf{x}\_3 \tag{9}
$$

$$\dot{\mathbf{x}}\_3 = -\frac{1}{R\_1 T\_{\mathcal{S}1}} \mathbf{x}\_1 - \frac{1}{T\_{\mathcal{S}1}} \mathbf{x}\_3 + \frac{1}{T\_{\mathcal{S}1}} \mathbf{x}\_9 \tag{10}$$

$$\dot{\mathbf{x}}\_4 = -\frac{1}{T\_{p2}}\mathbf{x}\_4 + \frac{\mathbf{K}\_{p2}}{T\_{p2}}\mathbf{x}\_5 + \frac{\mathbf{K}\_{p2}}{T\_{p2}}\mathbf{x}\_8 - \frac{\mathbf{K}\_{p2}}{T\_{p2}}\mathbf{u}\_2 \tag{11}$$

$$\dot{\mathbf{x}}\_{5} = -\frac{2T\_{2}}{R\_{2}T\_{1}T\_{3}}\mathbf{x}\_{4} - \frac{2}{T\_{w}}\mathbf{x}\_{5} + \left(\frac{2}{T\_{w}} + \frac{2}{T\_{3}}\right)\mathbf{x}\_{6} + \left(\frac{2T\_{2}}{T\_{1}T\_{3}} - \frac{2}{T\_{3}}\right)\mathbf{x}\_{7} - \frac{2T\_{2}}{T\_{1}T\_{3}}\mathbf{x}\_{10} \tag{12}$$

$$\dot{\mathbf{x}}\_6 = -\frac{T\_2}{R\_2 T\_1 T\_3} \mathbf{x}\_4 - \frac{1}{T\_3} \mathbf{x}\_6 + \left(\frac{1}{T\_3} - \frac{T\_2}{T\_1 T\_3}\right) \mathbf{x}\_7 + \frac{T\_2}{T\_1 T\_3} \mathbf{x}\_{10} \tag{13}$$

$$
\dot{\mathbf{x}}\_7 = -\frac{1}{R\_2 T\_1} \mathbf{x}\_4 - \frac{1}{T\_1} \mathbf{x}\_7 + \frac{1}{T\_1} \mathbf{x}\_{10} \tag{14}
$$

$$
\dot{\mathbf{x}}\_8 = T\_s \mathbf{x}\_1 - T\_s \mathbf{x}\_4 \tag{15}
$$

$$\dot{\mathbf{x}}\_{9} = \left(\frac{\mathbf{B}\_{1}\mathbf{K}\_{\text{pr1}}}{T\_{p1}} - \mathbf{K}\_{\text{pr1}}T\_{\text{s}} - \mathbf{K}\_{\text{il}}\mathbf{B}\_{1}\right)\mathbf{x}\_{1} - \frac{\mathbf{B}\_{1}\mathbf{K}\_{\text{pr1}}\mathbf{K}\_{\text{p1}}}{T\_{p1}}\mathbf{x}\_{2} + \mathbf{K}\_{\text{pr1}}T\_{\text{s}}\mathbf{x}\_{4} + \left(\frac{\mathbf{B}\_{1}\mathbf{K}\_{\text{pr1}}\mathbf{K}\_{\text{p1}}}{T\_{p1}} - \mathbf{K}\_{\text{il}}\right)\mathbf{x}\_{8} + \frac{\mathbf{B}\_{1}\mathbf{K}\_{\text{pr1}}\mathbf{K}\_{\text{p1}}}{T\_{p1}}\mathbf{u}\_{1} \tag{16}$$

.

$$\mathbf{x}\_{10}^{\cdot} = \mathbf{K}\_{p2}\mathbf{T}\_{\mathbf{x}}\mathbf{x}\_{1} + \left(\frac{\mathbf{B}\_{2}\mathbf{K}\_{p2}}{T\_{p2}} - \mathbf{K}\_{p2}\mathbf{T}\_{\mathbf{x}} - \mathbf{K}\_{\mathbf{Z}}\mathbf{B}\_{2}\right)\mathbf{x}\_{4} - \frac{\mathbf{B}\_{2}\mathbf{K}\_{p2}\mathbf{K}\_{p2}}{T\_{p2}}\mathbf{x}\_{5} + \left(-\frac{\mathbf{B}\_{2}\mathbf{K}\_{p2}\mathbf{K}\_{p2}}{T\_{p2}} + \mathbf{K}\_{\mathbf{Z}}\right)\mathbf{x}\_{8} + \frac{\mathbf{B}\_{2}\mathbf{K}\_{p2}\mathbf{K}\_{p2}}{T\_{p2}}\mathbf{u}\_{2} \tag{17}$$

Finally, the state matrices are obtained by incorporating the above differential equations. These are as shown below:

*A* = − <sup>1</sup> *TP*<sup>1</sup> *KP*<sup>1</sup> *TP*<sup>1</sup> 0 0 0 0 0 − *KP*<sup>1</sup> *TP*<sup>1</sup> 0 0 0 − <sup>1</sup> *Tt*<sup>1</sup> 1 *Tt*<sup>1</sup> 0 0 0 0 0 0 0 − <sup>1</sup> *R*1*Tg*<sup>1</sup> 0 1 *Tg*<sup>1</sup> 0 0 0 0 0 <sup>1</sup> *Tg*<sup>1</sup> 0 0 0 0 − <sup>1</sup> *TP*<sup>2</sup> *KP*<sup>2</sup> *TP*<sup>2</sup> 0 0 *<sup>K</sup>P*<sup>2</sup> *TP*<sup>2</sup> 0 0 0 0 0 <sup>2</sup>*T*<sup>2</sup> *R*2*T*1*T*<sup>2</sup> − <sup>2</sup> *Tw* ( 2 *Tw* + <sup>2</sup> *T*3 ) ( <sup>2</sup>*T*<sup>2</sup> *T*1*T*<sup>3</sup> − <sup>2</sup> *T*3 ) 0 0 − 2*T*<sup>2</sup> *T*1*T*<sup>3</sup> 0 0 0 − *T*2 *R*2*T*1*T*<sup>2</sup> 0 − <sup>1</sup> *T*3 ( 1 *T*3 − *T*2 *T*1*T*<sup>3</sup> ) 0 0 *<sup>T</sup>*<sup>2</sup> *T*1*T*<sup>3</sup> 0 0 0 − <sup>1</sup> *R*2*T*<sup>1</sup> 0 0 − <sup>1</sup> *T*1 0 0 <sup>1</sup> *T*1 *T<sup>s</sup>* 0 0 −*T* 0 0 0 0 0 0 *A*9,1 − *B*1*Kpr*1*KP*<sup>1</sup> *TP*<sup>1</sup> 0 *Kpr*1*T<sup>s</sup>* 0 0 0 ( *B*1*Kpr*1*KP*<sup>1</sup> *TP*<sup>1</sup> − *Ki*1) 0 0 *Kpr*2*T<sup>s</sup>* 0 0 *A*10,4 − *B*2*Kpr*2*KP*<sup>2</sup> *TP*<sup>2</sup> 0 0 (− *B*2*Kpr*2*KP*<sup>2</sup> *TP*<sup>2</sup> + *Ki*2) 0 0 

where,

$$\begin{aligned} A\_{9,1} &= (\frac{B\_1 K\_{pr1}}{T\_{P1}} - K\_{pr1} T\_s - K\_{l1} B\_1), \ A\_{10,4} = (\frac{B\_2 K\_{pr2}}{T\_{P2}} - K\_{pr2} T\_s - K\_{l2} B\_2) \\ B &= \begin{bmatrix} -\frac{K\_{pr1}}{T\_{P1}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{B\_1 K\_{pr1} K\_{P1}}{T\_{P1}} & 0 \\ 0 & 0 & 0 & -\frac{K\_{P2}}{T\_{P2}} & 0 & 0 & 0 & 0 & 0 & \frac{B\_2 K\_{p2} K\_{P2}}{T\_{P2}} \end{bmatrix}^T \\ C &= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix} D = 0 \end{aligned}$$

The transfer function (*T*.*F*) using matrices *A*, *B*, *C*, and *D* can be obtained using the relation:

$$T.F = \mathcal{C}[SI - A]^{-1}B + D$$

Without loss of generality, the transfer function will be a matrix of the dimension (3 × 2). The elements of the transfer function are as follows.

$$
\Delta f\_1 = \frac{\Delta f\_1(s)}{\Delta P\_{L1}(s)} \Delta P\_{L1} + \frac{\Delta f\_1(s)}{\Delta P\_{L2}(s)} \Delta P\_{L2} \tag{18}
$$

$$
\Delta f\_2 = \frac{\Delta f\_2(s)}{\Delta P\_{L1}(s)} \Delta P\_{L1} + \frac{\Delta f\_2(s)}{\Delta P\_{L2}(s)} \Delta P\_{L2} \tag{19}
$$

$$
\Delta P\_{12} = \frac{\Delta P\_{12}(s)}{\Delta P\_{L1}(s)} \Delta P\_{L1} + \frac{\Delta P\_{12}(s)}{\Delta P\_{L2}(s)} \Delta P\_{L2} \tag{20}
$$

#### **3. Hybrid Gravitational–Firefly Algorithm**

The hybrid gravitational–firefly algorithm uses two well-known optimization techniques, the GSA and the FA [35–38]. Newton's gravitational law is the working principle for the GSA. In this, every object is treated as a candidate solution. The masses of each of these variables or the candidate are used to evaluate their performance depending on the value of the selected objective function [21]. On the other hand, the working of the FA is similar to the flashing behavior of fireflies, which they use to achieve communication amongst themselves. The hybrid gravitational–firefly algorithm takes the properties of both these algorithms and updates the values of its objects (candidate solutions) using the updated equations.

#### *3.1. GSA*

The position (*x*) of each agent (*i*) out of *N* agents is represented by *x<sup>i</sup>* . If we assume *m* dimensional space, then the position of each agent *x d i* is given using the following equation:

$$\mathbf{x}\_{i} = (\mathbf{x}\_{i}^{1}, \dots, \mathbf{x}\_{i}^{d}, \dots \mathbf{x}\_{i}^{m}) \\ \forall i = \{1, 2, \dots, N\} \tag{21}$$

The mass '*j*' applies force on mass '*i*', which is given as,

$$F\_{ij}^d(t) = G(t) \times \frac{\left(\mathbf{x}\_j^d(t) - \mathbf{x}\_i^d(t)\right)}{R\_{ij}(t) + \varepsilon} \left(M\_{pi}(t) \* M\_{aj}(t)\right) \tag{22}$$

where *G*(*t*) denotes the gravitational constant, *Rij*(*t*) denotes the Euclidian distance between agent '*i*' and agent '*j*', and *Mpi*(*t*) and *Mai*(*t*) represent the passive gravitational and active gravitational masses of the agent '*i*' and the agent '*j*', respectively. The Euclidian and the gravitational constants are given as below.

$$\begin{array}{c} \mathcal{R}\_{ij}(t) = \left| \left| \left. x\_i(t), x\_j(t) \right| \right|\_2 \right. \\ G(t) = G(G\_0, t) \end{array} \tag{23}$$

If *Mii*(*t*) indicates the inertial mass of the *i*th agent, the total force and the total acceleration acting on an agent '*i*' due to other agents in the *d*-dimensional space is given as:

$$F\_i^d(t) = \sum\_{j=1, j \neq i}^{N} rand\_j F\_{ij}^d(t),\ a\_i^d(t) = \frac{F\_i^d(t)}{M\_{ii}(t)}\tag{24}$$

In each round, Equations (25) and (26) are used to update the velocity as well as the mass of each agent *i*.

$$m\_l(t) = \frac{\text{fit}\_l(t) - \text{worst}(t)}{\text{best}(t) - \text{worst}(t)} \tag{25}$$

$$M\_i(t) = \frac{m\_i(t)}{\sum\_{j=1}^{N} m\_j(t)}\tag{26}$$

#### *3.2. Firefly Algorithm*

The Euclidian distance between the firefly '*i*' and firefly '*j*' for given position *x<sup>i</sup>* and *x<sup>j</sup>* can be evaluated using the equation given below,

$$r\_{ij} = \sqrt{\sum\_{k=1}^{d} \left(\varkappa\_{i,k} - \varkappa\_{j,k}\right)^2} \tag{27}$$

where *k* indicates the *k*th element of the spatial coordinates.

The attractiveness parameter of every firefly basically can be given by the following equation:

$$
\beta = \beta\_0 e^{-\gamma r^2} \tag{28}
$$

where *γ* represents the coefficient of absorption, which is used to control the light intensity. We can define the movement of each firefly as follows.

$$v\_i^d(t+1) = \text{rand}\_i \times v\_i^d(t) + a\_i^d(t) \tag{29}$$

$$\mathbf{x}\_{l}^{d} = \mathbf{x}\_{l}^{d} + \beta\_{0}e^{-\gamma r\_{ij}^{2}}(\mathbf{x}\_{l} - \mathbf{x}\_{l}) + v\_{l}^{d}(t+1) + a\varepsilon \tag{30}$$

Here, *x<sup>i</sup>* is the instantaneous position of an object and *αε* indicates the random behavior of a firefly if no brighter firefly is detected.

#### **4. Methodology and Simulation Results**

This section presents a methodology to design and configure two-area hydrothermal power systems using simulation to evaluate the performance of the proposed optimization technique.

#### *4.1. Simulation Methodology*

To evalute the efficacy of the proposed hybrid gravitational–firefly algorithm, the test system shown in Figure 1 is designed. The test system comprises a two-area power system, one with a thermal power plant and the other with a hydropower plant. To design and configure the two-area power system, MATLAB Simulink is used. Further, the proposed algorithm is coded in a MATLAB script file and then interfaced to the Simulink model for testing its efficacy on the test system.

The various simulation configuration parameters are as follows: MATLAB (R2016a) software is used along with the Simulink tool. The system used for simulation has an i5-6200 CPU@ processor running at 2.30 GHz frequency and having 8 GB RAM. The proposed algorithm and other optimization techniques are written as MATLAB scripts, which are interfaced to the test power system through the Simulink. A few of the other key parameters are included in the Appendix A section.

#### *4.2. Simulation Results*

The simulation was performed to evaluate the load frequencies and tie-line parameters of the test system. The simulations were carried out by tuning the parameters of the PI controller using the novel hybrid gravitational search with firefly algorithm as an optimization technique. To understand the efficacy of this algorithm, the results are compared with PSO, GA, GSA, and FA, which can be used as a benchmark for performance evaluation. In all the optimization techniques, the cost function is ITAE and remains unaltered.

The performance of the proposed algorithm on hydrothermal power systems is evaluated by considering two case studies, which are discussed below.

#### 4.2.1. Case Study-I

To observe performance under a step load change, the load in the area having the thermal power plant (area-1) is incremented up to 20%. The load in the area having the hydropower plant (area-2) is unchanged. The test system is simulated using different optimization techniques, i.e., PSO, GA, GSA, FA, and hGFA. ITAE is used as the objective function. The parameters of the PI controller tuned using these optimization techniques are summarized in Table 1.


**Table 1.** Optimized parameters of PI controller for case study-I.

The performance of the test system based on LF as well as TLP is analyzed by using the optimized parameters shown in Table 1. The simulation results corresponding to the load frequency in area-1 (∆*f*1), the load frequency in area-2 (∆*f*2), and the tie-line power flow of the test system are shown in Figures 2–4, respectively.

∆<sup>1</sup> ∆<sup>2</sup> –

∆<sup>1</sup> ∆<sup>2</sup> –

**Figure 2.** Perturbation in load frequency response in area-1 of HTPS for case study-I.

**Figure 3.** Perturbation in load frequency response in area-2 of HTPS for case study **Figure 3.** Perturbation in load frequency response in area-2 of HTPS for case study-I.

**Figure 4.** Perturbation in the TLP flow response of HTPS for case study-I.

#### 4.2.2. Case Study-II

To observe performance under a step load change, the load in the area with a thermal power plant (area-1) is incremented up to 5% and the load in area-2 is subjected to a 1% change. The tuned parameters of the PI controller using different optimization techniques are summarized in Table 2.



∆<sup>1</sup> ∆<sup>2</sup> – The simulation results corresponding to the load frequency in area-1 (∆*f*1), the load frequency in area-2 (∆*f*2), and the tie-line power flow of the test system are shown in Figures 5–7, respectively. ∆<sup>1</sup> ∆<sup>2</sup> –

**Figure 5.** Perturbation in load frequency response in area-1 of HTPS for case study-II.

**Figure 6.** Perturbation in load frequency response in area-2 of HTPS for ca **Figure 6.** Perturbation in load frequency response in area-2 of HTPS for case study-II.

**Figure 7.** Perturbation in the TLP flow response of HTPS for case study-II.

#### **5. Analysis and Discussion**

The simulated results of the test system are analyzed in this section. To perform a comparative analysis, the results achieved using the hGFA are compared with those obtained using the PSO, the GA, the GSA, and the FA. The comparisons are in terms of the overshoot/undershoot, and settling time for the load frequency response and the TLP.

#### *5.1. Case Study-I*

For the first case study, for the perturbation in step load change, the load frequency responses in area-1 and area-2 have been shown in Figures 2 and 3, respectively. Figure 4 represents the TLP for the interconnected power system belonging to case study-I. It can be inferred from these responses that the hGFA optimization technique outperforms other optimization techniques. The overshoot/undershoot, and settling time of these results are summarized in Table 3.


**Table 3.** Comparison of optimization techniques based on overshoot and settling time for case study-I.

From Table 3, it can be seen that the overshoot/undershoot for the system variable ∆f<sup>1</sup> is −0.46, −0.47, −0.5, −0.46, and −0.47 for the hGFA, the PSO, the GA, the GSA, and the FA, respectively. Further, it can also be seen that the overshoot/undershoot for the system variable ∆f<sup>2</sup> is −0.50, −0.55, −0.6, −0.52, and −0.51, corresponding to the hGFA, the PSO, the GA, the GSA, and the FA, respectively. Furthermore, the system variable ∆P<sup>12</sup> is −0.11, −0.12, −0.12, −0.12, and −0.11, corresponding to the hGFA, the PSO, the GA, the GSA, and the FA, respectively. Thus, it can be inferred that the hGFA outperforms other well-known optimization techniques such as PSO, GA, GSA, and FA in terms of overshoot/undershoot for all the system variables of the HTPS. Similarly, it can be inferred by analyzing the results in Table 3 that the hGFA (with *T<sup>s</sup>* = 35 s, 33 s, 35 s, respectively, for ∆f1, ∆f2, ∆P12) performs better in terms of settling time (*Ts*) compared to its competitor optimization techniques such as PSO (with *T<sup>s</sup>* = 40 s, 38 s, 39 s, respectively, for ∆f1, ∆f2, ∆P12), GA (with *T<sup>s</sup>* = 40 s, 41 s, 45 s, respectively, for ∆f1, ∆f2, ∆P12), GSA (with *T<sup>s</sup>* = 40 s, 40 s, 43 s, respectively, for ∆f1, ∆f2, ∆P12), and FA (with *T<sup>s</sup>* = 45 s, 43 s, 42 s, respectively, for ∆f1, ∆f2, ∆P12). Hence, it can be concluded that, for case study-I, the performance of the proposed hybrid optimization technique is better than other optimization techniques.

#### *5.2. Case Study-II*

From the perturbation in responses of frequencies in both the area as well as tie-line power, illustrated in Figures 5–7, it can be inferred that the hGFA optimization technique performs better compared to other techniques. The overshoot/undershoot and settling time for different system variables are summarized in Table 4.

From Table 4, it can be seen that the overshoot/undershoot for the system variable ∆f<sup>1</sup> is −0.12, −0.158, −0.13, −0.15, and −0.14, corresponding to hGFA, PSO, GA, GSA, and FA, respectively. Further, it can also be seen that the overshoot/undershoot for the system variable ∆f<sup>2</sup> is −0.17, −0.2, −0.17, −0.15, and −0.12, corresponding to hGFA, PSO, GA, GSA, and FA, respectively. The system variable ∆P<sup>12</sup> is −0.02, −0.021, −0.02, −0.018, and −0.014, corresponding to hGFA, PSO, GA, GSA, and FA, respectively. Similarly, by analyzing Table 4, it can be seen that the hGFA (with *T<sup>s</sup>* = 35 s, 30 s, 30 s, respectively, for ∆f1, ∆f2, ∆P12) performs better in terms of settling time (*Ts*) compared to its competitor optimization techniques such as PSO (with *T<sup>s</sup>* = 70 s, 60 s, 45 s, respectively, for ∆f1, ∆f2, ∆P12), GA (with *T<sup>s</sup>* = 75 s, 70 s, 70 s, respectively, for ∆f1, ∆f2, ∆P12), GSA (with *T<sup>s</sup>* = 65 s, 65 s, 70 s, respectively, for ∆f1, ∆f2, ∆P12), and FA (with *T<sup>s</sup>* = 70 s, 60 s, 70 s, respectively, for ∆f1, ∆f2, ∆P12). These algorithms offer similar overshoot/undershoot characteristics, but the settling time characteristics vary. The proposed algorithm reduces the settling time considerably, by approximately 50% when compared to other techniques.


**Table 4.** Comparison of optimization techniques based on overshoot/undershoot and settling time for case study-II.

#### **6. Conclusions**

The TAS consideration of hydrothermal systems is considered in this paper for LFC as well as TLP analysis. Each area has been considered with a proportional integral controller for analysis. The parameters of each PI controller are tuned using several optimization techniques with ITAE as an objective function. From the in-depth simulation presented in Section 4, followed by the analysis in Section 5, it has been observed that the proposed hGFA performs better than the other optimization techniques. Further, the hGFA improves the response of the interconnected power system, resulting in low overshoot and reduced settling time for LF as well as TLP flow. The main findings of this study can be summarized as follows:


Hence, it can be inferred that the hGFA outperforms the other optimization techniques (PSO, GA, GSA, and FA). Therefore, it can be successfully applied for AGC of a hydrothermal power system.

**Author Contributions:** Conceptualization, D.K.G., A.K.S. and B.A.; methodology, D.K.G. and A.K.S.; software, D.K.G.; validation, A.V.J.; investigation, D.K.G.; resources, B.A.; data curation, A.V.J.; writing—original draft preparation, A.V.J. and A.S.; supervision, A.S. and N.B.; project administration, A.S.; formal analysis: N.B.; funding acquisition: N.B. and P.T.; visualization: P.T.; writing—review

and editing: N.B., P.T.; figure and table, S.K.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the International Research Partnerships: Electrical Engineering Thai–French Research Center (EE-TFRC) between King Mongkut's University of Technology North Bangkok and the University of Lorraine under grant KMUTNB−BasicR−64−17.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Some of the standards values of the system variables are as follows:

*Tp*<sup>1</sup> = *Tp*<sup>2</sup> = 20 s, *T<sup>g</sup>* = 0.08 s, *T<sup>t</sup>* = 0.3 s, *B*<sup>1</sup> = *B*<sup>2</sup> = 0.4249, *T*<sup>12</sup> = 0.0866, *T*<sup>1</sup> = 48.709, *T*<sup>2</sup> = 0.51308, *T*<sup>3</sup> = 10, *T<sup>w</sup>* = 1, *R*<sup>1</sup> = *R*<sup>2</sup> = 2.4 Hz/p.u. MW, *Pr*<sup>1</sup> = *Pr*<sup>2</sup> = 1200 MW, *Kp*<sup>1</sup> = *Kp*<sup>2</sup> = 0.120 Hz/p.u. MW, *D*<sup>1</sup> = *D*<sup>2</sup> = 0.00833 p.u. MW/Hz.

Hybrid gravitational firefly algorithm parameters: number of objects = 70, number of iterations = 15, α = 20, *G*<sup>0</sup> = 100, β<sup>0</sup> = 0.2, and γ = 1.

GSA parameters: number of populations = 70, number of iterations = 15, G<sup>0</sup> = 100, and α = 20.

FA parameters: number of fireflies = 70, number of iterations = 15, β<sup>0</sup> = 0.2, γ = 1, and α = 0.5.

#### **References**

