Parameter Identification of a Governing System in a Pumped Storage Unit Based on an Improved Artificial Hummingbird Algorithm

Abstract

Parameter identification is an important method to establish the governing system of a pumped storage unit. Appropriate parameters will make the governing system obtain better control performance. Therefore, in this study, an improved artificial hummingbird algorithm (IAHA) is proposed for the parameter identification of the governing system in a pumped storage unit. The algorithm integrates two key strategies to improve the optimization ability of the algorithm. First, the Chebyshev chaotic map is employed to initialize the artificial hummingbirds, which in turn increases and enhances the global search capability of the initial population. Second, the Levy flight is introduced in the guided foraging phase to expand the search space and avoid premature convergence. The performance of the proposed IAHA algorithm is compared with that of four other algorithms on 23 standard test functions, and the results show that IAHA has higher accuracy and faster convergence than the other four algorithms. Finally, IAHA was applied to the parameter identification of the governing system of a pumped storage unit to verify the effectiveness of the algorithm in tracking real-world problems.

1. Introduction

A pumped storage power station is an important energy storage technology. It can effectively alleviate the impact of intermittent fluctuation energy, such as scenery on the power system, and improve the absorption capacity of clean energy. To solve the security and stability of clean energy entering the power grid on a large scale [1]. Currently, the pumped storage power stations built or being built in China are developing in an intelligent direction. The complex water diversion system of a pumped storage power station brings great challenges to unit frequency governing and power grid security and stability. Therefore, it is necessary to carry out parameter identification of the regulation system of pumped storage units to ensure the safe and stable operation of units.

In modeling studies of governing systems, due to the diversity of unit and governor characteristics and their operating conditions, it is often difficult to derive accurate model parameters directly from their basic operating principles [2,3]; it is difficult to establish a complete system simulation model that can be used practically in power system simulation software or in the performance evaluation of governing systems. In order to solve this problem, the system identification method was widely used [4].

System identification is the determination of a mathematical model that describes the behavior of a system based on its input and output time functions. It consists of two basic parts: structure identification and parameter identification. The governing system is a complex dynamic system; the mathematical model structure can be determined by applying the mechanism analysis method, so the parameters of the model are determined by the parameter identification method [5]. In recent decades, the traditional methods for parameter identification were the least squares method [6], input response method [7], and maximum likelihood estimation method [8], but these methods have many limitations, for example, the least squares method requires sufficient system inputs, and the maximum likelihood estimation method can easily fall into local optimality, etc. In recent years, identification methods based on meta-heuristic algorithms have been developed that treat the parameter identification problem as an optimization problem. Since metaheuristic algorithms are a global optimization method, they can establish an objective function for parameter identification by optimizing the output error between the original system and the identification system, comparing with traditional identification methods, metaheuristic algorithms are more suitable for the parameter identification of complex systems, and the identification performances of these algorithms depend on optimization capability.

Many metaheuristic algorithms have been presented and successfully applied to different areas: particle swarm optimization algorithm ((PSO) simulates the foraging behavior of birds in a group [9,10]; the chimpanzee optimization algorithm (ChOA) simulates the cooperative hunting behaviors of attack, drive, intercept, and chase chimpanzees [11]; the artificial rabbit optimization (ARO) algorithm is inspired by the survival strategies of rabbits found in nature, including detour foraging and random hiding [12]; the genetic algorithm (GA) simulates biological evolutionary mechanisms in natural environments [13]; the artificial ecosystem optimization algorithm (AEO) simulates the energy flow process in the earth ecosystem [14]; the gravitational search algorithm (GSA) involves the optimization of populations based on the law of gravity and Newton’s second law [15]; the ant colony optimization (ACO) algorithm is a simulation of the way ants find paths in nature [16]; the Black Widow algorithm (BWO) simulates its entire life cycle [17]; the Tom search algorithm (ASO) simulates the displacement of atoms in a molecular system composed of atoms due to their mutual force and system constraint [18]; the gray wolf optimization (GWA) algorithm simulates gray wolf prey predation activities [19]; the artificial fish swarming algorithm (AFSA) simulates the foraging, clustering, and tail-chasing behaviors of fish [20]. The algorithm for the ant lion optimization algorithm (ALO) simulates the hunting mechanism of an ant lion hunting ants [21]; the whale optimization algorithm (WOA) is based on encircling prey, bubble netting prey, and searching for whale prey [22]; the manta ray foraging optimization algorithm (MRFO) simulates the foraging process of manta rays in the ocean [23]. Strong robustness, adaptability, and randomness are characteristics of these optimization algorithms. Although these metaheuristics outperform conventional numerical approaches in handling challenging engineering problems, they also have some drawbacks and nevertheless hold enormous promise for performance in optimization.

Some authors have made improvements to metaheuristic algorithms, for example, Zhongqiang Wu et al. [24] proposed an improved ant lion optimization algorithm to identify the parameters of the solar cell model by adding chaotic sequences. M Ali et al. proposed an algorithm to identify the parameters of a polymer electrolyte membrane fuel cell model using the gray wolf optimization algorithm [25]. Xiao Zhang et al. applied an elite backward learning particle swarm algorithm for the identification of PV cell parameters [26]. The metaheuristic algorithms have some limitations due to local optimum and premature phenomena, despite the fact that these algorithms have been successfully used in a variety of fields for parameter identification problems.

In 2022, Weiguo Zhao et al. proposed the artificial hummingbird algorithm (AHA) [27], which was inspired by simulating hummingbirds, special flight abilities, and their intelligent foraging strategies. The method’s advantages—it has few parameters, a fast speed, and performs well in solving optimization problems. This paper proposes an improved metaheuristic algorithm, named the improved artificial hummingbird algorithm (IAHA). In order to make the initialization more uniform and rich, IAHA added the Chebyshev chaotic map to initialize the artificial hummingbird and the Levy flight to improve the search efficiency when guiding foraging.

2. Artificial Hummingbird Algorithm (AHA)

2.1. Brief Introduction of AHA

AHA is a population-based metaheuristic algorithm that mainly simulates three foraging behaviors of hummingbirds: guided foraging, territorial foraging, and migratory foraging. During the foraging process, three flight skills are modeled—omnidirectional, diagonal, and axial flight. At the same time, an access table simulating the hummingbird’s extraordinary memory ability is constructed to guide the hummingbird to perform global optimization in the algorithm.

The three flying skills are defined as follows:

The flight skill simulation is extended to the d-D space with axial flight defined as follows:

$$D^{\left(i\right)}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}i=randi\left(\left[1,d\right]\right)i=1,\cdots ,d\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(1)

Diagonal flight is defined as follows:

$$D^{\left(i\right)}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}i=p\left(j\right)P=randperm\left(k\right),k\in [2,\lceil {\mathrm{r}}_1(d-2)\rceil +1]\hfill \\ 0,\hfill & \mathrm{else}.\hfill \end{array}\right.$$

(2)

Omnidirectional flight is defined as follows:

$$D^{\left(i\right)}=1i=1,\cdots ,d$$

(3)

where $randi\left(\right[1,d\left]\right)$ generates a random integer from 1 to d, randperm(k) creates a random permutation of integers from 1 to k, and $r_1$ is a random number in (0, 1]. The diagonal flight in a d-D space is inside a hyperrectangle.

The AHA first initialization a set of random solutions and a visit table. In each iteration, guided or territorial foraging is performed 50% of the time. Hummingbirds can migrate toward their intended food sources using guided foraging, which is based on nectar filling rates and a visit table. Territorial foraging allows hummingbirds to easily move to neighboring regions within their own territory and find new food sources as candidates. Migration foraging is performed every two iterations. Up until the stop rule is reached, all operations and calculations are performed interactively. Finally, the food source with the highest rate of nectar-refilling replenishment is returned as an approximate global optimum.

(1) A population of n hummingbirds are randomly initialized to n food sources as follows:

$$x_i=Low+r\times \left(Up-Low\right)i=1,\cdots \cdots n$$

(4)

