Fractional-Order Load Frequency Control of an Interconnected Power System with a Hydrogen Energy-Storage Unit

Abstract

Modern power systems are confronted with widespread concern on the frequency stability issue due to the widespread integration of randomly fluctuating renewable resources. To address the above concern, this work introduces a load-frequency-control (LFC) scheme based on a parameter tuning strategy for fractional-order proportional–integral–derivative (FOPID) controller. Firstly, a two-area interconnected power system (IPS) model, including thermal, hydro, solar, wind, and gas power generator and a hydrogen-based energy-storage unit, is established. Then, a FOPID controller is designed for this IPS model, and an improved gradient-based optimizer (IGBO) is developed to adaptively regulate the parameters of the FOPID controllers. Finally, the effectiveness of the offered LFC scheme is tested through load disturbance and renewable energy fluctuations test scenarios and provides a comparison and robustness analysis among different schemes. The test results validated that the offered LFC scheme can effectively suppress the frequency fluctuations of the IPS and has excellent robustness.

1. Introduction

Development and application of renewable energy generation technology could be one critical enabler of the global transition to low-carbon emissions economies [1]. Nevertheless, solar energy, wind energy, tidal energy, and other renewable energies have random and fluctuating characteristics, and renewable energy generators have almost no inertia contribution to power grids since these generators are isolated to the grids by power electronic converters [2,3]. As a result, the IPSs with a substantial integration of renewable energy generators always face a frequency stability issue, and a suitable design of the LFC method is commonly required for these power systems.

Over the past few years, various intelligent control technologies such as model predictive control [4], sliding mode control [5], reinforcement learning [6], active disturbance rejection control [7], and fuzzy logic control [8] have been utilized in the design of load frequency control (LFC). The LFC scheme developed by the above technologies are usually effective in different test scenarios, but the major challenge of employing these methods is their heavy reliance on model accuracy and engineering feasibility, while the traditional proportional–integral–derivative (PID) controller-based LFC design, in fact, remains the most extensively used approach in practice due to its advantages of simple structure and low engineering implementation difficulty [9,10]. However, the interconnected structure and complicated module expressions of modern IPSs make the empirical parameter tuning methods time consuming and ineffective when using a PID controller-based LFC scheme. In addition, by using traditional PID controller-based LFC design and empirical parameter tuning methods, the robustness and dynamic properties of the system may also become worse in highly volatile, renewable-energy power generation facilities and highly flexible loads.

To address the above concerns, various parameter tuning and optimization techniques were developed for optimizing the PID controller-based LFC. These techniques include the whale optimization algorithm (WOA) [11], fuzzy rule (FL) [12], artificial neural network (ANN) [13], and gray wolf optimizer (GWO) [14]. According to these works, one can find that by introducing optimization strategies to regulate the parameters, the frequency stability of IPSs can be improved significantly. Nonetheless, traditional PID controllers only have three adjustable parameters, their parameter settings are limited, and their performance deteriorate with system nonlinearities and uncertainties. Therefore, the FOPID controller is proposed. The FOPID controller introduces fractional-order calculus operations to the traditional PID controller, which not only retains the simple structure of the traditional PID controller but also increases the adjustable dimensions of the controller parameters, thereby improving the disturbance rejection capability and robustness of the system in response to variations in the controlled system parameters [15,16,17]. Recently, such controllers have also been used for the LFC design. For example, in [18,19], a FOPID controller is utilized to mitigate the frequency fluctuations in multi-area IPSs. Moreover, the gases Brownian motion optimization (GBMO) and chaos game optimization (CGO) are exploited to tune the parameters of the controller, yielding satisfactory results. However, the investigated multi-area IPSs only involved conventional generating units, without taking into account renewable energy generation facilities. Studies [20,21] used a FOPID controller-based LFC scheme to enhance frequency regulation performance of the IPS with renewable energy power plants, and the controller parameters are tuned using Jellyfish Search Optimizer (JSO) and movable damped wave algorithm (MDWA), while renewable energy units were taken into account in their IPSs, the stochastic characteristics of the renewable energy sources were not considered. Study [22] proposed an LFC scheme based on improved particle swarm optimization (IPSO) to tune the parameters of a FOPID controller for a distributed generation system considering the randomness of renewable energy sources. The results of these studies confirm that in the application of LFC in IPSs, FOPID controllers exhibit superior dynamic performance compared to traditional PID controllers. However, previous studies have rarely considered the randomness of the renewable energy output when designing multi-area IPSs.

In addition, with the continuous construction of large-scale renewable energy generation facilities for various types of power sources, the integration of energy storage facilities has emerged as a prevalent supplementary method to enhance the frequency characteristics of IPSs with renewable resources [23,24]. Among various energy storage technologies such as pumped storage, flywheel energy storage, and electrochemical energy storage, hydrogen energy storage exhibits significant advantages in energy storage capacity and economic feasibility [25,26,27]. Over the last decade, various research groups have explored the utilization of hydrogen energy storage units in contributing to the frequency regulation of conventional power systems primarily powered by thermal power plants [28,29]. Building upon existing research, this study enhances the frequency response capability of an IPS considering the stochasticity of renewable energy sources by integrating hydrogen-based energy storage units.

In recent years, the gradient-based optimizer (GBO) has emerged as a promising solution for addressing optimization problems across diverse fields, primarily due to its simple yet powerful structure, ease of implementation, robust exploration capability, and ability to effectively evade local optima [30]. However, the GBO has some limitations when tackling large-scale and complex engineering problems, such as slow convergence speed and low precision in the later stages of the algorithm [31]. To address these challenges, this study proposes an improved GBO (IGBO) for obtaining optimal parameters for the FOPID controller in LFC issues within a two-area IPS integrated with hydrogen energy storage. Firstly, this study proposes a two-area IPS that considers the randomness of renewable resource, incorporating hydro, thermal, wind, photovoltaic, and gas turbine units, and hydrogen energy storage units. Subsequently, an improved Oustaloup approximation algorithm is employed to design FOPID controllers for the LFC system in the two-area IPS. Furthermore, an IGBO algorithm was developed to adaptively adjust the parameters of the FOPID controllers. Finally, a series of test scenarios is conducted to assess the effectiveness and superiority of the provided LFC scheme, and the robustness of the system was analyzed.

2. LFC Model of a Two-Area IPS

In general, modern power systems are formed by multiple interconnecting areas, with each area consisting of a variety of different types of power sources and load units. In this study, we focus on an IPS model consisting of two areas, which is depicted in Figure 1.

In the two-area IPS, the nominal frequency is 50 Hz, and the rated power of both area 1 and area 2 is 2000 MW. Area 1 comprises a non-reheat thermal plant, hydro power generator, solar photovoltaic cell, and wind power generator, while Area 2 consists of a thermal plant, non-reheat thermal plant, gas power generator, and wind power generator. Furthermore, hydrogen storage units are incorporated in both areas to enhance system performance. The mathematical models for each unit can be described as follows.

2.1. Mathematical Model of Thermal Power Plant

In a thermal power plant, the primary components are a speed governor and a steam turbine. The transfer function for the governor and steam turbine can be governed as

$$G_{\mathrm{G}}=\frac{1}{1+s{T}_{\mathrm{G}}}$$

(1)

$$G_{\mathrm{T}}=\frac{1}{1+s{T}_{\mathrm{T}}}$$

(2)

where $T_{\mathrm{G}}$ and $T_{\mathrm{T}}$ denote the time constants of the governor and turbine, respectively.

The transfer function of a reheater can be governed by

$$G_{\mathrm{T}}=\frac{1+s{K}_{\mathrm{r}}{T}_{\mathrm{r}}}{1+s{T}_{\mathrm{r}}}$$

(3)

where $K_{\mathrm{r}}$ and $T_{\mathrm{r}}$ are the gain and time constants of the reheater, respectively.

2.2. Mathematical Model of Hydro Power Generator

The hydro power generator comprises a hydraulic turbine and a hydraulic governor, with its transfer function represented as

$$G_{\mathrm{HT}}=\frac{1-s{T}_{\mathrm{W}}}{1+0.5s{T}_{\mathrm{W}}}$$

(4)

$$G_{\mathrm{HG}}=\frac{1+s{T}_{\mathrm{R}}}{\left(1+s{T}_{\mathrm{GH}}\right)\left(1+s(R_{\mathrm{T}}/R)T_{\mathrm{R}}\right)}$$

(5)

where $T_{\mathrm{W}}$ and $T_{\mathrm{GH}}$ represent the time constants of the turbine and the governor, respectively. $T_{\mathrm{R}}$ represents the reset time constant, R denotes the adjustment coefficient of hydro power generator, and $R_{\mathrm{T}}$ is the transient fall rate.

2.3. Mathematical Model of Solar Photovoltaic Cell

The power output of a solar photovoltaic cell is [22]

$$P_{\mathrm{PV}}=\eta S\mathsf{\Phi }\left\{1-0.005\left(T_{\mathrm{a}}+25\right)\right\}$$

(6)

where the coefficient $\eta$ represents the conversion efficiency of the photovoltaic unit. The coefficient S represents the surface area of the photovoltaic array, $T_{\mathrm{a}}$ denotes the working temperature, and $\mathsf{\Phi }$ denotes the solar radiation on the photovoltaic array, which can be expressed as

$$\mathsf{\Phi }=0.054H\left(t\right)+0.052H\left(t-18\right)+0.028H\left(t-37\right)-0.047H\left(t-58\right)-0.042H\left(t-74\right)+{\mathsf{\Phi }}_{\mathrm{n}}\left(t\right)$$

(7)

where $H\left(t\right)$ represents the step response, simulating sudden fluctuations, ${\mathsf{\Phi }}_{\mathrm{n}}\left(t\right)$ falls within the range [−0.015, 0.015], utilized to simulate stochastic aberrations. The mathematical model of photovoltaic system is

$$G_{\mathrm{PV}}=\frac{K_{\mathrm{PV}}}{1+s{T}_{\mathrm{PV}}}$$

(8)

where $K_{\mathrm{PV}}$ and $T_{\mathrm{PV}}$ are the gain and time constants of the photovoltaic unit, respectively.

2.4. Mathematical Model of Wind Power Generator

The modeling of wind power generators consists of two parts: a nonlinear model and a transfer function model. The former, as described in [22], does not contribute to frequency regulation. The power output of a wind power generator is

$$P_{\mathrm{WTG}}=\frac{1}{2}\rho A_{\mathrm{r}}{C}_{\mathrm{p}}{V}_{\mathrm{W}}^3$$

(9)

where the coefficient $\rho$ refers to air density, and $A_{\mathrm{r}}$ represents the swept area of blades. The coefficient $C_{\mathrm{p}}$ denotes the power coefficient and $V_{\mathrm{W}}$ denotes wind speed, which can be expressed as

$$V_{\mathrm{W}}=V_{\mathrm{WB}}+V_{\mathrm{WN}}$$

(10)

$$V_{\mathrm{WB}}=7.2H\left(t\right)-1.4H\left(t-15\right)-1.1H\left(t-33\right)+1.6H\left(t-57\right)+1.4H\left(t-76\right)$$

(11)

$$V_{\mathrm{WN}}=2{\sigma }^2{\displaystyle \sum _{i=1}^{50}}\sqrt{V_{\mathrm{V}}\left({\omega }_i\right)\Delta \omega }cos\left({\omega }_i-{\phi }_i\right)$$

(12)

where $V_{\mathrm{WB}}$ is the base fluctuation to simulate the base wind component. $V_{\mathrm{WN}}$ is small stochastic component, utilized to simulate noise wind component. ${\sigma }^2=200$ represents the variance of the noise component, and $V_{\mathrm{V}}\left({\omega }_i\right)$ represents the spectral density function. $\Delta \omega$ represents the frequency step, and ${\phi }_i$ belongs to the range $[0,2\pi ]$. The dynamic model of the wind power generator is

$$G_{\mathrm{W}}=\frac{K_{\mathrm{WTG}}}{1+s{T}_{\mathrm{WTG}}}$$

(13)

where $K_{\mathrm{WTG}}$ and $T_{\mathrm{WTG}}$ represent the gain and time constants of the wind power generator, respectively.

The transfer function model focuses on the active standby of wind turbines and takes wind power as an active resource for frequency modulation [5], which can be mathematically expressed as

$$G_{\mathrm{wind}}=\frac{1}{1+s{T}_{\mathrm{pw}2}}\frac{K_{\mathrm{pw}1}(1+s{T}_{\mathrm{pw}1})}{1+s}\frac{K_{\mathrm{pw}2}}{1+s}$$