where $Low$ and $Up$ are the lower and upper boundaries for a d-dimensional problem, respectively, r is a random vector in $[0,1]$, and $x_i$ represents the position of the ith food source.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill V{T}_i{,}_j=\left\{\begin{array}{cc}0\hfill & ifi\ne j\hfill \\ null\hfill & i=j\hfill \end{array}\right.\end{array}\end{array}$$

(5)

where $i=j,V{T}_i{,}_j=null$ indicates that a hummingbird is taking food at its specific food source; $i\ne j,V{T}_i{,}_j=0$ indicates that the jth food source has just been visited by the ith hummingbird in the current iteration.

(2) Guided foraging: With the aforementioned flight capabilities, the hummingbird can access its target food source to obtain candidate food sources, so the mathematical equation for simulating guiding foraging behavior and candidate food sources is as follows:

$$v_i\left(t+1\right)=x_{i,tar}\left(t\right)+a\times D\times (x_i\left(t\right)-x_{i,tar}\left(t\right))$$

(6)

$$a\sim N\left(0,1\right)$$

(7)

where $x_i\left(t\right)$ is the position of the ith hummingbird food source in time t, $x_{i,tar}\left(t\right)$ is the position of the ith hummingbird target food source in time t, a is normally distributed, with mean = 0 and a standard deviation of 1.

The position update of the ith food source is as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill x_i\left(t+1\right)=\left\{\begin{array}{cc}{x}_i\left(t\right)\hfill & f\left(x_i\left(t\right)\right)\le f{(v}_i\left(t+1\right))\hfill \\ v_i\left(t+1\right)\hfill & f\left(x_i\left(t\right)\right)>f{(v}_i\left(t+1\right))\hfill \end{array}\right.\end{array}\end{array}$$

(8)

where f (·) indicates the function fitness value. Equation (8) shows that if the nectar refilling rate of the candidate food source is better than that of the current one, the hummingbird abandons the current food source and stays at the candidate one resulting from Equation (6) for feeding.

(3) Territorial foraging: After reaching a target food source where nectar is eaten, hummingbirds may seek new food sources. Therefore, a hummingbird can easily move to a neighboring region within its territory, where new food sources can be found that may be better candidate solutions. The mathematical equation for simulating the local search of hummingbirds in territorial foraging strategies and candidate food sources is as follows:

$$v_i\left(t+1\right)=x_i\left(t\right)+b\times D\times x_i\left(t\right)$$

(9)

$$b\sim N\left(0,1\right)$$

(10)

where b is normally distributed, with mean = 0 and a standard deviation of 1.

(4) When food is often scarce in a territory frequented by hummingbirds, the bird often migrates to more distant food sources to forage. In the AHA algorithm, a migration coefficient is defined. The hummingbirds in the food source with the weakest filling rate will, at random, migrate to another new food source throughout the entire search space when the number of iterations exceeds the predefined value of the migration coefficient. At this point, the hummingbird will abandon the original source and stay at the new source foraging, and then the migration foraging of the hummingbird from the source with the worst nectar-filling rate to the randomly generated new source can be given as follows:

$$x_{wor}\left(t+1\right)=Low+r\times \left(Up-Low\right)$$

(11)

where $x_{wor}$ is the food source with the worst nectar refilling rate in the population.

2.2. Improved Artificial Hummingbird Algorithm Based on Chebyshev Chaotic Map and Levy Flight (IAHA)

2.2.1. Chebyshev Chaotic Map

Chaos is a random state in a deterministic system, a phenomenon of the evolution of a nonlinear system. It is caused by deterministic rules and is very sensitive to the initial conditions, which is long-term behavior without a fixed period. Chaos has an ergodic, and searching using chaotic variables is obviously superior to an unordered search [28].

The Chebyshev chaotic map has a wide and more uniform distribution range, and it can be distributed in the interval [−1, 1]. When $k\ge 2$ (k is the order), no matter how close the initial value is chosen, the iterated sequence is uncorrelated and chaotic and ergodic within this range. The equation is shown as follows:

$$x_{n+1}=cos\left(karccos{x}_n\right)x_n\in [-1,1]$$

(12)

This paper uses the equation to generate uniformly distributed points to initialize the position of artificial hummingbirds, improve the global search ability of the initial population, and improve the solution accuracy of the algorithm.

2.2.2. Levy Flight

In the 1930s, the Levy distribution was a probability distribution proposed by Levy, a French mathematician. After that, many researchers have carried out many studies on the Levy distribution. So far, it has been proven that the foraging trajectory of many animals in nature follows the Levy distribution. The Levy flight follows the principle of the Levy distribution of many random phenomena, such as Brownian motion, random walk, etc. [29] At present, Levy flight is widely used in intelligent optimization. For example, the Cuckoo algorithm adopts Levy flight to update the position [30]. Levy flight can expand the search space, so it is easier to avoid premature convergence by introducing Levy flight into the AHA algorithm.

Levy flight position updated to:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill x_i\left(t+1\right)=\left\{\begin{array}{cc}{x}_i\left(t\right)+\alpha \oplus Levy\left(\lambda \right)\hfill & f\left(x_i\left(t\right)\right)\le f{(v}_i\left(t+1\right))\hfill \\ v_i\left(t+1\right)\hfill & f\left(x_i\left(t\right)\right)>f{(v}_i\left(t+1\right))\hfill \end{array}\right.\end{array}\end{array}$$

(13)

where $x_i^t$ is the tth generation position of $x_i$, ⊕ is the dot multiplication, $\alpha$ is the step size control parameter, and Levy$\left(\lambda \right)$ is the random search path, which satisfies:

$$Levy\sim u=t^{-\lambda },1<\lambda \le 3$$

(14)

Its step size obeys the Levy distribution, and step size s is calculated as [29]:

$$s={\displaystyle \frac{\mu }{{\left|v\right|}^{1/\beta }}}$$

(15)

where $\mu ,v$ are normally distributed, defined as:

$$\mu \sim N(0,{\sigma }_{\mu }^2)$$

$$v\sim N(0,{\sigma }_v^2)$$

where

$${\sigma }_{\mu }={\displaystyle \frac{(1+\beta )(sin{\displaystyle \frac{\pi \beta }{2}})}{{\displaystyle \frac{1+\beta }{2}}{\beta }^2{\displaystyle \frac{\beta -1}{2}}}}$$

(16)

$${\sigma }_v=1$$

where $\beta$ is usually a constant of 1.5.

According to the above, the flow chart of the IAHA algorithm is shown in Figure 1. Its basic steps are as follows:

Step 1: Set the parameters: a maximum number of iterations of T, a number of artificial hummingbirds of Pop, and a dimension of the fitness function of $Dim$.

Step 2: The Chebyshev chaotic map was used to initialize the food source location of artificial hummingbirds; the corresponding fitness function was calculated and the optimal value was recorded.

Step 3: Introducing Levy flight-guided foraging or territorial foraging to update the food source location of artificial hummingbirds. The probability of both forages was 50%.

Step 4: In the worst-case scenario of a food source, hummingbirds may visit the same food source as their target source after two repetitions, in which case, migratory foraging is performed.

Step 5: Determine if the result has reached the maximum number of iterations of the algorithm, if it has reached the maximum accuracy, then the optimal food source location is obtained, otherwise, transpose step 3.

3. Performance Testing and Analysis

In order to verify the effectiveness of the IAHA algorithm. We compare the IAHA method with four well-established optimizers, including the artificial hummingbird algorithm (AHA), particle swarm optimization algorithm (PSO), ant lion optimization algorithm (ALO), and gravitational search algorithm (GSA). In addition, the number of iterations of all test functions in the five algorithms was 500, and the population size was 30. The rest of the parameters are shown in

Table 1. Parameter settings for each algorithm in the test function.

Algorithm	Parameter Settings
PSO	$c_1=2;c_2=2;{\omega }_{max}=0.8;{\omega }_{min}=0.2$$\omega$ decreases linearly from 0.8 to 0.2
GSA	$G_0=100;a=20$

.

The metaheuristic algorithm is a kind of random search algorithm. For the same optimization problem, the results are usually not identical if the same algorithm is used many times. To avoid excessive errors arising from randomness in a single result, the algorithm was repeated 20 times and the statistical results for each algorithm are shown in

Table 2. Statistical results of unimodal test functions for each algorithm.

Function	Index	IAHA	AHA	PSO	ALO	GSA
F1	Mean	$9.48\times {10}^{-141}$	$1.72\times {10}^{-136}$	$3.05\times {10}^{-2}$	$1.34\times {10}^{-3}$	$2.27\times {10}^{-16}$
	Std	$4.08\times {10}^{-140}$	$7.67\times {10}^{-136}$	$7.57\times {10}^{-2}$	$7.96\times {10}^{-4}$	$9.07\times {10}^{-17}$
	Best	$8.48\times {10}^{-163}$	$1.19\times {10}^{-158}$	$3.50\times {10}^{-4}$	$2.38\times {10}^{-4}$	$7.15\times {10}^{-17}$
	Worst	$1.83\times {10}^{-139}$	$3.43\times {10}^{-135}$	$3.31\times {10}^{-1}$	$3.06\times {10}^{-3}$	$4.85\times {10}^{-16}$
F2	Mean	$9.65\times {10}^{-74}$	$1.85\times {10}^{-72}$	$4.28\times {10}^{-2}$	$38.00$	$1.83\times {10}^{-4}$
	Std	$2.96\times {10}^{-73}$	$8.18\times {10}^{-72}$	$7.35\times {10}^{-2}$	$40.06$	$4.97\times {10}^{-4}$
	Best	$1.36\times {10}^{-83}$	$4.7\times {10}^{-82}$	$3.44\times {10}^{-3}$	$4.77$	$3.66\times {10}^{-8}$
	Worst	$1.07\times {10}^{-72}$	$1.16\times {10}^{-71}$	$2.59\times {10}^{-1}$	1.19 $\times {10}^2$	$1.70\times {10}^{-3}$
F3	Mean	$1.05\times {10}^{-129}$	$6.04\times {10}^{-128}$	$4.10\times {10}^4$	$4.80\times {10}^4$	$9.90\times {10}^3$
	Std	$4.70\times {10}^{-129}$	$2.19\times {10}^{-127}$	$1.36\times {10}^3$	$2.72\times {10}^3$	$3.37\times {10}^2$
	Best	$2.99\times {10}^{-145}$	$2.26\times {10}^{-148}$	$1.67\times {10}^3$	$8.02\times {10}^2$	$4.31\times {10}^2$
	Worst	$2.10\times {10}^{-128}$	$9.63\times {10}^{-127}$	$7.66\times {10}^3$	$1.05\times {10}^4$	$1.78\times {10}^3$
F4	Mean	$7.15\times {10}^{-65}$	$4.80\times {10}^{-63}$	$19.51$	$17.64$	$7.13$
	Std	$1.95\times {10}^{-64}$	$1.94\times {10}^{-62}$	$4.18$	$4.58$	$1.91$
	Best	$1.31\times {10}^{-73}$	$1.45\times {10}^{-70}$	$12.50$	$11.90$	$3.35$
	Worst	$7.98\times {10}^{-64}$	$8.67\times {10}^{-62}$	$31.18$	$24.64$	$10.76$
F5	Mean	$26.64$	$26.80$	$1.80\times {10}^2$	$2.93\times {10}^2$	$62.28$
	Std	$3.17\times {10}^{-1}$	$3.69\times {10}^{-1}$	$1.26\times {10}^2$	$3.85\times {10}^2$	$56.67$
	Best	$26.04$	$26.10$	$52.86$	$27.15$	$27.09$
	Worst	$27.29$	$27.46$	$4.91\times {10}^2$	$1.22\times {10}^3$	$2.25\times {10}^2$
F6	Mean	0	0	$1.8$	$1.84\times {10}^{-3}$	$6.7$
	Std	0	0	$1.77$	$1.55\times {10}^{-3}$	$8.77$
	Best	0	0	0	$3.80\times {10}^{-4}$	0
	Worst	0	0	6	$5.27\times {10}^{-3}$	30
F7	Mean	$1.90\times {10}^{-4}$	$2.48\times {10}^{-4}$	$9.16\times {10}^{-2}$	$2.53\times {10}^{-1}$	$8.68\times {10}^{-2}$
	Std	$1.61\times {10}^{-4}$	$1.84\times {10}^{-4}$	$3.21\times {10}^{-2}$	$6.79\times {10}^{-2}$	$3.82\times {10}^{-2}$
	Best	$1.27\times {10}^{-5}$	$3.79\times {10}^{-5}$	$2.39\times {10}^{-2}$	$1.15\times {10}^{-1}$	$2.95\times {10}^{-2}$
	Worst	$5.23\times {10}^{-4}$	$6.66\times {10}^{-4}$	$1.56\times {10}^{-1}$	$3.51\times {10}^{-1}$	$1.65\times {10}^{-1}$

,

Table 3. Statistical results of high-dimensional multi-peak test functions for each algorithm.

Function	Index	IAHA	AHA	PSO	ALO	GSA
F8	Mean	$-1.15\times {10}^4$	$-1.16\times {10}^4$	$-5.78\times {10}^3$	$-5.66\times {10}^3$	$-2.61\times {10}^3$
	Std	$3.73\times {10}^2$	$3.50\times {10}^2$	$9.04\times {10}^2$	$8.60\times {10}^2$	$4.75\times {10}^2$
	Best	$-1.23\times {10}^4$	$-1.22\times {10}^4$	$-7.35\times {10}^3$	$-9.29\times {10}^3$	$-4.06\times {10}^3$
	Worst	$-1.07\times {10}^4$	$-1.08\times {10}^4$	$-4.30\times {10}^3$	$-5.42\times {10}^3$	$-1.87\times {10}^3$
F9	Mean	0	0	$41.99$	$75.67$	$31.44$
	Std	0	0	$11.33$	$21.79$	$6.73$
	Best	0	0	$25.87$	$30.84$	$18.90$
	Worst	0	0	$72.75$	$1.16\times {10}^2$	$43.78$
F10	Mean	$8.88\times {10}^{-16}$	$8.88\times {10}^{-16}$	$1.4$1	$5.68$	$1.22\times {10}^{-8}$
	Std	0	0	$8.32\times {10}^{-1}$	$3.48$	$3.71\times {10}^{-9}$
	Best	$8.88\times {10}^{-16}$	$8.88\times {10}^{-16}$	$1.24\times {10}^{-2}$	$2.32$	$4.95\times {10}^{-9}$
	Worst	$8.88\times {10}^{-16}$	$8.88\times {10}^{-16}$	$2.50$	$12.67$	$1.88\times {10}^{-8}$
F11	Mean	0	0	$6.77\times {10}^{-2}$	$6.83\times {10}^{-2}$	$25.62$
	Std	0	0	$1.23\times {10}^{-1}$	$4.03\times {10}^{-2}$	$5.87$
	Best	0	0	$2.37\times {10}^{-3}$	$1.85\times {10}^{-2}$	$15.83$
	Worst	0	0	$5.73\times {10}^{-1}$	$1.87\times {10}^{-1}$	$36.56$
F12	Mean	$6.78\times {10}^{-4}$	$8.67\times {10}^{-4}$	$2.21$	$12.36$	$1.97$
	Std	$1.54\times {10}^{-3}$	${1.99}^{-3}$	$1.42$	$4.69$	$1.07$
	Best	$3.71\times {10}^{-5}$	$7.02\times {10}^{-5}$	$1.99\times {10}^{-1}$	$5.59$	$1.91\times {10}^{-1}$
	Worst	$6.67\times {10}^{-3}$	$7.34\times {10}^{-3}$	$5.63$	$23.23$	$4.21$
F13	Mean	$1.82$	$1.96$	$3.19$	$30.74$	$7.37$
	Std	$5.16\times {10}^{-1}$	$3.39\times {10}^{-1}$	$4.33$	$17.89$	$5.53$
	Best	$5.25\times {10}^{-1}$	$1.40$	$2.74\times {10}^{-1}$	$2.60\times {10}^{-2}$	$1.08\times {10}^{-1}$
	Worst	$2.48$	$2.53$	$16.02$	$62.73$	$22.86$

and

Table 4. Statistical results of low-dimensional multi-peak test functions of each algorithm.

Function	Index	IAHA	AHA	PSO	ALO	GSA
F14	Mean	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$3.017$	$6.38$
	Std	0	0	$5.09\times {10}^{-17}$	$2.76$	$3.84$
	Best	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$1.06$
	Worst	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$9.98\times {10}^{-1}$	$10.76$	$13.32$
F15	Mean	$3.07\times {10}^{-4}$	$3.07\times {10}^{-4}$	$4.97\times {10}^{-4}$	$2.27\times {10}^{-3}$	$6.03\times {10}^{-3}$
	Std	$4.21\times {10}^{-11}$	$1.54\times {10}^{-9}$	$1.32\times {10}^{-4}$	$4.47\times {10}^{-3}$	$4.21\times {10}^{-3}$
	Best	$3.07\times {10}^{-4}$	$3.07\times {10}^{-4}$	$3.07\times {10}^{-4}$	$6.14\times {10}^{-4}$	$8.10\times {10}^{-4}$
	Worst	$3.07\times {10}^{-4}$	$3.07\times {10}^{-4}$	$6.84\times {10}^{-4}$	$2.06\times {10}^{-2}$	$1.68\times {10}^{-2}$
F16	Mean	$-1.03$	$-1.03$	$-1.03$	$-1.03$	$-1.03$
	Std	$1.76\times {10}^{-16}$	$1.69\times {10}^{-16}$	$2.22\times {10}^{-16}$	$1.04\times {10}^{-13}$	$1.02\times {10}^{-16}$
	Best	$-1.03$	$-1.03$	$-1.03$	$-1.03$	$-1.03$
	Worst	$-1.03$	$-1.03$	$-1.03$	$-1.03$	$-1.03$
F17	Mean	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$
	Std	0	0	$6.06\times {10}^{-10}$	$5.26\times {10}^{-14}$	0
	Best	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$
	Worst	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$	$3.98\times {10}^{-1}$
F18	Mean	3	3	3	3	3
	Std	$2.04\times {10}^{-16}$	$3.95\times {10}^{-16}$	$1.28\times {10}^{-15}$	$6.49\times {10}^{-13}$	$4.34\times {10}^{-15}$
	Best	3	3	3	3	3
	Worst	3	3	3	3	3
F19	Mean	$-3.86$	$-3.86$	$-3.86$	$-3.86$	$-3.86$
	Std	$2.24\times {10}^{-15}$	$2.26\times {10}^{-15}$	$2.22\times {10}^{-15}$	$1.06\times {10}^{-12}$	$1.87\times {10}^{-15}$
	Best	$-3.86$	$-3.86$	$-3.86$	$-3.86$	$-3.86$
	Worst	$-3.86$	$-3.86$	$-3.86$	$-3.86$	$-3.86$
F20	Mean	$-3.32$	$-3.30$	$-3.26$	$-3.27$	$-3.32$
	Std	$4.20\times {10}^{-12}$	$4.36\times {10}^{-2}$	$6.10\times {10}^{-2}$	$6.02\times {10}^{-2}$	$2.67\times {10}^{-2}$
	Best	$-3.32$	$-3.32$	$-3.32$	$-3.32$	$-3.32$
	Worst	$-3.32$	$-3.30$	$-3.30$	$-3.30$	$-3.30$
F21	Mean	$-10.15$	$-10.15$	$-7.50$	$-6.49$	$-6.37$
	Std	$2.31\times {10}^{-6}$	$6.22\times {10}^{-6}$	$3.19$	$2.88$	$3.66$
	Best	$-10.15$	$-10.15$	$-10.15$	$-10.15$	$-10.15$
	Worst	$-10.15$	$-10.15$	$-2.68$	$-2.63$	$-2.68$
F22	Mean	$-10.40$	$-10.40$	$-7.54$	$-6.65$	$-9.86$
	Std	$1.19\times {10}^{-9}$	$1.55\times {10}^{-5}$	$3.58$	$3.57$	$1.82$
	Best	$-10.40$	$-10.40$	$-10.40$	$-10.40$	$-10.40$
	Worst	$-10.40$	$-10.40$	$-2.75$	$-2.75$	$-2.75$
F23	Mean	$-10.54$	$-10.54$	$-9.05$	$-7.07$	$-10.13$
	Std	$6.20\times {10}^{-9}$	$9.08\times {10}^{-4}$	$3.06$	$3.62$	$1.81$
	Best	$-10.54$	$-10.54$	$-10.54$	$-10.54$	$-10.54$
	Worst	$-10.54$	$-10.54$	$-2.81$	$-2.42$	$-2.43$

. The tables include the mean (mean), standard deviation (Std), best (best), and worst (worst) values of the 23 test functions (see

Table A1. Standard test function.

Function	D	Range
$f_1\left(x\right)={\sum }_{i=1}^n{x}_i^2$	30	${[-100,100]}^D$
$f_2\left(x\right)={\sum }_{i=1}^n|x_i|+{\prod }_{i=1}^n\left|x_i\right|$	30	${[-100,100]}^D$
$f_3\left(x\right)={\sum }_{i=1}^n({\sum }_{j=1}^i{x_j)}^2$	30	${[-100,100]}^D$
$f_4\left(x\right)={max}_i\left\{|x_i|,1\le i\le n\right\}$	30	${[-100,100]}^D$
$f_5\left(x\right)={\sum }_{i=1}^{n-1}(100{(x_{i+1}-x_i)}^2)+\left({x_i-1)}^2\right)$	30	${[-30,30]}^D$
$f_6\left(x\right)={\sum }_{i=1}^n{{(x}_i+0.5)}^2$	30	${[-100,100]}^D$
$f_7\left(x\right)={\sum }_{i=1}^{n}i{x}_i^4+random[0,1)$	30	${[-1.28,1.28]}^D$
$f_8\left(x\right)=-{\sum }_{i=1}^n(x_{i}sin(\sqrt{|x_i|})$	30	${[-500,500]}^D$
$f_9\left(x\right)={\sum }_{i=1}^n{(x_i^2-10cos\left(2\pi x_i\right)+10)}^2$	30	${[-5.12,5.12]}^D$
$f_{10}\left(x\right)=-20exp\left(-0.2\sqrt{\frac{1}{n}{\sum }_{i=1}^n{x}_i^2}\right)-exp(\frac{1}{n}{\sum }_{i=1}^{n}cos2\pi x_i)+20+e$	30	${[-32,32]}^D$
$f_{11}\left(x\right)=\frac{1}{4000}{\sum }_{i=1}^n{(x_i-100)}^2-{\prod }_{i=1}^{n}cos\left(\frac{x_i-100}{\sqrt{i}}\right)+1$	30	${[-600,600]}^D$
$f_{12}\left(x\right)=\frac{\pi }{n}\left\{10{sin}^2\left(\pi y_1\right)+{\sum }_{i=1}^{n-1}{{(y}_1-1)}^2\left[1+10{sin}^2\left(\pi y_1+1\right)\right]+({y_n-1)}^2\right\}+{\sum }_{i=1}^{30}u(x_i,10,100,4)$	30	${[-50,50]}^D$
$f_{13}\left(x\right)=0.1\left\{si{n}^2\left(3\pi x_1\right)+{\sum }_{i=1}^{29}{(x_i-1)}^{2}p[1+10si{n}^2\left(3\pi x_{i+1}\right)]+{(x_n-1)}^2[1+si{n}^2\left(2\pi x_{30}\right)]\right\}+{\sum }_{i=1}^{30}u(x_i,5,10,4)$	30	${[-50,50]}^D$
$f_{14}\left(x\right)={[{\displaystyle \frac{1}{500}}+{\sum }_{j=1}^{25}{\displaystyle \frac{1}{j+{\sum }_{j=1}^2{(x_i-a_i)}^6}}]}^{-1}$	2	${[-65.536,65.536]}^D$
$f_{15}\left(x\right)={{\sum }_{i=1}^{11}|a_i-\frac{x_1(b_i^2+b_i{x}_2)}{b_i^2+b_i{x}_3+x_4}|}^2$	4	${[-5,5]}^D$
$f_{16}\left(x\right)=4x_1^2+2.1x_1^4+\frac{1}{3}{x}_1^6+x_1{x}_2-4x_2^2+4x_2^4$	2	${[-5,5]}^D$
$f_{17}\left(x\right)={(x_2-\frac{5.1}{4{\pi }^2}{x}_1^2+\frac{5}{\pi }{x}_1-6)}^2+10(1-\frac{1}{8\pi }cos{x}_1)+10$	2	$\left[-5,10\right]\times [0,15]$
$f_{18}\left(x\right)=\left[1+{\left(x_1+x_2+1\right)}^2\left(19-14x_1+3x_1^2-14x_2+6x_1{x}_2+3x_2^2\right)\right]\times [30+{\left(2x_1+1-3x_2\right)}^2\left(18-32x_1+12x_1^2+48x_2-36x_1{x}_2+27x_2^2\right)]$	2	${[-2,2]}^D$
$f_{19}\left(x\right)=-{\sum }_{i=1}^{4}exp[-{\sum }_{j=1}^3{a}_{ij}{(x_j-p_{ij})}^2]$	3	${[0,1]}^D$
$f_{20}\left(x\right)=-{\sum }_{i=1}^{4}exp[-{\sum }_{j=1}^6{a}_{ij}{(x_j-p_{ij})}^2]$	6	${[0,1]}^D$
$f_{21}\left(x\right)=-{\sum }_{i=1}^5|{\left(x_i-a_i\right)\left(x_i-a_i\right)}^T+c_i{|}^{-1}$	4	${[0,10]}^D$
$f_{22}\left(x\right)=-{\sum }_{i=1}^7|{\left(x_i-a_i\right)\left(x_i-a_i\right)}^T+c_i{|}^{-1}$	4	${[0,10]}^D$
$f_{23}\left(x\right)=-{\sum }_{i=1}^{10}|{\left(x_i-a_i\right)\left(x_i-a_i\right)}^T+c_i{|}^{-1}$	4	${[0,10]}^D$

in “Appendix A”), in which the bold are the best values of the mean, standard deviation, best, and worst values of the different test functions for the five algorithms.

3.1. Unimodal Test Functions

F1–F7 are unimodal functions, which mainly test the capability of the development stage. Moreover, these types of functions do not have local extrema, but only global optimal values, which are easy to optimize. Based on this, the convergence speed may be more important than the global optimal value of the algorithm. The statistical results of the IAHA algorithm proposed compared with the other four algorithms are shown in Table 2 and Figure 2. Among the five intelligent optimization algorithms, the IAHA algorithm has a significant advantage in terms of mean and standard deviation, except for F6, where both the IAHA algorithm and the AHA algorithm achieve optimal values. As a result, the IAHA algorithm is superior in development ability, convergence, and stability.

3.2. High-Dimensional Multi-Peak Test Functions

F8–F23 are multi-peak functions, which mainly demonstrate the ability of the exploration phase. Among these, F8–F13 are high-dimensional multi-peak functions and F14–F23 are low-dimensional multi-peak functions. The multi-peak test functions have multiple local extremes and global optima, which are more difficult to optimize compared to unimodal functions. As a result, the multi-peak function algorithm’s capability of global optimization is given additional focus. The statistical results of the high-dimensional multi-peak function test for each algorithm are shown in Table 3 and Figure 3. In F8, the mean value and standard deviation of IAHA rank second, second only to AHA. In F9–F11, the mean and standard deviation of IAHA and AHA are the same. The mean value of IAHA (in F12 and F13) is the best. Although the standard deviation of F13 is not the smallest, its best value and the worst value are optimal compared with other algorithms. Thus, it can be seen that IAHA has a better advantage in the global search.

3.3. Low-Dimensional Multi-Peak Test Functions

F14–F23 are low-dimensional multi-peaked functions, which have fewer local extremes than higher-dimensional multi-peaked functions and are relatively easier to optimize. Table 4 and Figure 4 show the statistical results of the low-dimensional multi-peak test functions of each algorithm. The mean values of the functions in F16–F19 tested by the IAHA algorithm and other algorithms are the same and reach optimal values. In addition, the standard deviation of the IAHA algorithm in F16 and F19 is second only to the GSA algorithm. In F20, the result of the IAHA algorithm is slightly better than that of other algorithms. The IAHA algorithm has a minimum standard deviation in F15, F18, F19, F21–F23. However, regarding the mean values of the AHA algorithm and IAHA algorithm, eight functions are the same and do not show better advantages. Therefore, because of the particularity of the function, the results obtained by the proposed IAHA algorithm and other algorithms are basically close to the theoretical global optimum.

4. Analysis of Parameter Identification under Frequency Disturbance Conditions

4.1. Mathematical Model of Governor

4.1.1. PID Controller Model

PID controllers (due to their simple structures, high stability, and convenient adjustments) are the main technologies of industrial control. At present, PID controllers are used in most hydropower stations in China, and their transfer functions are as follows [31]:

$$G_{PID}\left(s\right)={\displaystyle \frac{K_p+\frac{K_i}{s}+\frac{K_d}{(T_{1v}s+1)}}{\frac{b_p{K}_i}{s}+1}}$$

(17)

where $K_p$ is the proportional parameter, $K_i$ is the integral parameter, $K_d$ is the derivative parameter, $T_{1v}$ is the derivative time constant, $b_p$ is the permanent speed droop, and s is the Laplace operator.

4.1.2. Electrohydraulic Servo System Model

The electrohydraulic servo system is the executive agency of the governor, which consists of the comprehensive amplification element, electrohydraulic converter, the main control valve, and the main servomotor [32]. Considering that the time constant of the main servomotor is much larger than that of the auxiliary servomotor, the time constant of the auxiliary servomotor is ignored. In order to prevent the guide vane from moving too fast, the main servomotor usually has a limited amplitude (in regard to achieving the speed limit of the main servomotor) [33]:

$$G_y\left(s\right)=\frac{1}{T_{y}s+1}$$

(18)

where $T_y$ is the servomotor time constant.

4.2. Mathematical Model of the Pressure Water Diversion System

The pressure diversion system is an important hydraulic construction of the pumped storage power station; it mainly consists of the upper reservoir, headrace tunnel, upper surge chamber, penstock, casing, turbine pump, draft tube, downstream surge chamber, tailwater tunnel, and lower reservoir [34]. The two types of water hammer models for the water diversion system are elastic and rigid. In this paper, the elastic water hammer model was chosen, and the hyperbolic tangent function was used as the mathematical model as follows:

$$G_h\left(s\right)=\frac{H\left(s\right)}{Q\left(s\right)}=-2h_w\frac{e^{\frac{T_r}{2}s}-e^{-\frac{T_r}{2}s}}{e^{\frac{T_r}{2}s}+e^{-\frac{T_r}{2}s}}$$

$$=-2h_w\frac{sh\left(\frac{T_r}{2}s\right)}{ch\left(\frac{T_r}{2}s\right)}$$

(19)

where $T_r$ is the water hammer pressure wave time constant, $h_w$ is the pipeline characteristic coefficient, Q is the flow rate, and H is the head.

The rigid water hammer model is unable to accurately represent the dynamic properties of water flow in the pipe because it ignores the frictional resistance of water flow as well as the elasticity of the water body and pipe wall. Therefore, this paper makes use of the fourth-order elastic water hammer model. The transfer function is as follows:

$$G_h\left(s\right)=\frac{H\left(s\right)}{Q\left(s\right)}=-h_w\frac{T_{r}s+\frac{1}{24}{T}_r^3{s}^3}{1+\frac{1}{8}{T}_r^2{s}^2+\frac{1}{384}{T}_r^4{s}^4}$$

(20)

4.3. Turbine Pump Mathematical Model

Under some operating conditions, the turbine pump operates with small fluctuations around the specified operating point, and the variation in the nonlinear characteristics of the turbine pump is not particularly significant within a small range. In this case, the full characteristic model of the turbine pump is simplified at a certain operating point to obtain a linear model of the turbine pump, which is calculated as follows [35]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left\{\begin{array}{c}{m}_t\left(t\right)=e_{x}x\left(t\right)+e_{y}y\left(t\right)+e_{h}h\left(t\right)\hfill \\ q\left(t\right)=e_{qx}x\left(t\right)+e_{qy}y\left(t\right)+e_{qh}h\left(t\right)\hfill \end{array}\right.\end{array}\end{array}$$

(21)

where x is speed, y is guide vane opening, h is working head, $e_x,e_y,e_h,e_{qx},e_{qy},e_{qh}$ are turbine transfer coefficient at a certain operating point.

4.4. Generator and Load Mathematical Model

The first-order model (reflecting the dynamic characteristics and self-regulating performance of the generator rotor) is usually used in research on the modeling and optimization control of the governing system of the pumped storage unit. Considering the balance between the torque of the turbine pump and the torque of the generator, the transfer function of the generator and load is [36]:

$$G_g\left(s\right)=\frac{1}{T_{a}s+e_g}$$

(22)

where $T_a$ is the inertial time constant of the generator and $e_g$ is the adjustment coefficient of the generator.

A six-parameter model was used for the turbine pump in the pumped storage unit, and the water diversion system adopted the four-order elastic water hammer model, where the model took into account the nonlinear link of the electrohydraulic follower system with limiting amplitude; the model structure of the governing system of the pumped storage unit is shown in Figure 5.

4.5. Parameter Identification of the Governing System

The parameter identification with a known model structure is usually transformed into an optimization problem [37]. The optimization objective function is the deviation between the system output variable and the identified system output variable, and then the integration is performed. Under the assumption that the identification system can reflect the real system, the optimal value can be found through unknown parameters to achieve the identification effect and make the output consistent with the real system.

4.5.1. Objective Function

In the governing system of the pumped storage units, although the PID controller parameters can be set on the turbine pump governor to verify the integrity and effectiveness of the proposed identification strategy, the parameter set to be identified includes the PID control parameters. Therefore, eight parameters are identified as follows: proportional parameter $K_p$, integral parameter $K_i$, derivative parameter $K_d$, servomotor time constant $T_y$, pipeline characteristic coefficient $h_w$, time constant of water hammer pressure wave $T_r$, generator inertia time constant $T_a$, and generator regulation coefficient $e_g$. In the model of the system, the square of the difference between the measurable output selected for the objective function and the reintegration within 50 s after the disturbance starts. The objective function is as follows:

$$C\left(\theta \right)=\sum _{i=1}^N{{(x}_i-\widehat{x_i})}^2+\sum _{i=1}^N{{(y}_i-\widehat{y_i})}^2+\sum _{i=1}^N{{(m}_t-\widehat{m_t})}^2$$

(23)

where N is the sample capacity, $x_i$, $y_i$, and $m_t$ are the speed, guide vane opening, and torque of the real system output; $\widehat{x_i}$, $\widehat{y_i}$, and $\widehat{m_t}$ are the speed, guide vane opening, and torque of the identification system output; $\theta =\left[K_p{K}_i{K}_d{T}_y{h}_w{T}_r{T}_a{e}_g\right]$ are the parameters to be identified.

4.5.2. Identification Strategy

The process of identification is illustrated in Figure 6. Firstly, the system is given an excitation signal to obtain the respective dynamic responses $[x,y,m_t]$. Secondly, the intelligent optimization algorithm is used to initialize the parameters and obtain the dynamic response of the identification system $\left[\widehat{x},\widehat{y},\widehat{m_t}\right]$, and then the fitness function is calculated. Finally, the fitness function is optimized according to the intelligent optimization algorithm to obtain the minimum value of the fitness function.

4.5.3. Analysis of the Simulation Test Results

To verify the accuracy of the improved artificial hummingbird algorithm for identifying the governing system in the pumped storage unit, this paper compares the following four optimization algorithms—the particle swarm optimization algorithm (PSO), ant lion optimization algorithm (ALO), gravitational search algorithm (GSA), artificial hummingbird algorithm (AHA). The governing system in the pumped storage unit established in the previous section was used as the object of study and parameter values $\left[K_p{K}_i{K}_d{T}_y{h}_w{T}_r{T}_a{e}_g\right]$ were set as the parameter values of the system, as shown in

Table 5. Algorithm identification results of 5% frequency disturbance.

$\mathit{\theta }$	Real Value	PSO		ALO		GSA		AHA		IAHA
		$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE
$K_p$	3.21	3.4725	0.0818	3.4089	0.0620	2.7837	0.1328	3.1871	0.0071	3.1822	0.0087
$K_i$	2.68	2.7284	0.0181	2.6924	0.0046	2.8566	0.0659	2.7007	0.0077	2.7007	0.0077
$K_d$	1.24	1.5805	0.2746	1.5197	0.2256	1.6001	0.2904	1.3074	0.0544	1.2964	0.0455
$T_y$	0.30	0.5283	0.761	0.5040	0.68	0.3615	0.205	0.3152	0.0507	0.3146	0.0487
$h_w$	1.00	1.0480	0.0480	1.0580	0.0580	0.8805	0.1195	0.9518	0.0482	0.9691	0.0309
$T_r$	1.50	1.5566	0.0377	1.5049	0.0033	1.6834	0.1223	1.5536	0.0357	1.5283	0.0189
$T_a$	8.86	9.1299	0.0305	8.8038	0.0063	7.7609	0.1241	8.8523	0.0009	8.8335	0.0030
$e_g$	1.50	1.4974	0.0017	1.5004	0.0003	1.5104	0.0069	1.4996	0.0003	1.5001	0.00007
APE		0.1567		0.1300		0.1334		0.0256		0.0204

. The output state variables for simulation under certain operating conditions were used as the output variables. In order to evaluate the accuracy of parameter identification, the parameter error (PE) and average parameter error (APE) were used to measure the accuracy of the model parameters. The calculation formulas of PE and APE are as follows:

$$PE=\frac{|{\theta }_i-\widehat{{\theta }_i}|}{{\theta }_i}\times 100\%i=1,2,3,\cdots ,m$$

(24)

$$APE=\frac{1}{m}\sum _{i=1}^m\frac{|{\theta }_i-\widehat{{\theta }_i}|}{{\theta }_i}\times 100\%$$

(25)

The experiment was carried out under the condition of frequency disturbance; the frequency of the unit step was applied in the stable state of the unit. The simulation duration was 50 s, the sampling period was 0.01 s, and the number of identification tests was 20. Because the identification of the optimization algorithm was random, the average value of the results was taken. The six parameters of the hydraulic turbine can be obtained as follows: $e_x=-1.925,e_y=0.7133,e_h=1.413,e_{qx}=-0.7,e_{qy}=0.5833,e_{qh}=0.8555$, according to the rotational speed–flow characteristic curve and rotational speed–torque characteristic curve [38].

Table 5,

Table 6. Algorithm identification results of 10% frequency disturbance.

$\mathit{\theta }$	Real Value	PSO		ALO		GSA		AHA		IAHA
		$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE
$K_p$	3.21	3.6584	0.1397	3.3515	0.0441	2.9985	0.0659	3.1554	0.0170	3.2077	0.0007
$K_i$	2.68	2.6507	0.0109	2.7020	0.0082	2.7975	0.0438	2.7050	0.0093	2.690	0.0049
$K_d$	1.24	1.3532	0.0913	1.5432	0.2445	1.6201	0.3065	1.3087	0.0554	1.2974	0.0463
$T_y$	0.30	0.5116	0.7053	0.4603	0.5343	0.3820	0.2733	0.3084	0.0280	0.3155	0.0517
$h_w$	1.00	0.9294	0.0706	0.9430	0.057	0.9297	0.0703	0.9657	0.0343	0.9769	0.0231
$T_r$	1.50	1.6058	0.0705	1.6508	0.1005	1.7389	0.1593	1.5303	0.0202	1.5269	0.0179
$T_a$	8.86	9.4549	0.0671	8.8340	0.0029	7.6352	0.1382	8.8043	0.0063	8.8533	0.0008
$e_g$	1.50	1.4962	0.0025	1.5001	0.00007	1.5108	0.0072	1.5009	0.0006	1.5007	0.0005
APE		0.1447		0.1240		0.1330		0.0214		0.0182

and

Table 7. Algorithm identification results of 15% frequency disturbance.

$\mathit{\theta }$	Real Value	PSO		ALO		GSA		AHA		IAHA
		$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE	$\widehat{\mathit{\theta }}$	PE
$K_p$	3.21	3.4168	0.0644	3.3238	0.0355	2.8265	0.1195	3.1421	0.0241	3.1651	0.0140
$K_i$	2.68	2.7092	0.0109	2.7251	0.0168	2.9182	0.0889	2.7138	0.0126	2.7067	0.0100
$K_d$	1.24	1.7219	0.3886	1.5115	0.2190	1.5179	0.2241	1.3120	0.0581	1.2893	0.0398
$T_y$	0.30	0.4489	0.4963	0.4500	0.500	0.3286	0.0953	0.3075	0.0250	0.3009	0.0030
$h_w$	1.00	1.1061	0.1061	1.0637	0.0637	0.9323	0.0677	0.9701	0.0299	0.9869	0.0131
$T_r$	1.50	1.4625	0.025	1.4911	0.0059	1.6037	0.0691	1.5165	0.0110	1.5032	0.0021
$T_a$	8.86	8.7925	0.0076	8.684	0.0199	7.8136	0.1181	8.7751	0.0096	8.7994	0.0068
$e_g$	1.50	1.4993	0.0005	1.4996	0.0003	1.5038	0.0025	1.5014	0.0009	1.4997	0.0002
APE		0.1406		0.1076		0.0982		0.0211		0.0128

show the results of the parameter identification under 5%, 10%, and 15% frequency disturbances by different optimization algorithms. From the PE index of the parameter error of each algorithm, the identification of individual parameters was not optimal compared to PSO, ALO, GSA, and AHA, but only next to the best. From the APE index of the average parameter identification accuracy, it can be seen that under the three frequency disturbance conditions, 15% of the APE was only 1.28%, and the IAHA algorithm was better than the other four algorithms. IAHA improves the identification accuracy of the AHA algorithm accordingly. This shows that the proposed Chebyshev chaotic map and Levy flight improvement strategy effectively improve the search and convergence of the algorithm.

Figure 7a–c, Figure 8a–c, and Figure 9a–c, respectively, show the comparison diagram of the dynamic response of the unit speed, torque output, and guide vane opening output of the original parameters and IAHA algorithm optimization identification parameters under different frequency disturbance conditions, and they were locally amplified. It can be seen from the figure that the IAHA identification system is consistent with the original system, with high accuracy. Moreover, d–f are the identification errors of different links. It can be seen that the error finally approaches zero.

Figure 10, Figure 11 and Figure 12 show the convergence curves of fitness functions of different optimization algorithms under 5%, 10%, and 15% frequency interferences, respectively. The fitness function of the PSO is always the largest in the three graphs. It falls into the local optimum when it iterates 100 times and converges ahead of time. When GSA and ALO iterate 120 times, the fitness function value tends to be stable. The fitness function of the AHA algorithm is second only to IAHA. The IAHA algorithm is always the smallest, and has a rapid downward trend in the end. It reflects that IAHA has the ability of global optimization.

Figure 13, Figure 14 and Figure 15 show the comparison charts of the dynamic response of speed, torque, and guide vane opening output of different algorithms under different frequency disturbance conditions. It can be seen intuitively from the figure that the consistency between the output curve identified by IAHA and the output curve of the real system is higher than that of the other four algorithms, indicating that the accuracy is also higher.

5. Conclusions

In this paper, an improved AHA (IAHA) is proposed. Compared with AHA, there are two improved strategies. First, the population initialization of the Chebyshev chaotic map was used to expand the search range and improve the accuracy. Second, Levy flight was added to guide foraging, which gives the algorithm better convergence and stability. In order to verify the optimization performance of the IAHA algorithm, the 23 standard functions, which include unimodal and multimodal functions, were evaluated. The calculated values in the statistical analyses were the mean value, standard deviation, optimal value, and the worst value. IAHA performs best in most functions. The practical application involves identifying the parameters of the governing system of pumped storage units and calculating the average value through 20 independent operation algorithms. The results show that the errors of the IAHA are only 2.04%, 1.82%, and 1.28% under 5%, 10%, and 15% frequency disturbance conditions. Compared with PSO, GSA, ALO, and AHA, the identification accuracy of the governing system is improved.