(14)

where $T_{\mathrm{pw}1}$, $T_{\mathrm{pw}2}$, $K_{\mathrm{pw}1}$, and $K_{\mathrm{pw}2}$ denote the time and gain constants of the wind power governor, respectively.

2.5. Mathematical Model of Gas Power Generator

The mathematical model of gas power generator can be stated as [32]

$$G_{\mathrm{GT}}=\frac{1}{c_{\mathrm{g}}+s{b}_{\mathrm{g}}}\frac{1+s{X}_{\mathrm{c}}}{1+s{Y}_{\mathrm{c}}}\frac{1+s{T}_{\mathrm{cr}}}{1+s{T}_{\mathrm{f}}}\frac{1}{1+s{T}_{\mathrm{cd}}}$$

(15)

where $c_{\mathrm{g}}$ and $b_{\mathrm{g}}$ represent the gain and time constants of the gas turbine, respectively. $X_{\mathrm{c}}$ and $Y_{\mathrm{c}}$ are complex parameters that reflect the impact of the gas turbine on the phase and amplitude of the input signal, $T_{\mathrm{cr}}$, $T_{\mathrm{f}}$, and $T_{\mathrm{cd}}$ are the time constants of the cooling system, fuel injection and entire gas power generator, respectively.

2.6. Mathematical Model of Hydrogen Storage Unit

The hydrogen storage system comprises an electrolyzer, a storage tank, and a hydrogen fuel cell. When surplus power is generated by the new energy generators within the system, the electrolyzer produces hydrogen, which is subsequently compressed and stored in the hydrogen storage tanks. During periods of heightened load demand and electricity scarcity, the hydrogen fuel cell converts the chemical energy into electrical power to meet the system’s load requirements. The mathematical model of the electrolyzer and the hydrogen fuel cell in the system can be represented as [28]

$$G_{\mathrm{HAE}}=\frac{K_{\mathrm{HAE}}}{1+s{T}_{\mathrm{HAE}}}$$

(16)

$$G_{\mathrm{FC}}=\frac{K_{\mathrm{FC}}}{1+s{T}_{\mathrm{FC}}}$$

(17)

where $K_{\mathrm{HAE}}$ and $K_{\mathrm{FC}}$ are the gain constants of the electrolyzer and the hydrogen fuel cell, respectively. $T_{\mathrm{HAE}}$ and $T_{\mathrm{FC}}$ represent the time constants of the electrolyzer and the hydrogen fuel cell, respectively.

Based on the mathematical models of each unit, as well as the hydrogen storage unit, an LFC model is constructed for the IPS consisting of two areas, as depicted in Figure 2. In Figure 2, the term B is the frequency offset coefficient, ${\mathrm{ACE}}_1$, and ${\mathrm{ACE}}_2$ represent the area control errors (ACE) for area 1 and area 2, respectively, and $R_{\mathrm{T}}$ and $R_{\mathrm{H}}$ represent the governor constants of the thermal plants and the hydro power generator, respectively. $R_{\mathrm{G}}$ and $R_{\mathrm{W}}$ represent the governor parameters of the gas power generator and wind power generator, respectively. ${\alpha }_{\mathrm{G}}$, ${\alpha }_{\mathrm{T}}$, ${\alpha }_{\mathrm{N}}$, ${\alpha }_{\mathrm{W}}$ and ${\alpha }_{\mathrm{H}}$ represent the power distribution ratio of gas, thermal, non-reheat thermal, wind power, and hydro power generators, respectively, $\Delta P_{\mathrm{L}1}$ and $\Delta P_{\mathrm{L}2}$ represent load disturbances in area 1 and area 2, respectively. $K_{\mathrm{PS}}$ and $T_{\mathrm{PS}}$ denote the gain and time constants of the IPS, respectively. $T_{12}$ denotes the synchronization coefficient of the tie line, and $\Delta f_1$ and $\Delta f_2$ represent the frequency deviation of the area 1 and area 2, respectively. $\Delta P_{12}$ denotes the power flow of the tie line.

3. Design of FOPID Controller-Based LFC

The LFC scheme in this study is based on a FOPID controller, which is two more parameters compared to the traditional PID controller. Basically, the transfer function of FOPID controller can be expressed by

$$G_{\mathrm{C}}\left(s\right)=K_{\mathrm{P}}+\frac{K_{\mathrm{I}}}{s^{\lambda }}+K_{\mathrm{D}}{s}^{\mu }$$

(18)

where $K_{\mathrm{P}}$, $K_{\mathrm{I}}$ and $K_{\mathrm{D}}$ represent the proportional coefficient, integral coefficient and differential coefficient, respectively. Additionally, $1/s^{\lambda }$ and $s^{\mu }$ denote the fractional integral operator and fractional differential operator, respectively. These ideal fractional-order operators cannot be implemented directly by integer-order rational fractions. Thus, French scholar Professor Oustaloup proposed a piecewise approximation method based on the Oustaloup filter, where all segments are composed of integer-order zeros and poles, resulting in the asymptotic magnitude-frequency response alternating between 0 and −20 dB/dec. The standard form of the Oustaloup filter is as follows:

$$G\left(s\right)={\omega }_{\mathrm{h}}^{\alpha }{\displaystyle \prod _{k=1}^{N_{\mathrm{f}}}}\frac{s+{\omega }_k^{\prime }}{s+{\omega }_k}$$

(19)

$${\omega }_k^{\prime }={\omega }_{\mathrm{b}}{\left(\sqrt{\frac{{\omega }_{\mathrm{h}}}{{\omega }_{\mathrm{b}}}}\right)}^{\frac{2k-1-\alpha }{N_{\mathrm{f}}}}$$

(20)

$${\omega }_k={\omega }_{\mathrm{b}}{\left(\sqrt{\frac{{\omega }_{\mathrm{h}}}{{\omega }_{\mathrm{b}}}}\right)}^{\frac{2k-1+\alpha }{N_{\mathrm{f}}}}$$

(21)

where $\alpha$ represents the order of the fractional-order operator, and $N_{\mathrm{f}}$ represents the order of the filter. ${\omega }_{\mathrm{b}}$ and ${\omega }_{\mathrm{h}}$ represent the lower and upper limits of the predetermined frequency range. ${\omega }_k^{\prime }$ and ${\omega }_k$ represent the constructed zeros and poles of the transfer function, respectively. However, the approximation performance of the Oustaloup filter is unsatisfactory at the boundary points of the selected frequency range, and the order of the numerator is the same as that of the denominator. Therefore, this study has improved the above filter, and the mathematical model of the improved Oustaloup filter is as follows [33]:

$$G\left(s\right)={\left(\frac{d{\omega }_{\mathrm{b}}}{b}\right)}^{\alpha }\frac{d{s}^2+b{\omega }_{\mathrm{h}}s}{d(1-\alpha )s^2+b{\omega }_{\mathrm{h}}s+d\alpha }{\displaystyle \prod _{k=-N_{\mathrm{f}}}^{N_{\mathrm{f}}}}\frac{s+{\omega }_k^{\prime }}{s+{\omega }_k}$$

(22)

$${\omega }_k^{\prime }={\omega }_{\mathrm{b}}{\left(\sqrt{\frac{{\omega }_{\mathrm{h}}}{{\omega }_{\mathrm{b}}}}\right)}^{\frac{k+N_{\mathrm{f}}+0.5(1-\alpha )}{2N_{\mathrm{f}}+1}}$$

(23)

$${\omega }_k={\omega }_{\mathrm{b}}{\left(\sqrt{\frac{{\omega }_{\mathrm{h}}}{{\omega }_{\mathrm{b}}}}\right)}^{\frac{k+N_{\mathrm{f}}+0.5(1+\alpha )}{2N_{\mathrm{f}}+1}}$$

(24)

where the coefficients b and d are preset coefficients with values of 9 and 10, respectively. The predetermined frequency band is set to [0.001 Hz, 1000 Hz], and $N_{\mathrm{f}}$ is set to 5. The comparison between the ideal fractional-order operators and approximate fractional-order operators is shown in Figure 3.

As depicted in Figure 3, it can be observed that within the predetermined frequency range, the amplitude–frequency characteristic curves of the approximate fractional-order operator and the ideal fractional-order operator coincide with each other. Additionally, the phase-frequency characteristic curves of both operators fall within the preset phase error range throughout the predetermined frequency range. Obviously, the fitting effect of the improved Oustaloup-based approximation method at the boundary points of the predetermined frequency band is better than that of the original Oustaloup-based approximation method. The results of the comparison confirm that the approximate fractional-order operators can be applied.

4. Parameter Tuning Process Based on IGBO

To ensure that the two-area IPS in Figure 2 can swiftly recover frequency stability under load or renewable resources fluctuations, this paper proposes a parameter tuning strategy by an IGBO algorithm, and the objective fitness function is formulated using the integral of time multiplied absolute error (ITAE) criterion based on the ACE, which can be written as

$$min{J}_{\mathrm{ITAE}}={\int }_0^{T_{\mathrm{sim}}}t\left(\left|\Delta f_1\right|+\left|\Delta f_2\right|+\left|\Delta P_{12}\right|\right)\mathrm{d}t$$

(25)

where $T_{\mathrm{sim}}$ is the total duration of the optimization process. With the objective of minimizing $J_{\mathrm{ITAE}}$, the IGBO algorithm is utilized to tune the parameters of the FOPID controller in the IPS, and the LFC controller structure based on IGBO optimizer is shown in Figure 4.

The proposed IGBO algorithm integrates the gradient search mechanism (GSM) based on the sorting algorithm and the improved local escape operator (ILEO) based on chaotic mapping during the optimization process [34]. The specific updating method of IGBO is as follows.

4.1. Initialization

The IGBO algorithm initializes the parameters of the FOPID controller based on

$$X_n=Low+rand\times (Up-Low)$$

(26)

where $n=1,2,\dots ,N$, N corresponds to the total number of IGBO vectors. $Low$ and $Up$ denote the lower and upper bounds of the parameter space. $rand$ denotes a randomly generated number within the range of 0 to 1.

4.2. Parameter Update Method Based on GSM

GSM is used for the global search in the parameter space, and its definition is given as

$$\mathrm{GSM}=randn\times {\rho }_1\times \frac{2\Delta X\times X_n}{y{p}_{\mathrm{n}}-y{q}_{\mathrm{n}}+\epsilon }$$

(27)

where $randn$ represents a random number that follows a normal distribution, and $\epsilon$ represents a small number in the range of 0 to 0.1. ${\rho }_1$ is a crucial parameter that determines the global exploration ability of the IGBO algorithm. To obtain a better solution for the FOPID controller parameters, a sorting algorithm is proposed to update the parameter ${\rho }_1$, which can be represented as

$${\rho }_1=\left\{\begin{array}{c}\begin{array}{cc}\hfill r1/N+0.15\times randn& ifrand<0.5\hfill \\ \hfill n/N+0.15\times randn& else\hfill \end{array}\hfill \end{array}\right.$$

(28)

where $r1$ is an integer randomly selected from 1 to N. Furthermore, the GSM definition is represented by Equation (29) for $\Delta X$, and $y{p}_{\mathrm{n}}$ and $y{q}_{\mathrm{n}}$ are represented by Equation (32) and Equation (33), respectively.

$$\Delta X=rand(1:D)\times \left|step\right|$$

(29)

$$step=\frac{(X_{\mathrm{best}}-X_{r1}^{it})+\delta }{2}$$

(30)

$$\delta =2\times rand\times \left(\left|\frac{X_{r1}^{it}+X_{r2}^{it}+X_{r3}^{it}+X_{r4}^{it}}{4}-X_n^{it}\right|\right)$$

(31)

where $rand(1:D)$ is a D-dimensional random vector, and D corresponds to the number of parameters of the FOPID controller to be optimized. $step$ is determined by the optimal solution ($X_{\mathrm{best}}$) and the random selection vector ($X_{r1}^{it}$) in the iterations. $it$ is the current number of iterations within the range 1 to $I{t}_{\max }$, and $I{t}_{\max }$ represents the maximum number of iterations. The variables $r2$, $r3$, and $r4$ represent distinct integers that are randomly selected from 1 to N.

$$y{p}_{\mathrm{n}}=rand\times \left(\frac{\left[Z_{n+1}+X_n\right]}{2}+rand\times \Delta X\right)$$

(32)

$$y{q}_{\mathrm{n}}=rand\times \left(\frac{\left[Z_{n+1}+X_n\right]}{2}-rand\times \Delta X\right)$$

(33)

$$Z_{n+1}=X_n-randn\times \frac{2\Delta X\times X_n}{X_{\mathrm{worst}}-X_{\mathrm{best}}+\epsilon }+DM$$

(34)

$$DM=rand\times {\rho }_2\times \left(X_{\mathrm{best}}-X_n^{it}\right)$$

(35)

where $X_{\mathrm{worst}}$ represents the worst vector during the iterations. $DM$ is the direction of movement, guiding the individual vector towards a more optimal solution. The update method for parameter ${\rho }_2$ is the same as ${\rho }_1$.

According to GSM, the IGBO generates two solutions with different capabilities using Equation (36) and Equation (37), respectively.

$$X{1}_n^{it}=X_n^{it}-randn\times {\rho }_1\times \frac{2\Delta X\times X_n^{it}}{y{p}_{\mathrm{n}}-y{q}_{\mathrm{n}}+\epsilon }+rand\times {\rho }_2\times \left(X_{\mathrm{best}}-X_n^{it}\right)$$

(36)

$$X{2}_n^{it}=X_{\mathrm{best}}-randn\times {\rho }_1\times \frac{2\Delta X\times X_n^{it}}{y{p}_{\mathrm{n}}-y{q}_{\mathrm{n}}+\epsilon }+rand\times {\rho }_2\times \left(X_{r1}^{it}-X_{r2}^{it}\right)$$

(37)

where Equation (36) is beneficial for the global search, and Equation (37) is beneficial for the local exploitation. Therefore, based on the solutions $X{1}_n^{it}$ and $X{2}_n^{it}$, the IGBO generates the next solution ($X_n^{it+1}$) which is defined as

$$X_{n,j}^{it+1}=\left\{\begin{array}{c}\left\{\begin{array}{c}X{1}_{n,j}^{it}ifran{d}_2<0.5\hfill \\ X{2}_{n,j}^{it}ifran{d}_2\ge 0.5\hfill \end{array}\right.ifran{d}_1<P_{cr}\hfill \\ X_{n,j}^{it}ifran{d}_1\ge P_{cr}\hfill \end{array}\right.j=1toD$$

(38)

where $ran{d}_1$ and $ran{d}_2$ are two different random numbers, $P_{cr}$ is the crossover probability, and its update method is the same as ${\rho }_1$.

4.3. ILEO

The ILEO is used in IGBO to avoid algorithms falling into local optimality. The generation method of the ILEO vector solution can be expressed as

$$X_{\mathrm{ILEO}}^{it+1}=\left\{\begin{array}{c}{X}_n^{it+1}+{\rho }_1\left(X_{\mathrm{best}2}+X_k^{it}\right)+{\rho }_2(X_{r1}^{it}-X_{r2}^{it})ifrand<\frac{(1-it)/I{t}_{max}}{2}\hfill \\ X_{\mathrm{best}}+{\rho }_1\left(X_{\mathrm{best}2}+X_k^{it}\right)+{\rho }_2(X_{r1}^{it}-X_{r2}^{it})else\hfill \end{array}\right.$$

(39)

where $X_{\mathrm{best}2}$ is the second best solution. $X_k^{it}$ represents a solution that is randomly selected from 1 to N. Therefore, the next generation solutions generated by the IGBO can be represented as

$$X_n^{it+1}=\left\{\begin{array}{c}{X}_{ILEO}^{it+1}ifrand<LC\hfill \\ X_n^{it+1}else\hfill \end{array}\right.$$

(40)

where the term $LC$ denotes a chaotic logistic map, and its update method is $LC=4LC(1-LC)$, and the initial LC is 0.67. The IGBO optimizes the parameters of two FOPID controllers in the following parameter space.

$$\left\{\begin{array}{c}{K}_{\mathrm{P}i}^{Low}\le K_{\mathrm{P}i}\le K_{\mathrm{P}i}^{Up}\hfill \\ K_{\mathrm{I}i}^{Low}\le K_{\mathrm{I}i}\le K_{\mathrm{I}i}^{Up}\hfill \\ K_{\mathrm{D}i}^{Low}\le K_{\mathrm{D}i}\le K_{\mathrm{D}i}^{Up}\hfill \\ {\lambda }_i^{Low}\le {\lambda }_i\le {\lambda }_i^{Up}\hfill \\ {\mu }_i^{Low}\le {\mu }_i\le {\mu }_i^{Up}\hfill \end{array}\right.$$

(41)

where $i=1,2$ indicates area 1 and area 2, respectively. The parameter tuning process based on IGBO is depicted in Figure 5.

5. Validation and Comparison

5.1. Performance Testing of IGBO

To validate the search capability of the proposed IGBO algorithm, PSO, WOA, MDWA, CGO, JSO, GBO, and IGBO were employed to minimize three congresses on evolutionary computation (CEC) benchmark functions ($F1\left(x\right)$, $F2\left(x\right)$, and $F3\left(x\right)$), separately. The simulations run on 64-bit Windows, which has 4 GB of memory. Set the number of iterations for all algorithms to 500. The test results are presented in

Table 1. Benchmark test function values for different algorithms.

Algorithm	$\mathit{F}1\left(\mathit{x}\right)$	$\mathit{F}2\left(\mathit{x}\right)$	$\mathit{F}3\left(\mathit{x}\right)$
PSO	1.52 $\times {10}^{-17}$	4.44 $\times {10}^{-15}$	−10.5364
WOA	0.05	4.44 $\times {10}^{-15}$	−10.87
MDWA	5.33 $\times {10}^{-13}$	4.44 $\times {10}^{-15}$	−3.84
CGO	8.54 $\times {10}^{-47}$	7.99 $\times {10}^{-15}$	−10.5359
JSO	1.71 $\times {10}^{-24}$	7.99 $\times {10}^{-15}$	−10.5364
GBO	4.37 $\times {10}^{-46}$	8.44 $\times {10}^{-16}$	−10.5364
IGBO	8.99 $\times {10}^{-73}$	8.44 $\times {10}^{-16}$	−10.5364

. The convergence curves of the three benchmark functions under different algorithms are provided in Figure 6.

$$F1\left(x\right)={\displaystyle \sum _{i=1}^n}{\left({\displaystyle \sum _{i=1}^n}{x}_i\right)}^2$$

(42)

$$F2\left(x\right)=-20exp\left(-0.2\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n}{x}_i^2}\right)-exp\left(-\frac{1}{n}{\displaystyle \sum _{i=1}^n}cos\left(2\pi x_i\right)\right)+20+e$$

(43)

$$F3\left(x\right)=-{\displaystyle \sum _{i=1}^{10}}{\left[C_i+{\left(x-A_i\right)}^T\left(x-A_i\right)\right]}^{-1}$$

(44)

where $F1\left(x\right)$ is a unimodal standard function used to measure the exploitability of an algorithm. $F2\left(x\right)$ is a multimodal function used to test exploration performance. $F3\left(x\right)$ is a fixed-dimensional function used to demonstrate the ability to explore in low dimensions.

Based on Table 1 and Figure 6, it can be observed that for the unimodal function $F1\left(x\right)$, the IGBO algorithm achieves a better minimum value compared to other algorithms. For multimodal function $F2\left(x\right)$, both IGBO and GBO obtain the same minimum value, but the convergence speed of IGBO is faster than that of GBO. For fixed-dimensional function $F3\left(x\right)$, IGBO converges to the minimum value of the benchmark function earliest. The test results indicate that the IGBO algorithm possesses high optimization capability and speed.

To test the effectiveness of the offered IGBO-FOPID-based LFC scheme in this study, a two-area IPS LFC model is constructed according to Figure 2 in this section, and the system parameter settings are depicted in Appendix A. The control performance of the IGBO-FOPID-based LFC scheme is compared with existing LFC methods, including the FOPID controller regulated through WOA (WOA-FOPID), MDWA (MDWA-FOPID), PSO (PSO-FOPID), CGO (CGO-FOPID), JSO(JSO-FOPID), and GBO tuned PID/FOPID controllers (GBO-PID/FOPID). The maximum number of iterations, $I{t}_{\max }$, for each control scheme algorithm is set to 30, and N is set to 50.

Table 2. Parameters of two FOPID controllers obtained by different LFC schemes.

LFC Scheme	${\mathit{K}}_{\mathbf{P}1}$	${\mathit{K}}_{\mathbf{I}1}$	${\mathit{K}}_{\mathbf{D}1}$	${\mathit{\lambda }}_1$	${\mathit{\mu }}_1$	${\mathit{K}}_{\mathbf{P}2}$	${\mathit{K}}_{\mathbf{I}2}$	${\mathit{K}}_{\mathbf{D}2}$	${\mathit{\lambda }}_2$	${\mathit{\mu }}_2$	Execution Time (s)
GBO-PID	68.22	385.33	25.42	-	-	61.77	430.56	16.58	-	-	3327.58
PSO-FOPID	583.84	513.51	31.38	0.97	1.19	461.51	512.46	22.13	1.15	1.21	3425.62
MDWA-FOPID	516.76	327.85	4.89	2	1.81	622.30	25.11	3.43	2	1.89	3404.82
WOA-FOPID	54.38	308.42	10.66	0.41	1.19	305.22	525.31	12.62	1.17	1.32	3583.06
CGO-FOPID	589.18	507.01	28.50	0.99	1.22	477.50	504.34	24.60	1.14	1.19	5496.67
JSO-FOPID	404.94	609.39	23.76	1.06	1.31	400.06	672.12	21.70	1.13	1.45	3037.6
GBO-FOPID	385.42	519.77	15.98	1.09	1.23	695.18	683.66	16.16	1.12	1.24	3673.43
IGBO-FOPID	545.39	642.31	11.48	1.02	1.62	697.38	698.66	22.18	1.04	1.47	3176.52

lists the controller parameters obtained by different LFC schemes.

Based on the execution times of each optimization algorithm presented in Table 2, it is evident that the convergence speed of the IGBO has been enhanced. The LFC scheme in Table 2 is initially tested using step load-disturbance scenarios. Subsequently, its performance is validated under scenarios involving random fluctuations of photovoltaic and wind power. Furthermore, the robustness of the IPS is evaluated by varying the inertia constant.

5.2. Step Load-Disturbance Test Scenarios

In this scenario, the output of renewable energy units is not considered. At 0 s, 20 MW step load disturbances, i.e., 0.01 puMW are introduced to area 1 and area 2 at the same time. In simulation, the performances of eight different LFC schemes are compared, in which the GBO-PID, PSO-FOPID, MDWA-FOPID, WOA-FOPID, CGO-FOPID, JSO-FOPID, GBO-FOPID, and IGBO-FOPID-based LFC scheme are used. The response of the IPS under step load disturbance with different LFC schemes are depicted in Figure 7.

The results presented in Figure 7 demonstrate the superiority of the FOPID controller-based LFC scheme over the PID controller-based LFC scheme. Implementation of the FOPID controller leads to a considerable enhancement in system frequency modulation performance. Specifically, the FOPID controller-based LFC scheme effectively reduces the frequency deviation overshoot and minimizes the frequency fluctuation range within the system. Furthermore, the proposed IGBO-FOPID-based LFC scheme outperforms FOPID controllers optimized using PSO, MDWA, WOA, CGO, JSO, and GBO algorithms. Notably, the IGBO-FOPID-based LFC scheme yields lower overshoot, as well as shorter settling time to achieve system stability. These outcomes affirm that the IGBO-FOPID-based LFC scheme significantly enhances the frequency response capability of the IPS under step load-disturbance scenarios.

To further evaluate the superiority of the proposed LFC scheme,

Table 3. Maximum US/Maximum OS/ ISE/ITSE/IAE/ITAE of system response under step load-disturbance scenarios.

LFC Scheme	Maximum US			Maximum OS			ISE	ITSE	IAE	ITAE
	$\Delta {\mathbf{f}}_{\mathbf{1}}$	$\Delta {\mathbf{f}}_{\mathbf{2}}$	$\Delta {\mathbf{P}}_{\mathrm{tie}\mathbf{12}}$	$\Delta {\mathbf{f}}_{\mathbf{1}}$	$\Delta {\mathbf{f}}_{\mathbf{2}}$	$\Delta {\mathbf{P}}_{\mathrm{tie}\mathbf{12}}$
GBO-PID	−1.10 $\times {10}^{-3}$	−1.70 $\times {10}^{-3}$	−5.15 $\times {10}^{-5}$	4.97 $\times {10}^{-4}$	9.94 $\times {10}^{-4}$	1.42 $\times {10}^{-4}$	4.86 $\times {10}^{-7}$	1.10 $\times {10}^{-7}$	8.54 $\times {10}^{-4}$	4.27 $\times {10}^{-4}$
PSO-FOPID	−5.62 $\times {10}^{-4}$	−9.19 $\times {10}^{-4}$	−1.89 $\times {10}^{-5}$	2.95 $\times {10}^{-4}$	5.23 $\times {10}^{-4}$	4.74 $\times {10}^{-5}$	6.86 $\times {10}^{-8}$	7.33 $\times {10}^{-9}$	2.31 $\times {10}^{-4}$	8.28 $\times {10}^{-5}$
MDWA-FOPID	−3.71 $\times {10}^{-4}$	−4.86 $\times {10}^{-4}$	−2.22 $\times {10}^{-5}$	1.34 $\times {10}^{-4}$	2.64 $\times {10}^{-4}$	3.07 $\times {10}^{-5}$	4.80 $\times {10}^{-8}$	2.08 $\times {10}^{-8}$	4.02 $\times {10}^{-4}$	4.47 $\times {10}^{-4}$
WOA-FOPID	−1 $\times {10}^{-2}$	−9.55 $\times {10}^{-4}$	−8.37 $\times {10}^{-6}$	9.34 $\times {10}^{-5}$	3.45 $\times {10}^{-5}$	5.91 $\times {10}^{-6}$	7.79 $\times {10}^{-8}$	7.75 $\times {10}^{-9}$	2.82 $\times {10}^{-4}$	4.12 $\times {10}^{-4}$
CGO-FOPID	−6.02 $\times {10}^{-4}$	−8.54 $\times {10}^{-4}$	−2.44 $\times {10}^{-6}$	2.87 $\times {10}^{-5}$	4.29 $\times {10}^{-4}$	3.72 $\times {10}^{-5}$	4.75 $\times {10}^{-8}$	3.87 $\times {10}^{-9}$	1.83 $\times {10}^{-4}$	9.24 $\times {10}^{-5}$
JSO-FOPID	−5.27 $\times {10}^{-4}$	−7.35 $\times {10}^{-4}$	−9.39 $\times {10}^{-7}$	8.76 $\times {10}^{-5}$	1.10 $\times {10}^{-4}$	2.56 $\times {10}^{-5}$	2.45 $\times {10}^{-8}$	1.01 $\times {10}^{-9}$	9.41 $\times {10}^{-5}$	5.93 $\times {10}^{-5}$
GBO-FOPID	−7.24 $\times {10}^{-4}$	−9.76 $\times {10}^{-4}$	−3.85 $\times {10}^{-5}$	2.98 $\times {10}^{-4}$	8.64 $\times {10}^{-4}$	3.60 $\times {10}^{-5}$	1.20 $\times {10}^{-7}$	1.66 $\times {10}^{-8}$	3.05 $\times {10}^{-4}$	9.93 $\times {10}^{-5}$
IGBO-FOPID	−3.09 $\times {10}^{-4}$	−4.72 $\times {10}^{-4}$	−6.21 $\times {10}^{-7}$	3.21 $\times {10}^{-5}$	2.06 $\times {10}^{-5}$	9.51 $\times {10}^{-6}$	9.79 $\times {10}^{-9}$	6.87 $\times {10}^{-10}$	8.76 $\times {10}^{-5}$	4.89 $\times {10}^{-5}$

lists some key performance indicators of the IPS with different LFC schemes, such as the maximum undershoot (US), maximum overshoot (OS), integral-squared error (ISE), integral time-multiplied squared error (ITSE), integral absolute error (IAE), and ITAE of $\Delta f_1$/$\Delta f_2$/$\Delta P_{\mathrm{tie}12}$.

The comparison in Table 3 showcases that the IGBO-FOPID-based LFC scheme achieves superior values of maximum US, maximum OS, ISE, ITSE, IAE, and ITAE compared to other optimization strategy. This further illustrate that the offered IGBO-FOPID LFC scheme outperforms GBO-PID/PSO-FOPID/MDWA-FOPID/WOA-FOPID/CGO-FOPID/JSO-FOPID/GBO-FOPID-based LFC schemes.

5.3. Photovoltaic and Wind Power Random Fluctuation Test Scenarios

The effectiveness of the IGBO-FOPID-based LFC scheme is verified under the scenario of random power output from renewable energy sources. The output fluctuations of photovoltaic power and wind power, as described in Section 2.3 and Section 2.4, are utilized in this section. The fluctuations are visualized in Figure 8.

We keep the control parameters in Table 2 unchanged and use different LFC schemes for simulation. The response of IPS under random fluctuations in photovoltaic and wind power using different LFC schemes are shown in Figure 9. The maximum US/maximum OS/ISE/ITSE/IAE/ITAE of $\Delta f_1$/$\Delta f_2$/$\Delta P_{\mathrm{tie}12}$ results are displayed in

Table 4. Maximum US/Maximum OS/ ISE/ITSE/IAE/ITAE of system response under photovoltaic and wind power random fluctuation scenarios.

LFC Scheme	Maximum US			Maximum OS			ISE	ITSE	IAE	ITAE
	$\Delta {\mathbf{f}}_{\mathbf{1}}$	$\Delta {\mathbf{f}}_{\mathbf{2}}$	$\Delta {\mathbf{P}}_{\mathrm{tie}\mathbf{12}}$	$\Delta {\mathbf{f}}_{\mathbf{1}}$	$\Delta {\mathbf{f}}_{\mathbf{2}}$	$\Delta {\mathbf{P}}_{\mathrm{tie}\mathbf{12}}$
GBO-PID	−4.45 $\times {10}^{-4}$	−2.67 $\times {10}^{-4}$	−8.04 $\times {10}^{-5}$	4.49 $\times {10}^{-4}$	5.36 $\times {10}^{-4}$	5.84 $\times {10}^{-5}$	2.14 $\times {10}^{-6}$	9.95 $\times {10}^{-5}$	1.58 $\times {10}^{-2}$	0.77
PSO-FOPID	−1.83 $\times {10}^{-4}$	−7.30 $\times {10}^{-5}$	−2.45 $\times {10}^{-5}$	1.96 $\times {10}^{-4}$	1.68 $\times {10}^{-4}$	1.63 $\times {10}^{-5}$	3.90 $\times {10}^{-7}$	1.88 $\times {10}^{-5}$	6.58 $\times {10}^{-3}$	0.32
MDWA-FOPID	−2.61 $\times {10}^{-4}$	−2.59 $\times {10}^{-4}$	−9.84 $\times {10}^{-5}$	2.67 $\times {10}^{-4}$	2.49 $\times {10}^{-4}$	9.43 $\times {10}^{-5}$	3.35 $\times {10}^{-6}$	1.68 $\times {10}^{-4}$	2.62 $\times {10}^{-2}$	1.31
WOA-FOPID	−2.81 $\times {10}^{-4}$	−1.12 $\times {10}^{-4}$	−3.91 $\times {10}^{-5}$	3.11 $\times {10}^{-4}$	2.17 $\times {10}^{-4}$	4.82 $\times {10}^{-5}$	1.45 $\times {10}^{-6}$	5.87 $\times {10}^{-5}$	1.51 $\times {10}^{-2}$	0.69
CGO-FOPID	−1.74 $\times {10}^{-4}$	−7.02 $\times {10}^{-5}$	−2.25 $\times {10}^{-5}$	2.03 $\times {10}^{-4}$	1.48 $\times {10}^{-4}$	1.60 $\times {10}^{-5}$	2.56 $\times {10}^{-6}$	1.18 $\times {10}^{-5}$	5.37 $\times {10}^{-3}$	0.26
JSO-FOPID	−1.08 $\times {10}^{-4}$	−5.72 $\times {10}^{-5}$	−2.33 $\times {10}^{-5}$	1.54 $\times {10}^{-4}$	1.20 $\times {10}^{-4}$	1.50 $\times {10}^{-5}$	1.15 $\times {10}^{-7}$	4.80 $\times {10}^{-5}$	3.71 $\times {10}^{-3}$	0.17
GBO-FOPID	−2.13 $\times {10}^{-4}$	−8.77 $\times {10}^{-5}$	−1.94 $\times {10}^{-5}$	1.81 $\times {10}^{-4}$	1.87 $\times {10}^{-4}$	1.39 $\times {10}^{-5}$	4.99 $\times {10}^{-7}$	2.36 $\times {10}^{-5}$	8.12 $\times {10}^{-3}$	0.39
IGBO-FOPID	−5.27 $\times {10}^{-5}$	−4.07 $\times {10}^{-5}$	−1.55 $\times {10}^{-5}$	8.96 $\times {10}^{-5}$	8.25 $\times {10}^{-5}$	1.25 $\times {10}^{-5}$	5.25 $\times {10}^{-8}$	1.91 $\times {10}^{-6}$	2.55 $\times {10}^{-3}$	0.11

.

In Figure 9, one can find that the effectiveness of the offered IGBO-FOPID-based LFC scheme in regulating the IPS under scenarios involving random fluctuations in photovoltaic and wind power. Notably, the system incorporating the IGBO-FOPID-based LFC scheme exhibits smaller frequency fluctuations in two areas when compared to GBO-PID/PSO-FOPID/MDWA-FOPID/WOA-FOPID/CGO-FOPID/JSO-FOPID/GBO-FOPID-based LFC schemes. Additionally, the comparison presented in Table 4 reveals that the IPS performance with the IGBO-FOPID-based LFC scheme outperforms that of GBO-PID/PSO-FOPID/MDWA-FOPID/WOA-FOPID/CGO-FOPID/JSO-FOPID/GBO-FOPID-based LFC schemes in terms of maximum US, maximum OS, and various performance indicators. These findings highlight the superior anti-jamming capabilities of the proposed IGBO-FOPID-based LFC scheme compared to other LFC methods. Thus, the proposed LFC scheme is highly suitable for multi-source IPS applications incorporating renewable energy sources.

5.4. Robustness Analysis

Robustness refers to the system’s ability to withstand disturbances such as uncertainties and parameter variations. In engineering control, robustness is utilized to evaluate the control method’s capability to maintain stability in the face of various uncertainties. To evaluate the robustness of the proposed IGBO-FOPID-based LFC scheme in this study, the inertia constant of the IPS model was reduced by 25% and 50% from its nominal value. The parameters of the system after the change in inertia are provided in Appendix B. The controller parameters in Table 2 remained unchanged, and GBO-PID, GBO-FOPID, and IGBO-FOPID-based LFC schemes were tested, with the results presented in Figure 10, Figure 11, and Figure 12, respectively.

From Figure 10, it can be observed that when the inertia constant of the IPS is reduced by 25%, the system response obtained from the LFC scheme based on the traditional PID controller can still maintain stability, albeit barely. However, when the inertia constant is reduced by 50%, the system response obtained from the GBO-PID-based LFC scheme becomes unstable. This indicates that the traditional PID controller-based LFC scheme is not robust when faced with significant changes in the parameters of the IPS.

From Figure 11 and Figure 12, one can see that the system response obtained from the LFC scheme based on the FOPID controller remains consistent before and after the change in inertia constant. This illustrated that the FOPID controller-based LFC scheme has excellent robustness when facing significant changes in the parameters of the IPS. Additionally, the proposed IGBO-FOPID-based LFC scheme exhibits better frequency stability compared to the GBO-FOPID-based LFC scheme, further demonstrating that the proposed IGBO-FOPID-based LFC scheme can enhance the system’s robustness.

6. Conclusions

This study presents a novel LFC scheme based on the IGBO-FOPID controller to address frequency fluctuation issues in an IPS involving renewable energy generation units. The LFC scheme is tested under two scenarios: step load-disturbance test scenarios and renewable energy sources random fluctuation test scenarios. The test results indicate that the IPS response obtained by the proposed LFC scheme outperforms the existing literature’s PSO/WOA/MDWA/CGO/JSO-optimized FOPID and GBO-optimized PID/FOPID controllers. Compared to the aforementioned LFC scheme, the IGBO-FOPID-based LFC scheme demonstrated significant advantages in achieving the minimum values of US/OS/ISE/ITSE/IAE/ITAE. Finally, the robustness of the proposed LFC scheme is analyzed by evaluating its performance with a 25% and 50% decrease in the inertia constant of the IPS. The evaluation results demonstrate that the proposed IGBO-FOPID-based LFC scheme maintains robust performance under nominal conditions, as well as with varying IPS inertia constants. In the future, this study may consider the latency of the system and adopt updated optimization techniques to design LFC schemes suitable for more complex modern interconnected power systems.