Population Forecast of China’s Rural Community Based on CFANGBM and Improved Aquila Optimizer Algorithm

Abstract

Rural community population forecasting has important guiding significance to rural construction and development. In this study, a novel grey Bernoulli model combined with an improved Aquila Optimizer (IAO) was used to forecast rural community population in China. Firstly, this study improved the Aquila Optimizer by combining quasi-opposition learning strategy and wavelet mutation strategy, and proposed the new IAO algorithm. By comparing with other algorithms on CEC2017 test functions, the proposed IAO algorithm has the advantages of faster convergence speed and higher convergence accuracy. Secondly, based on the data of China’s rural community population from 1990 to 2019, a consistent fractional accumulation nonhomogeneous grey Bernoulli model called CFANGBM(1, 1, b, c) was established for rural population forecasting. The proposed IAO algorithm was used to optimize the parameters of the model, and then the rural population of China was predicted. Four error measures were used to evaluate the model, and by comparing with other forecasting models, the experimental results show that the proposed model had the smallest error between the forecasted value and the real value, which illustrates the effectiveness of using the IAO algorithm to solve CFANGBM(1, 1, b, c). At the end of this paper, the forecast data of China’s rural population from 2020 to 2024 are given for reference.

1. Introduction

Since the 1990s, China’s rural areas have experienced a drastic change, and the decline of rural areas is an indisputable objective fact. With the transfer of a large number of rural surplus laborers to cities, the trend of population urbanization is also taking place and developing in China [1]. Grasping the trend of rural population transfer and guiding the orderly transfer of rural populations is key to promoting urbanization construction [2]. The change in population affects the development of social economy to a certain extent, and is also the basis and reference for relevant government departments to formulate policies. Therefore, a more accurate forecasting of the development trend of rural population is of great significance to both rural revitalization strategy and rural economic development. This study aimed to establish a prediction model to more accurately predict China’s rural population.

Currently, available population forecasting models are generally divided into the following categories: logistic model, neuronet model, machine-learning-based model, and grey model, etc. At present, there is little research on the prediction of the rural population in China. Lv et al. [3] predicted and analyzed the total and structure of the rural population in Heilongjiang Province from 2011 to 2060 by constructing a revised Leslie rural population prediction model. Guan et al. [4] established an ARIMA model to forecast the rural population in China from 1970 to 2015, published in the open database of the World Bank. Many scholars have conducted related studies on population forecasting. Xuan et al. [5] obtained an improved forecasting model based on the logistic model, and the model was applied to forecast China’s total population in the next 30 years; Zhang et al. [6] proposed and investigated a weights and structure determination (WASD) neuronet model to forecast the UK population; Wang et al. [7] proposed a machine-learning-based method to forecast and analyze regional population objectively; Fernandes et al. [8] used a novel approach by combining Micro–Macro (MicMac) models into an agent-based perspective to simulate and forecast the Portuguese population; Gao et al. [9] predicted the population of Anhui Province in China based on a GM(1, 1) model; Zeng [10] built a differential equation model of population growth based on the grey forecasting model, and the model was applied to predict China’s population. The grey model uses the grey theory of “small sample, poor information”, and uses the accumulative operator, inverse accumulative operator, and least-squares estimation to establish the model, which is easier to solve. GM(1, 1) [11], as a classic grey prediction model, has good prediction ability even with only a small amount of data, and is widely used in the prediction of small sample data. Since the traditional GM(1, 1) model was proposed, many new grey models have been proposed successively. Zhou et al. [12] presented a trigonometric grey prediction model (TRGM) by combining the GM(1, 1) with the trigonometric residual modification technique for forecasting electricity demand. This approach helps to improve the forecasting accuracy of the GM(1, 1). Xie et al. [13] proposed a novel discrete grey forecasting model termed the DGM model. This model enhanced the tendency-catching ability and increased the prediction accuracy of the GM(1, 1) model. On the basis of the fractional grey model, Ma et al. [14] proposed a new fractional time-delay grey model (FTDGM) considering the time-delay effect, and applied it to the prediction of coal and natural gas consumption in Chongqing. In this study, a consistent fractional accumulation nonhomogeneous grey Bernoulli model named CFANGBM(1, 1, b, c) was established, which is also an improvement of GM(1, 1). This model improves the prediction accuracy by combining the non-equal weight property of the uniform fractional accumulation operator with the advantage of the non-uniform exponential model in parameters. As the CFANGBM(1, 1, b, c) has multiple model parameters and nonlinear characteristics, it is difficult to be solved by conventional methods. Therefore, this study adopted the swarm intelligence algorithms, popular in recent years, to solve this model.

Swarm intelligence algorithms typically mimic the collective behavior of different populations of organisms in nature. One of the most classic swarm intelligence algorithms is particle swarm optimization (PSO) [15], which is inspired by the collective behavior of birds and fish. Ant colony optimization (ACO) [16] is inspired by the foraging behavior of ants. In recent years, new swarm intelligence algorithms have been proposed, including the bat algorithm (BA) [17], moth-flame optimization algorithm (MFO) [18], whale optimization algorithm (WOA) [19], grey wolf optimizer (GWO) [20], firefly algorithm (FA) [21], grasshopper optimization algorithm (GOA) [22], Harris hawk optimization (HHO) [23], barnacles mating optimizer (BMO) [24], salp swarm algorithm (SSA) [25], manta rays foraging optimization (MRFO) [26], marine predators algorithm (MPA) [27], chimp optimization algorithm(CHOA) [28], slime mould algorithm (SMA) [29], side-blotched lizard algorithm (SBLA) [30], African vultures optimization algorithm (AVOA) [31], artificial gorilla troops optimizer (GTO) [32], and Aquila Optimizer (AO) [33], etc. The specific development process of swarm intelligence algorithms can be found in the review article of Brezočnik et al. [34]. Aquila Optimizer is a new meta-heuristic optimization algorithm proposed by Abualigah et al., which is inspired by Aquila’s hunting behavior. In the process of experiment, the AO algorithm also suffers the defects of slow convergence rate and ease of falling into local optimum. Therefore, quasi-opposition learning strategy [35] and the wavelet mutation strategy [36] were introduced to improve the Aquila Optimizer, and an improved Aquila Optimizer (IAO) with higher computational accuracy and faster convergence speed was proposed. In this study, the proposed improved Aquila Optimizer algorithm was used to solve the Chinese rural population forecasting model (CFANGBM(1, 1, b, c)). In addition, the CFANGBM(1, 1, b, c) was compared with other forecasting models and used to predict China’s rural population in the next five years. To sum up, the main contributions of this study can be summarized as follows:

• In this study, an improved Aquila Optimizer (namely, IAO) was proposed, which combines quasi-opposition learning and wavelet mutation strategy to improve the solution accuracy and convergence speed of the algorithm. The performance of the IAO was tested on the CEC2017 test set.
• A consistent fractional accumulation nonhomogeneous grey Bernoulli model named the CFANGBM(1, 1, b, c) for predicting rural population in China was established. The proposed IAO algorithm was used to solve the model parameters. The fitting error of the CFANGBM(1, 1, b, c) on population data was compared with other grey prediction models: GM(1, 1), DGM(1, 1), TRGM, and FTDGM. The rural population of China in 2020–2024 was forecast.

The rest of the paper is arranged as follows. In Section 2, a novel improved AO algorithm is proposed, and its specific algorithm steps are given. Section 3 verifies the superiority of IAO by comparing it with eight other intelligent algorithms on CEC2017 test functions. Section 4 uses the proposed IAO algorithm to solve the established CFANGBM(1, 1, b, c) model to fit and forecast the rural population in China. Section 5 makes a simple summary of the whole study.

2. Improved Aquila Optimizer

The AO algorithm is a new swarm-intelligence-based optimization method proposed by Abualigah et al. [33] in 2021. This method is inspired by the behavior of the Aquila in the process of capturing prey in nature. From the statistical results of benchmark functions, the AO algorithm showed great performance in searching the optimal solution. However, similar to other classical optimization algorithms, there are still some opportunities to improve the ability of AO algorithm in exploration and exploitation.

2.1. Aquila Optimizer

The Aquila has four hunting methods: selecting the search space by high soar with the vertical stoop, exploring within a divergent search space by contour flight with short glide attack, exploiting within a convergent search space by low flight with slow descent attack, and swooping by walk and grabbing prey. These four methods correspond to the expanded exploration, narrowed exploration, expanded exploitation, and narrowed exploitation stages of the AO algorithm, respectively. The AO algorithm can use different behaviors to transfer from the exploration step to the exploitation step. When $t\le (2/3)T$, the exploration step will be stimulated; otherwise, the exploitation step will be implemented. The mathematical model of the AO algorithm is given below.

2.1.1. The Process of Initialization

The Aquila Optimizer uses the method of random initialization to generate the initial population. The specific expression is as follows:

$$X_{ij}=L{B}_j+rand\times (U{B}_j-L{B}_j), i=1,2,\cdots ,N, j=1,2,\cdots ,D,$$

(1)

where rand is a random number between 0 and 1, N represents the population size, and D represents the problem dimension. UB_j and LB_j, respectively, represent the upper and lower bounds of the search agent in the j-th dimension search space.

2.1.2. Expanded Exploration (X₁)

In the expanded exploration phase, the Aquila uses the first method to hunt prey, namely, the Aquila identifies the prey area and selects the best hunting area by high soar with the vertical stoop. The AO widely explores from high soar to determine the area of search space where the prey is located. The mathematical expression of this behavior is as follows:

$$X_1(t+1)=X_{best}(t)\times (1-\frac{t}{T})+(X_M(t)-X_{best}(t)\times rand),$$

(2)

where X₁(t + 1) is the solution of the next iteration of t generated by the first search method. X_best(t) is the best solution obtained thus far, which reflects the approximate position of prey. $1-t/T$ is used to control the extended exploration by the number of iterations. X_M(t) represents the average position of all current solutions during the t-th iteration, which is calculated by Equation (3). rand is a random value between 0 and 1. t is the current number of iterations and T is the maximum number of iterations.

$$X_M(t)=\frac{1}{N}{\displaystyle \sum _{i=1}^N{X}_i(t)}, \forall j=1,2,\cdots ,D,$$

(3)

where D is the dimension of the problem and N is the population number.

2.1.3. Narrowed Exploration (X₂)

At this phase, the second hunting method is performed. This method is called contour flight of a short glide attack. When the prey area is spotted from a high altitude, the Aquila hovers over the target prey, prepares to land, and then attacks. The Aquila narrowly explores a selected area of the prey in preparation for an attack. This behavior is expressed mathematically as follows:

$$X_2(t+1)=X_{best}\left(t\right)\times Levy\left(D\right)+X_R\left(t\right)+(y-x)\times rand,$$

(4)

where X₂(t + 1) is the solution of the next iteration of t generated by the second search method. D is the dimension of the problem. Levy(D) is the Levy flight distribution function, calculated by Equation (5). X_R(t) is a random solution chosen from N solutions in the t-th iteration.

$$Levy(D)=s\times \frac{u\times \sigma }{{\left|v\right|}^{\frac{1}{\beta }}},$$

(5)

where s is a constant value fixed at 0.01, and u and v are random numbers between 0 and 1. $\sigma$ is calculated by the following formula:

$$\sigma =\left(\frac{\Gamma (1+\beta )\times \sin (\frac{\pi \beta }{2})}{\Gamma (\frac{1+\beta }{2})\times \beta \times 2^{\frac{\beta -1}{2}}}\right),$$

(6)

where $\beta$ is a constant value fixed at 1.5.

In Equation (4), y and x are used to represent the spiral shape in the search process, and the calculation is as follows:

$$y=r\times \cos (\theta ), x=r\times \sin (\theta ),$$

(7)

where:

$$r=r_1+U\times D_1,$$

(8)

$$\theta =-\omega \times D_1+{\theta }_1,$$

(9)

$${\theta }_1=\frac{3\times \pi }{2},$$

(10)

where r₁ takes a value between 1 and 20 to fix the number of search cycles, and U is a fixed value of 0.00565. D₁ is an integer from 1 to the length of the search space (D), and ω is a fixed value of 0.005.

2.1.4. Expanded Exploitation (X₃)

During the extended exploitation phase, the Aquila uses the third method to hunt prey. At this time, the prey area is precisely designated and the Aquila is ready to land and attack. The Aquila vertically descends and makes an initial attack to see how the prey reacts. This approach is called a low flight with slow descent attack. Here, the Aquila uses a selected area of the target to approach prey and attack. This behavior is mathematically shown in Equation (11).

$$X_3(t+1)=\left(X_{best}\right(t)-X_M(t\left)\right)\times \alpha -rand+\left(\right(UB-LB)\times rand+LB)\times \delta ,$$

(11)

where X₃(t + 1) is the solution of the next iteration of t, which is generated by the third search method. rand is a random value between 0 and 1. α and δ are exploitation adjustment parameters fixed at 0.1. UB and LB represent the upper and lower bounds of the problem, respectively.

2.1.5. Narrowed Exploitation (X₄)

During the narrowed exploitation phase, the fourth hunting method is carried out, that is, when the Aquila approaches the prey, it attacks the prey on land based on its random movement. This method is called walking and grabbing. Finally, the Aquila attacks the prey in the last position. This behavior can be mathematically expressed as follows:

$$X_4(t+1)=QF\times X_{best}\left(t\right)-(G_1\times X(t)\times rand)-G_2\times Levy\left(D\right)+rand\times G_1,$$

(12)

where X₄(t + 1) is the solution of the fourth search method. QF represents the quality function used to balance the search strategy, which is calculated by Equation (13). G₁ represents the various movements of the Aquila tracking prey, which is generated using Equation (14). G₂ represents the decreasing value from 2 to 0, which is generated by Equation (15). X(t) is the current solution in the t-th iteration.

$$QF(t)=t^{\frac{2\times rand-1}{{(1-T)}^2}},$$

(13)

$$G_1=2\times rand-1,$$

(14)

$$G_2=2\times (1-\frac{t}{T}),$$

(15)

where QF(t) is the value of the quality function in the t-th iteration, and rand is a random value between 0 and 1.

In conclusion, the AO’s search strategy explores the reasonable position of near-optimal solution or optimal solution obtained thus far. Each solution updates its location based on the best solution obtained during the AO’s optimization process. To emphasize the balance between exploration and exploitation in AO, there are four different exploration and exploitation search strategies, that are the four stages described above. The pseudo-code of the original AO algorithm is given in Algorithm 1 below.

2.2. The Proposed Improved Aquila Optimizer

The core idea of the improved algorithm is to enhance the accuracy, keep the balance between global search and local exploitation ability, avoid premature convergence, and improve the convergence speed. Thus, two improvement strategies were introduced to the original AO algorithm.

Algorithm 1: Aquila Optimizer

Input: Aquila population X and related parameters (i.e., $\delta$, $\omega$, etc.)

Output: The optimal value fit(Xbest)

While (t < T)

Calculate the fitness values of population, and record the best solution (Xbest)

Update the parameters such as x, y, QF, G1, G2, etc.

if $t\le (\frac{2}{3})\ast T$

if $rand\le 0.5$

Expanded exploration (X1):

Update the current solution based on Equation (2)

When fit(X1(t + 1)) < fit(X(t)), replace X(t) by X1(t + 1)

else

Narrowed exploration (X2):

Update the current solution based on Equation (4)

When fit(X2(t + 1)) < fit(X(t)), replace X(t) by X2(t + 1)

end if

else

if $rand\le 0.5$

Expanded exploitation (X3):

Update the current solution based on Equation (11)

When fit(X3(t + 1)) < fit(X(t)), replace X(t) by X3(t + 1)

else

Narrowed exploitation (X4):

Update the current solution based on Equation (12)

When fit(X4(t + 1)) < fit(X(t)), replace X(t) by X4(t + 1)

end if

end if

t = t + 1

End While

2.2.1. Quasi-Opposition Learning Strategy

Quasi-opposition learning (QOL) [35] strategy was proposed on the basis of opposition-based learning (OBL) strategy, which can enhance the diversity of population, improve the quality of solution, and then improve the performance of the algorithm. Opposition-based learning strategy expands the search range by calculating the opposition solution of the current solution in the search space, with the formula of:

$$X_i^{\prime }=LB+UB-X_i$$

(16)

where X_i is the current solution, and UB and LB represent the upper and lower bounds of the current solution in the search space, respectively.

The quasi-opposition solution obtained by quasi-opposition learning strategy is located in the middle of the midpoint of the upper and lower bounds and the opposition solution of the current solution. At this time, the probability that the quasi-opposition solution is closer to the unknown optimal solution than the opposition solution is greater than 1/2, that is, the quasi-opposition solution is closer to the optimal solution than the opposition solution. The calculation formula of the quasi-opposition solution is:

$$X_i^q=\left\{\begin{array}{l}m+(m-X_i)\times r_1, if X_i<m\\ m-(X_i-m)\times r_2, else\end{array}\right.$$

(17)

where $m=\frac{LB+UB}{2}$ is the midpoint of the current search space, and r₁ and r₂ are random numbers between 0 and 1.

After randomly initializing the positions of N individuals, new N individuals are generated through quasi-opposition learning strategy, and then the fitness values of the current solution and quasi-opposition solution are calculated. The 2N individuals are sorted in ascending order according to the fitness values, and the first N individuals with better fitness are selected as the new population. Then, the exploration and exploitation processes of the AO algorithm are executed.

2.2.2. Wavelet Mutation Strategy

The quasi-opposition learning strategy can enhance the diversity of a population and improve the quality of the solution, but the problem that the algorithm easily converges prematurely and falls into the local optimum is not solved. Mutation is an important method to help the algorithm jump out of local optimum. At this stage, wavelet mutation strategy [36] was introduced to improve the performance of the AO algorithm. By setting the mutation probability P, each individual after the exploration and exploitation stage of the algorithm has the chance of mutation through wavelet mutation strategy. When $rand<P$, the individual performs Morlet wavelet mutation. The specific formula of mutation is:

$$X_i^{new}(t)=\left\{\begin{array}{l}{X}_i(t)+\sigma (UB-X_i(t)),\begin{array}{c}\end{array}rand<0.5\\ X_i(t)+\sigma (X_i(t)-LB),\begin{array}{c}\end{array}rand\ge 0.5\end{array}\right.$$

(18)

where $X_i(t) (i=1,2,\cdots ,N)$ is the position of the i-th individual in t-th generation, and UB and LB are the upper and lower bounds of the current search space, respectively. $\sigma$ is the wavelet mutation coefficient. Its calculation formula is as follows:

$$\sigma =\frac{1}{\sqrt{\alpha }}\psi (\frac{\phi }{\alpha })$$

(19)

where $\psi (\phi /\alpha )=e^{-{(\phi /\alpha )}^2/2}·\cos (5\phi /\alpha )$ is the Morlet wavelet function, and 99% of its energy is concentrated between −2.5 and 2.5, so $\phi$ is the random number between −2.5α and 2.5α. The a is the scaling parameter, and its expression is:

$$\alpha =s·{(\frac{1}{s})}^{(1-\frac{t}{t_{\max }})}$$

(20)

where s is a given constant.

By integrating wavelet mutation into the AO algorithm, the mutation space is dynamically adjusted by controlling the scaling parameter of the Morlet wavelet function, so as to improve the stability of solution. After the completion of the wavelet mutation operation, the mutant individual $X_i^{new}$ is obtained, and greedy selection is made between the mutant individual $X_i^{new}$ and the original individual $X_i$, i.e.,:

$$X_i(t+1)=\left\{\begin{array}{ll}{X}_i^{new}(t),& f(X_i^{new}(t))\le f(X_i(t))\\ X_i(t),& f(X_i^{new}(t))>f(X_i(t))\end{array}\right.$$

(21)

This process ensures that individuals with better fitness value enter the next iteration, thus improving the optimization ability and convergence speed of the algorithm.

2.2.3. Overview of Improved Aquila Optimizer

In this study, the original AO algorithm was improved by introducing quasi-opposition learning strategy and wavelet mutation strategy. Quasi-opposition learning strategy can enhance population diversity and global exploration ability, while wavelet mutation can improve the ability of the algorithm to jump out of the local optimum. The improved Aquila Optimizer is called IAO for short. The detailed steps of the IAO algorithm are as follows:

Algorithm 2: The Proposed IAO

Initialize the population X randomly and set parameters such as $\delta$, $\omega$, etc.

Calculate quasi-opposition individual $X_i^q$ of current individual Xi based on Equation (17)

Select the N individuals with better fitness value from $X\cup X^q$ as the current population

While (t < T)

Judge whether the individual position is beyond the boundary

Calculate the fitness values of population, and update the best solution (Xbest)

Update the parameters such as QF, G1, G2, etc.

for i = 1, 2, …, N do

if $t\le (\frac{2}{3})\times T$

if $\mathit{r}\mathit{a}\mathit{n}\mathit{d}\le 0.5$ //Expanded exploration

$X_1=X_{best}\times (1-\frac{t}{T})+\left(X_M\right(t)-X_{best}(t)\times rand),$

if fit(X1) < fit(X)

X1 = X

end if

else //Narrowed exploration

$X_2=X_{best}\times Levy(D)+X_R+(y-x)\times rand,$

if fit(X2) < fit(X)

X2 = X

end if

end if

else

if $rand\le 0.5$ //Expanded exploitation

$X_3=(X_{best}-X_M)\times \alpha -rand+((UB-LB)\times rand+LB)\times \delta ,$

if fit(X3) < fit(X)

X3 = X

end if

else //Narrowed exploitation

$X_4=QF\times X_{best}-(G_1\times X\times rand)-G_2\times Levy(D)+rand\times G_1,$

if fit(X4) < fit(X)

X4 = X

end if

end if

end if

if rand < 0.5 //wavelet mutation

$X_i^{new}=\left\{\begin{array}{l}{X}_i+\sigma (UB-X_i),\begin{array}{c}\end{array}rand<0.5\\ X_i+\sigma (X_i-LB),\begin{array}{c}\end{array}rand\ge 0.5\end{array}\right.$

if fit($X_i^{new}$) < fit($X_i$)

$X_i$ = $X_i^{new}$

end if

end if

end for

t = t + 1

End While

Return: The optimal value fit(Xbest)

Step1 Initialize population size N, variable dimension D, maximum iteration T, and parameters $\delta$, $\omega$, etc.;

Step2 Initialize the population according to Equation (1); calculate quasi-opposition individual $X_i^q$ of current individual X_i based on Equation (17); and then select the N individuals with better fitness value from $X\cup X^q$ as the current population;

Step3 When t < T, calculate the fitness value of each individual, record the best solution X_best, update the parameters x, y, QF, G₁, and G₂;

Step4 When $t\le (2/3)\times T$, perform the exploration operation of the algorithm. If the random number $rand\le 0.5$, perform the expanded exploration and update the current solution by Equation (2); if the fitness value of X₁(t + 1) is less than the fitness value of X(t), X₁(t + 1) is used to update X(t). If the random number $rand>0.5$, perform the narrowed exploration and update the current solution by Equation (4); if the fitness value of X₂(t + 1) is less than the fitness value of X(t), X₂(t + 1) is used to update X(t);

Step5 When $t>(2/3)\times T$, perform the exploitation operation of the algorithm. If the random number $rand\le 0.5$, perform the expanded exploitation and update the current solution by Equation (11); if the fitness value of X₃(t + 1) is less than the fitness value of X(t), X₃(t + 1) is used to update X(t). If the random number $rand>0.5$, perform the narrowed exploitation and update the current solution by Equation (12); if the fitness value of X₄(t + 1) is less than the fitness value of X(t), X₄(t + 1) is used to update X(t);

Step6 When $rand\le 0.5$, execute wavelet mutation according to Equation (18);

Step7 Let t = t + 1, if t < T, return Step3; otherwise, output the optimal value fit(X_best) and the optimal individual X_best.

In addition, Algorithm 2 gives the pseudo-code of IAO.

2.3. Computational Complexity of the Improved Aquila Optimizer

The computational complexity of IAO depends on three rules: initial solution, calculating fitness function, and updating solution. Let N be the number of solutions, T be the number of max iterations, and D be the problem dimension. Firstly, N solutions are randomly generated and their quasi-opposition solutions are calculated during the initialization of the population, so the computational complexity of the initialization process of the solution is O(N × D). In the original AO algorithm, the computational complexity of the fitness function of the solution is O(T × N), and that of the updated solution is O(T × N × D). As the computational complexity of the added wavelet mutation strategy is O(T ×N × D), the total computational complexity of the proposed IAO is O(N × (D + T + T × D)).

3. Comparison of Improved Aquila Optimizer with Other Algorithms

In this section, the validity of the proposed algorithm is verified by the CEC2017 test functions, which contain four types of functions: unimodal, multimodal, hybrid, and composite. The detailed description of these functions is presented in

Table 1. CEC2017 test functions.

Function	No.	Function Name	Optimal Value
Unimodal Functions	F1	Shifted and Rotated Bent Cigar Function	100
	F3	Shifted and Rotated Zakharov Function	300
Multimodal Functions	F4	Shifted and Rotated Rosenbrock’s Function	400
	F5	Shifted and Rotated Rastrigin’s Function	500
	F6	Shifted and Rotated Expanded Scaffer’s F6 Function	600
	F7	Shifted and Rotated Lunacek Bi_Rastrigin Function	700
	F8	Shifted and Rotated Non-Continuous Rastrigin’s Function	800
	F9	Shifted and Rotated Levy Function	900
	F10	Shifted and Rotated Schwefel’s Function	1000
Hybrid Functions	F11	Hybrid Function of Zakharov, Rosenbrock, and Rastrigin’s	1100
	F12	Hybrid Function of High Conditioned Elliptic, Modified Schwefel, and Bent Cigar	1200
	F13	Hybrid Function of Bent Cigar, Rosenbrock, and Lunacek Bi_Rastrigin	1300
	F14	Hybrid Function of Elliptic, Ackley, Schaffer, and Rastrigin	1400
	F15	Hybrid Function of Bent Cigar, HGBat, Rastrigin, and Rosenbrock	1500
	F16	Hybrid Function of Expanded Schaffer, HGBat, Rosenbrock, and Modified Schwefel	1600
	F17	Hybrid Function of Katsuura, Ackley, Expanded Griewank plus Rosenbrock, Modified Schwefel, and Rastrigin	1700
	F18	Hybrid Function of high conditioned Elliptic, Ackley, Rastrigin, HGBat, and Discus	1800
	F19	Hybrid Function of Bent Cigar, Rastrigin, Expanded Griewank plus Rosenbrock, Weierstrass, and expanded Schaffer	1900
	F20	Hybrid Function of HappyCat, Katsuura, Ackley, Rastrigin, Modified Schwefel, and Schaffer	2000
Composition Functions	F21	Composition Function of Rosenbrock, High Conditioned Elliptic, and Rastrigin	2100
	F22	Composition Function of Rastrigin, Griewank, and Modified Schwefel	2200
	F23	Composition Function of Rosenbrock, Ackley, Modified Schwefel, and Rastrigin	2300
	F24	Composition Function of Ackley, High Conditioned Elliptic, Griewank, and Rastrigin	2400
	F25	Composition Function of Rastrigin, HappyCat, Ackley, Discus, and Rosenbrock	2500
	F26	Composition Function of Expanded Scaffer, Modified Schwefel, Griewank, Rosenbrock, and Rastrigin	2600
	F27	Composition Function of HGBat, Rastrigin, Modified Schwefel, Bent-Cigar, High Conditioned Elliptic, and Expanded Scaffer	2700
	F28	Composition Function of Ackley, Griewank, Discus, Rosenbrock, HappyCat, and Expanded Scaffer	2800
	F29	Composition Function of shifted and rotated Rastrigin, Expanded Scaffer, and Lunacek Bi_Rastrigin	2900
	F30	Composition Function of shifted and rotated Rastrigin, Non-Continuous Rastrigin, and Levy Function	3000
Search Range: [−100, 100]D

[37], and the specific expression of the functions’ formula is shown in [38]. In order to evaluate the advantages and disadvantages of the proposed IAO, it was compared with the original Aquila Optimizer (AO) [33], Archimedes optimization algorithm (AOA) [39], particle swarm optimization (PSO) [15], Harris hawks optimization (HHO) [23], sine–cosine algorithm (SCA) [40], whale optimization algorithm (WOA) [19], ant colony optimization (ACO) [16], and Genetic algorithm (GA) [41], which are rather classical and have been widely used in practical engineering optimization problems. Related specific parameters of these algorithms were set as

Table 2. The parameter settings of all comparison algorithms.

Algorithm	Parameter Value
AO	Exploitation adjustment parameters: α 0.1, δ = 0.1
AOA	Constant parameters: C1 = 2, C2 = 6, C3 = 1, C4 = 2
PSO	Velocity range: 0.5 times the variable range
	Cognitive and social factors: c1 = 2, c2 = 2.5
HHO	The value of E0 = [−1,1]
SCA	Constant value $a$ = 2
WOA	$\overrightarrow{a}$ variable linearly decreases from 2 to 0
ACO	Pheromone evaporation coefficient: Rho = 0.95 Pheromones increase intensity: Q = 1 Speed of ant: Lambda = 0.5
GA	Crossover probability: Pc = 0.95 Mutation probability: Pm = 0.001 Select probability: Er = 0.2
IAO	Mutation probability is 0.5, constant s = 10,000

shows. To obtain an unbiased experimental result, all the experiments were carried out on the same computer, the algorithms run on Windows 10(64-bit), Intel(R) Core(TM) i7-1165G7 (Santa Clara, CA, USA), CPU @2.80GHz (Santa Clara, CA, USA), MATLAB 2020a (Natick, MA, USA). In order to compare algorithms more fairly, the population size of all algorithms was 50, dimension D = 10, and the maximum number of iterations T = 500.

Table 3. Results of the comparative algorithms on CEC2017 test functions.

	Result	Algorithms
		IAO	AO	AOA	PSO	HHO	SCA	WOA	ACO	GA
F1	Best	2.05E+05	6.55E+05	3.83E+09	2.91E+06	1.48E+05	3.20E+06	2.38E+06	1.62E+09	1.43E+09
	Mean	4.25E+05	2.48E+06	7.97E+09	3.62E+07	5.32E+05	6.07E+07	1.46E+07	3.99E+09	2.61E+09
	Worst	7.26E+05	1.04E+07	1.45E+10	7.68E+07	1.52E+06	4.06E+08	8.05E+07	7.45E+09	4.46E+09
	Std	1.52E+05	2.10E+06	2.85E+09	2.03E+07	3.39E+05	1.19E+08	1.91E+07	1.58E+09	7.32E+08
	Rank	1	3	9	5	2	6	4	8	7
F3	Best	3.01E+02	3.26E+02	1.49E+03	5.23E+02	3.04E+02	7.66E+03	6.87E+02	1.07E+05	1.17E+04
	Mean	3.08E+02	1.19E+03	3.94E+03	1.64E+03	3.54E+02	5.27E+04	3.49E+03	4.17E+05	2.95E+05
	Worst	3.28E+02	3.05E+03	9.15E+03	2.78E+03	5.77E+02	2.03E+05	1.46E+04	1.08E+06	1.83E+06
	Std	8.61E+00	7.11E+02	1.98E+03	6.00E+02	7.12E+01	4.90E+04	4.23E+03	2.23E+05	5.91E+05
	Rank	1	3	6	4	2	7	5	9	8
F4	Best	4.03E+02	4.04E+02	5.81E+02	4.07E+02	4.00E+02	4.05E+02	4.06E+02	5.13E+02	5.40E+02
	Mean	4.05E+02	4.32E+02	7.79E+02	4.16E+02	4.19E+02	4.14E+02	4.56E+02	7.33E+02	7.67E+02
	Worst	4.08E+02	4.85E+02	1.10E+03	4.89E+02	4.89E+02	4.47E+02	5.78E+02	1.26E+03	1.23E+03
	Std	1.20E+00	3.46E+01	1.60E+02	1.80E+01	2.69E+01	9.81E+00	5.67E+01	1.95E+02	1.72E+02
	Rank	1	5	9	3	4	2	6	7	8
F5	Best	5.05E+02	5.12E+02	5.58E+02	5.22E+02	5.17E+02	5.16E+02	5.30E+02	5.63E+02	5.20E+02
	Mean	5.13E+02	5.31E+02	5.70E+02	5.31E+02	5.52E+02	5.27E+02	5.51E+02	5.85E+02	5.32E+02
	Worst	5.28E+02	5.46E+02	5.81E+02	5.44E+02	5.83E+02	5.40E+02	6.02E+02	6.04E+02	5.53E+02
	Std	6.44E+00	9.69E+00	7.18E+00	5.35E+00	1.66E+01	8.28E+00	1.81E+01	1.08E+01	7.74E+00
	Rank	1	4	8	3	7	2	6	9	5
F6	Best	6.00E+02	6.07E+02	6.29E+02	6.02E+02	6.14E+02	6.02E+02	6.11E+02	6.34E+02	6.02E+02
	Mean	6.02E+02	6.17E+02	6.37E+02	6.04E+02	6.38E+02	6.05E+02	6.35E+02	6.49E+02	6.08E+02
	Worst	6.06E+02	6.37E+02	6.47E+02	6.06E+02	6.60E+02	6.11E+02	6.70E+02	6.67E+02	6.19E+02
	Std	1.63E+00	7.39E+00	5.63E+00	1.26E+00	1.29E+01	2.47E+00	1.57E+01	9.59E+00	4.85E+00
	Rank	1	5	7	2	8	3	6	9	4
F7	Best	7.21E+02	7.34E+02	7.58E+02	7.42E+02	7.42E+02	7.26E+02	7.31E+02	8.32E+02	7.27E+02
	Mean	7.34E+02	7.57E+02	7.82E+02	7.62E+02	7.84E+02	7.38E+02	7.77E+02	9.05E+02	7.40E+02
	Worst	7.59E+02	7.86E+02	8.16E+02	7.87E+02	8.27E+02	7.55E+02	8.29E+02	9.64E+02	7.59E+02
	Std	1.15E+01	1.46E+01	1.50E+01	1.42E+01	2.55E+01	8.63E+00	2.51E+01	3.19E+01	8.44E+00
	Rank	1	4	7	5	8	2	6	9	3
F8	Best	8.07E+02	8.11E+02	8.26E+02	8.24E+02	8.13E+02	8.11E+02	8.17E+02	8.72E+02	8.12E+02
	Mean	8.18E+02	8.23E+02	8.39E+02	8.34E+02	8.28E+02	8.20E+02	8.39E+02	8.93E+02	8.20E+02
	Worst	8.40E+02	8.36E+02	8.50E+02	8.51E+02	8.42E+02	8.27E+02	8.70E+02	9.06E+02	8.29E+02
	Std	9.02E+00	6.83E+00	6.85E+00	7.08E+00	9.24E+00	4.72E+00	1.47E+01	8.24E+00	4.68E+00
	Rank	1	4	8	6	5	2	7	9	3
F9	Best	9.00E+02	9.17E+02	9.98E+02	9.12E+02	1.02E+03	9.02E+02	1.00E+03	1.80E+03	1.02E+03
	Mean	9.00E+02	1.03E+03	1.19E+03	9.55E+02	1.50E+03	9.39E+02	1.52E+03	2.48E+03	1.38E+03
	Worst	9.01E+02	1.26E+03	1.39E+03	1.04E+03	1.90E+03	1.12E+03	3.72E+03	3.53E+03	2.31E+03
	Std	2.50E-01	9.56E+01	9.51E+01	3.83E+01	2.45E+02	5.16E+01	6.03E+02	5.27E+02	3.34E+02
	Rank	1	4	5	3	7	2	8	9	6
F10	Best	1.24E+03	1.53E+03	1.99E+03	1.70E+03	1.48E+03	1.50E+03	1.32E+03	2.25E+03	1.55E+03
	Mean	1.52E+03	1.97E+03	2.50E+03	2.26E+03	2.03E+03	2.15E+03	2.06E+03	2.46E+03	1.97E+03
	Worst	1.82E+03	2.55E+03	2.86E+03	2.60E+03	2.46E+03	2.89E+03	2.81E+03	3.32E+03	2.52E+03
	Std	1.93E+02	2.73E+02	2.15E+02	2.60E+02	2.91E+02	3.80E+02	4.46E+02	1.28E+03	2.70E+02
	Rank	1	2	9	7	4	6	5	8	3
F11	Best	1.11E+03	1.13E+03	1.20E+03	1.12E+03	1.13E+03	1.12E+03	1.13E+03	3.89E+03	2.62E+03
	Mean	1.13E+03	1.18E+03	1.63E+03	1.14E+03	1.19E+03	1.16E+03	1.23E+03	9.97E+03	1.76E+04
	Worst	1.16E+03	1.34E+03	3.44E+03	1.36E+03	1.31E+03	1.33E+03	1.52E+03	2.55E+04	6.86E+04
	Std	1.50E+01	6.16E+01	5.66E+02	5.20E+01	5.18E+01	4.66E+01	9.11E+01	5.53E+03	1.48E+04
	Rank	1	4	7	2	5	3	6	8	9
F12	Best	1.70E+04	8.54E+04	4.25E+05	2.71E+05	8.13E+03	2.10E+04	8.48E+04	2.53E+07	2.99E+07
	Mean	3.92E+05	3.83E+06	1.21E+07	2.20E+06	3.53E+06	1.82E+06	6.73E+06	1.83E+08	2.37E+08
	Worst	2.50E+06	1.69E+07	5.72E+07	8.49E+06	1.31E+07	1.49E+07	1.92E+07	4.76E+08	4.78E+08
	Std	6.11E+05	4.42E+06	1.64E+07	2.03E+06	3.40E+06	3.67E+06	5.98E+06	1.23E+08	1.30E+08
	Rank	1	5	7	3	4	2	6	8	9
F13	Best	2.34E+03	3.60E+03	8.60E+03	1.69E+03	1.83E+03	5.30E+03	1.84E+03	1.69E+05	1.56E+06
	Mean	7.70E+03	1.69E+04	8.92E+04	9.69E+03	2.39E+04	9.98E+04	1.58E+04	3.49E+06	6.33E+07
	Worst	1.98E+04	4.86E+04	9.36E+05	2.78E+04	8.49E+04	2.98E+05	3.89E+04	9.28E+06	2.98E+08
	Std	4.31E+03	1.23E+04	2.09E+05	8.43E+03	1.78E+04	7.66E+04	1.19E+04	2.60E+06	6.74E+07
	Rank	1	4	6	2	5	7	3	8	9
F14	Best	1.47E+03	1.49E+03	1.46E+03	1.43E+03	1.47E+03	1.72E+03	1.47E+03	1.61E+03	2.40E+03
	Mean	1.53E+03	2.72E+03	1.85E+03	1.57E+03	1.59E+03	8.68E+03	1.88E+03	2.42E+03	6.93E+05
	Worst	1.75E+03	7.25E+03	5.26E+03	2.44E+03	1.90E+03	3.55E+04	2.58E+03	4.48E+03	3.17E+06
	Std	6.32E+01	1.40E+03	8.61E+02	2.27E+02	1.16E+02	8.83E+03	3.71E+02	8.62E+02	9.86E+05
	Rank	1	7	4	2	3	8	5	6	9
F15	Best	1.62E+03	2.06E+03	1.70E+03	1.58E+03	1.72E+03	8.67E+03	1.63E+03	3.09E+03	2.51E+04
	Mean	3.84E+03	9.12E+03	3.42E+03	2.16E+03	6.15E+03	5.11E+04	7.96E+03	9.81E+03	2.06E+06
	Worst	7.53E+03	3.43E+04	7.43E+03	3.36E+03	1.08E+04	4.04E+05	1.57E+04	2.68E+04	1.84E+07
	Std	2.02E+03	6.82E+03	1.65E+03	5.40E+02	3.23E+03	9.75E+04	3.97E+03	6.31E+03	4.07E+06
	Rank	3	6	2	1	4	8	5	7	9
F16	Best	1.60E+03	1.64E+03	1.73E+03	1.60E+03	1.77E+03	1.63E+03	1.65E+03	1.83E+03	2.05E+03
	Mean	1.70E+03	1.80E+03	1.95E+03	1.66E+03	1.96E+03	1.76E+03	1.92E+03	2.19E+03	2.37E+03
	Worst	1.97E+03	2.04E+03	2.17E+03	1.76E+03	2.15E+03	1.93E+03	2.22E+03	2.51E+03	2.76E+03
	Std	1.32E+02	1.19E+02	1.24E+02	6.15E+01	1.06E+02	9.19E+01	1.67E+02	1.73E+02	1.85E+02
	Rank	2	4	6	1	7	3	5	8	9
F17	Best	1.72E+03	1.73E+03	1.76E+03	1.73E+03	1.75E+03	1.74E+03	1.76E+03	1.88E+03	1.87E+03
	Mean	1.76E+03	1.77E+03	1.80E+03	1.77E+03	1.77E+03	1.77E+03	1.83E+03	2.02E+03	2.33E+03
	Worst	1.81E+03	1.85E+03	1.85E+03	1.86E+03	1.82E+03	1.89E+03	1.99E+03	2.25E+03	2.81E+03
	Std	2.05E+01	2.90E+01	2.42E+01	3.91E+01	2.16E+01	3.08E+01	6.87E+01	1.04E+02	2.61E+02
	Rank	1	5	6	2	4	3	7	8	9
F18	Best	2.88E+03	5.32E+03	3.29E+03	3.10E+03	4.27E+03	5.66E+03	2.26E+03	4.35E+05	8.26E+06
	Mean	1.09E+04	2.69E+04	3.83E+06	1.15E+04	1.64E+04	8.98E+04	1.58E+04	3.83E+06	3.12E+08
	Worst	2.29E+04	7.29E+04	4.55E+07	2.68E+04	3.76E+04	5.92E+05	3.90E+04	7.79E+06	8.05E+08
	Std	6.28E+03	1.50E+04	1.11E+07	8.09E+03	1.05E+04	1.31E+05	1.27E+04	2.11E+06	2.99E+08
	Rank	1	5	7	2	4	6	3	8	9
F19	Best	1.93E+03	2.13E+03	2.15E+03	1.92E+03	2.38E+03	2.04E+03	2.17E+03	1.22E+04	5.22E+04
	Mean	3.60E+03	2.24E+04	4.80E+04	3.00E+03	1.11E+04	2.24E+04	7.51E+04	2.44E+05	1.88E+07
	Worst	1.59E+04	2.22E+05	4.18E+05	9.67E+03	3.81E+04	1.43E+05	4.69E+05	1.38E+06	6.09E+07
	Std	3.14E+03	4.78E+04	9.84E+04	1.92E+03	1.06E+04	3.40E+04	1.21E+05	3.32E+05	1.88E+07
	Rank	2	4	6	1	3	5	7	8	9
F20	Best	2.02E+03	2.04E+03	2.09E+03	2.02E+03	2.05E+03	2.04E+03	2.07E+03	2.14E+03	2.01E+03
	Mean	2.04E+03	2.12E+03	2.14E+03	2.05E+03	2.17E+03	2.08E+03	2.18E+03	2.25E+03	2.06E+03
	Worst	2.06E+03	2.20E+03	2.24E+03	2.19E+03	2.34E+03	2.20E+03	2.34E+03	2.33E+03	2.16E+03
	Std	1.18E+01	4.80E+01	4.13E+01	4.39E+01	7.71E+01	4.46E+01	6.90E+01	5.59E+01	5.07E+01
	Rank	1	5	6	2	7	4	8	9	3
F21	Best	2.20E+03	2.20E+03	2.22E+03	2.20E+03	2.20E+03	2.21E+03	2.22E+03	2.37E+03	2.21E+03
	Mean	2.26E+03	2.32E+03	2.32E+03	2.32E+03	2.32E+03	2.30E+03	2.33E+03	2.39E+03	2.29E+03
	Worst	2.33E+03	2.35E+03	2.38E+03	2.34E+03	2.39E+03	2.33E+03	2.38E+03	2.40E+03	2.35E+03
	Std	5.88E+01	3.71E+01	5.62E+01	3.94E+01	6.13E+01	4.28E+01	4.92E+01	7.09E+00	5.00E+01
	Rank	1	5	4	6	7	3	8	9	2
F22	Best	2.30E+03	2.27E+03	2.35E+03	2.31E+03	2.31E+03	2.30E+03	2.31E+03	2.31E+03	2.31E+03
	Mean	2.31E+03	2.31E+03	2.54E+03	2.32E+03	2.32E+03	2.73E+03	2.32E+03	3.52E+03	2.32E+03
	Worst	2.31E+03	2.32E+03	2.88E+03	2.33E+03	2.33E+03	4.20E+03	2.34E+03	4.69E+03	2.34E+03
	Std	1.22E+00	1.11E+01	1.44E+02	5.09E+00	5.05E+00	6.30E+02	9.05E+00	1.81E+03	8.82E+00
	Rank	1	2	7	4	3	8	6	9	5
F23	Best	2.61E+03	2.62E+03	2.67E+03	2.62E+03	2.62E+03	2.62E+03	2.62E+03	2.70E+03	2.62E+03
	Mean	2.61E+03	2.65E+03	2.71E+03	2.63E+03	2.67E+03	2.63E+03	2.66E+03	2.78E+03	2.66E+03
	Worst	2.62E+03	2.67E+03	2.75E+03	2.64E+03	2.71E+03	2.66E+03	2.72E+03	2.86E+03	2.74E+03
	Std	4.33E+00	1.49E+01	2.35E+01	7.11E+00	3.01E+01	9.60E+00	2.87E+01	3.39E+01	2.89E+01
	Rank	1	4	8	2	7	3	5	9	6
F24	Best	2.50E+03	2.50E+03	2.59E+03	2.53E+03	2.50E+03	2.74E+03	2.56E+03	2.74E+03	2.54E+03
	Mean	2.68E+03	2.75E+03	2.77E+03	2.74E+03	2.81E+03	2.76E+03	2.79E+03	2.84E+03	2.71E+03
	Worst	2.75E+03	2.81E+03	2.89E+03	2.76E+03	2.91E+03	2.78E+03	2.83E+03	2.87E+03	2.83E+03
	Std	1.07E+02	8.65E+01	8.92E+01	6.78E+01	8.02E+01	1.03E+01	5.67E+01	2.59E+01	1.13E+02
	Rank	1	4	6	3	8	5	7	9	2
F25	Best	2.90E+03	2.90E+03	3.06E+03	2.90E+03	2.90E+03	2.91E+03	2.90E+03	3.00E+03	2.92E+03
	Mean	2.91E+03	2.93E+03	3.22E+03	2.94E+03	2.94E+03	2.93E+03	2.95E+03	3.17E+03	2.96E+03
	Worst	2.94E+03	2.96E+03	3.43E+03	2.95E+03	3.03E+03	2.95E+03	3.03E+03	3.30E+03	3.03E+03
	Std	2.12E+01	2.28E+01	1.03E+02	1.98E+01	2.95E+01	1.09E+01	3.13E+01	7.29E+01	2.34E+01
	Rank	1	3	9	4	5	2	6	8	7
F26	Best	2.82E+03	2.62E+03	3.27E+03	2.92E+03	2.82E+03	2.84E+03	2.96E+03	3.33E+03	2.87E+03
	Mean	2.89E+03	3.06E+03	3.58E+03	3.00E+03	3.25E+03	3.17E+03	3.34E+03	3.70E+03	3.25E+03
	Worst	2.95E+03	3.75E+03	3.93E+03	4.13E+03	4.15E+03	3.47E+03	4.46E+03	4.09E+03	4.11E+03
	Std	3.13E+01	2.16E+02	1.93E+02	2.66E+02	3.68E+02	1.76E+02	4.36E+02	1.81E+02	3.19E+02
	Rank	1	3	8	2	5	4	7	9	6
F27	Best	3.09E+03	3.09E+03	3.14E+03	3.10E+03	3.10E+03	3.08E+03	3.10E+03	3.15E+03	3.11E+03
	Mean	3.09E+03	3.11E+03	3.21E+03	3.11E+03	3.17E+03	3.18E+03	3.14E+03	3.18E+03	3.16E+03
	Worst	3.10E+03	3.12E+03	3.32E+03	3.17E+03	3.27E+03	3.20E+03	3.22E+03	3.21E+03	3.23E+03
	Std	2.37E+00	7.78E+00	3.54E+01	1.50E+01	4.08E+01	3.92E+01	4.37E+01	1.73E+01	3.24E+01
	Rank	1	2	9	3	6	8	4	7	5
F28	Best	3.10E+03	3.10E+03	3.28E+03	3.13E+03	3.10E+03	3.27E+03	3.18E+03	3.47E+03	3.12E+03
	Mean	3.26E+03	3.37E+03	3.39E+03	3.31E+03	3.38E+03	3.29E+03	3.38E+03	3.54E+03	3.48E+03
	Worst	3.41E+03	3.48E+03	3.77E+03	3.45E+03	3.75E+03	3.30E+03	3.74E+03	3.63E+03	3.72E+03
	Std	1.50E+02	1.18E+02	1.74E+02	1.34E+02	1.79E+02	9.88E+00	1.37E+02	4.45E+01	2.02E+02
	Rank	1	4	7	3	5	2	6	9	8
F29	Best	3.15E+03	3.18E+03	3.26E+03	3.17E+03	3.24E+03	3.15E+03	3.22E+03	3.35E+03	3.55E+03
	Mean	3.18E+03	3.26E+03	3.34E+03	3.21E+03	3.35E+03	3.26E+03	3.34E+03	3.46E+03	3.90E+03
	Worst	3.24E+03	3.38E+03	3.49E+03	3.35E+03	3.52E+03	3.36E+03	3.53E+03	3.56E+03	4.36E+03
	Std	2.17E+01	5.39E+01	5.67E+01	4.11E+01	7.74E+01	4.75E+01	8.17E+01	6.86E+01	2.06E+02
	Rank	1	4	5	2	7	3	6	8	9
F30	Best	6.82E+03	4.65E+03	1.12E+04	2.56E+04	7.54E+03	3.26E+03	1.67E+04	4.94E+05	1.19E+07
	Mean	1.37E+05	1.37E+06	1.46E+06	1.07E+06	1.28E+06	1.53E+04	1.51E+06	1.39E+06	6.63E+07
	Worst	1.45E+06	3.66E+06	4.30E+06	2.95E+06	5.13E+06	7.42E+04	7.17E+06	2.20E+06	1.42E+08
	Std	3.59E+05	1.28E+06	1.19E+06	1.14E+06	1.48E+06	1.77E+04	1.81E+06	5.67E+05	3.75E+07
	Rank	2	5	7	3	4	1	8	6	9
Mean Rank		1.1724	4.1034	6.7241	3.0345	5.1724	4.1379	5.8966	8.2069	6.5517
Runtime(s)		0.4300	0.2980	0.1404	0.0551	0.3092	0.0848	0.1146	1.3843	15.2474

shows the results of IAO and other optimization algorithms running 20 times for solving the 10-dimensional CEC2017 test set, including the best value, the mean value, the worst value, and the standard deviation. Meanwhile, according to the mean value on each function, the ranking of different algorithms was given to display the position of the IAO in all compared algorithms. As can be seen from Table 3, on the unimodal functions F1 and F3, the IAO algorithm was significantly better than other comparison algorithms, mainly because of the introduction of quasi-opposition learning strategy, which improves the solution accuracy of the algorithm. Comparing with the original population of the AO, the new initial population included more solutions closer to the optimum solution by selecting N elite individuals from the 2N individuals. That is, higher quality of the initial population ensured the accuracy of the algorithm on the unimodal functions. Thus, the IAO’s results ranked first among the nine comparison algorithms on F1 and F3.

For the multimodal functions F4–F10 with lots of local optimum, the wavelet mutation strategy was adopted to improve the ability of IAO to jump out of local optimum, which makes the algorithm more likely to find the global optimal solution. IAO ranked first on all multimodal functions, which confirmed the above hypothesis. Other traditional algorithms such as GA and ACO performed less well because of the disadvantage of falling into the local optimum easily.

For the hybrid functions F11–F20 and composite functions F21–F30, the algorithm must converge to the optimal solution in a more complex search space quickly, which is a challenge to the comprehensive performance of algorithms, and the IAO combining the quasi-opposition and wavelet mutation strategy achieved a better performance on F11–F30. Expect F15, F16, F19, and F30, the mean values of the IAO ranked first on all test functions. Overall, the IAO algorithm ranked first on 26 test functions, indicating that the IAO algorithm significantly outperformed other optimization algorithms on the CEC2017 test set. That is, by applying quasi-opposition learning and wavelet mutation strategy to the original AO algorithm, the IAO kept a better balance between exploration and exploitation, so as to improve the calculation accuracy of the original algorithm.

In the last two rows of Table 3, the average ranking and average running time are given. It can be seen that, for solving the CEC2017 test problem, the AO algorithm ranks third, and the IAO algorithm performs best. At the same time, the standard deviations of the IAO algorithm on 17 test functions were the smallest, which shows that the IAO algorithm was more stable in solving CEC2017 test functions. For the time cost, the time cost of IAO algorithm increased slightly due to the introduction of improvement strategies. The time cost of the IAO ranks seventh and the original AO ranks fifth. Thus, compared with the original AO algorithm, the slight increase in the cost of time is acceptable.

Moreover, the Wilcoxon rank-sum test was used to illustrate whether there was significant difference between the IAO and other compared algorithms. The statistical data are summarized in

Table 4. Wilcoxon rank-sum test p-values of comparison algorithms.

	AO	AOA	PSO	HHO	SCA	WOA	ACO	GA
F1	9.17E-08/+	6.80E-08/+	6.80E-08/+	4.25E-01/=	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+
F3	9.17E-08/+	6.80E-08/+	6.80E-08/+	9.28E-05/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+
F4	8.36E-04/+	6.80E-08/+	9.17E-08/+	6.01E-02/=	9.13E-07/+	1.20E-06/+	6.80E-08/+	6.80E-08/+
F5	3.07E-06/+	6.80E-08/+	2.22E-07/+	1.66E-07/+	1.41E-05/+	6.80E-08/+	6.80E-08/+	2.56E-07/+
F6	6.80E-08/+	6.80E-08/+	4.68E-05/+	6.80E-08/+	1.41E-05/+	6.80E-08/+	6.80E-08/+	2.69E-06/+
F7	1.41E-05/+	7.90E-08/+	1.80E-06/+	3.94E-07/+	1.26E-01/=	2.06E-06/+	6.80E-08/+	3.85E-02/+
F8	1.33E-02/+	1.20E-06/+	8.60E-06/+	2.34E-03/+	1.26E-01/=	1.81E-05/+	6.80E-08/+	1.14E-01/=
F9	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	6.80E-08/+
F10	1.60E-05/+	6.80E-08/+	2.22E-07/+	2.69E-06/+	1.80E-06/+	1.29E-04/+	1.22E-03/+	3.50E-06/+
F11	2.75E-04/+	6.80E-08/+	5.79E-01/=	1.10E-05/+	3.34E-03/+	2.06E-06/+	6.80E-08/+	6.80E-08/+
F12	5.87E-06/+	2.96E-07/+	1.25E-05/+	7.41E-05/+	2.75E-02/+	2.69E-06/+	6.80E-08/+	6.80E-08/+
F13	1.78E-03/+	3.07E-06/+	7.76E-01/=	7.41E-05/+	3.50E-06/+	6.79E-02/=	6.80E-08/+	6.80E-08/+
F14	1.20E-06/+	1.48E-01/=	1.64E-01/=	4.39E-02/+	7.90E-08/+	2.47E-04/+	1.92E-07/+	6.80E-08/+
F15	1.44E-04/+	7.56E-01/=	9.05E-03/−	2.23E-02/+	6.80E-08/+	3.75E-04/+	9.28E-05/+	6.80E-08/+
F16	4.70E-03/+	4.68E-05/+	3.51E-01/=	4.54E-06/+	1.33E-02/+	1.29E-04/+	1.92E-07/+	6.80E-08/+
F17	2.85E-01/=	3.99E-06/+	8.60E-01/=	9.62E-02/=	2.18E-01/=	2.36E-06/+	6.80E-08/+	6.80E-08/+
F18	5.90E-05/+	3.85E-02/+	9.46E-01/=	8.59E-02/=	9.75E-06/+	5.25E-01/=	6.80E-08/+	6.80E-08/+
F19	4.68E-05/+	9.28E-05/+	2.62E-01/=	8.36E-04/+	5.25E-05/+	2.36E-06/+	7.90E-08/+	6.80E-08/+
F20	2.06E-06/+	6.80E-08/+	8.82E-01/=	2.96E-07/+	1.79E-04/+	6.80E-08/+	6.80E-08/+	8.60E-01/=
F21	6.04E-03/+	1.04E-04/+	6.87E-04/+	7.41E-05/+	8.36E-04/+	1.67E-02/+	8.35E-04/+	2.56E-03/+
F22	2.00E-04/+	6.80E-08/+	6.80E-08/+	6.92E-07/+	4.60E-04/+	1.06E-07/+	1.22E-03/+	1.38E-06/+
F23	1.66E-07/+	6.80E-08/+	5.23E-07/+	6.80E-08/+	2.22E-07/+	1.92E-07/+	6.80E-08/+	1.43E-07/+
F24	1.41E-05/+	2.75E-02/+	6.67E-06/+	6.92E-07/+	6.01E-07/+	6.01E-07/+	1.23E-07/+	3.15E-02/+
F25	5.90E-05/+	6.80E-08/+	2.04E-05/+	4.68E-05/+	3.34E-03/+	4.54E-06/+	6.80E-08/+	9.13E-07/+
F26	8.29E-05/+	6.80E-08/+	6.92E-07/+	2.04E-05/+	9.13E-07/+	6.80E-08/+	6.80E-08/+	9.13E-07/+
F27	6.92E-07/+	6.80E-08/+	6.80E-08/+	6.80E-08/+	4.54E-06/+	6.80E-08/+	6.80E-08/+	6.80E-08/+
F28	1.12E-03/+	2.85E-01/=	4.11E-02/+	2.94E-02/+	1.00E+00/=	3.97E-03/+	6.80E-08/+	2.75E-04/+
F29	6.67E-06/+	6.80E-08/+	3.15E-02/+	7.90E-08/+	3.07E-06/+	7.90E-08/+	6.80E-08/+	6.80E-08/+
F30	8.29E-05/+	1.10E-05/+	5.17E-06/+	1.81E-05/+	2.80E-03/+	5.17E-06/+	5.23E-07/+	6.80E-08/+
+/=/−	28/1/0	26/3/0	20/8/1	25/4/0	25/4/0	27/2/0	29/0/0	27/2/0

, where ‘−’ represents that the IAO algorithm was inferior to the comparison algorithm at the 95% significance level (α = 0.05), ‘=’ represents that there was no significant difference between the IAO algorithm and the comparison algorithm, and ‘+’ indicates that the IAO algorithm was significantly better than the comparison algorithm. According to the statistical result of the original AO algorithm 28/1/0, there was significant difference between the IAO and AO in the results on 28 test functions, which shows that the IAO algorithm had significantly improved computational accuracy compared with the original algorithm. This difference is also evident in comparison with other algorithms. According to the statistical results of the last line in Table 4, PSO had a better effect among all the comparison algorithms. It was superior to IAO in 8 test functions, but inferior to IAO in 20 test functions. Therefore, overall, IAO had a significant advantage in at least 20 test functions.

In order to compare the performance of the algorithms more intuitively, Figure 1, Figure 2, Figure 3 and Figure 4 show the average convergence curves of each algorithm on CEC2017 test functions. Due to the introduction of quasi-opposition learning strategy, the proposed IAO produced more excellent individuals in the initial population stage. Therefore, the IAO was able to close to the theoretical optimal solution with quicker speed at the beginning of the search on unimodal functions F1 and F3. PSO and WOA also showed excellent performance at the beginning, but with the increase in the number of iterations, the performance of other algorithms’ convergence speed gradually deteriorated. In functions F5, F10, F12, F13, F21, F24, F26, F28, and F30, IAO continued to search quickly and efficiently in the later stage of the search, while the convergence speed of other comparison algorithms gradually slowed down. That is, the improvement strategies in IAO overcame the defect of premature convergence which caused the terrible performance of other algorithms, especially for the traditional GA and ACO algorithms.

Therefore, the IAO algorithm may show its advantages in parameter optimization problems, combinatorial optimization problems, or path optimization and other aspects. Excellent performance in convergence speed, accuracy, and stability may help it ascertain competitiveness in practical application. Meanwhile, the convergence speed of PSO and WOA were also competitive in the early iteration stage such as F1, F18, and F19. Therefore, it would be an interesting study how to learn from the initial search methods of these two algorithms to further improve the convergence speed of IAO.

4. The Solving of China’s Rural Community Population Forecast Model Using Improved Aquila Optimizer

In recent years, with the development of population urbanization in China, and a large number of urban and rural population flows, the rural community population shows a trend of decreasing year by year. The effective prediction of China’s rural population is of great significance to the rational planning of rural construction. In this section, a mathematical model is established to forecast the rural population.

Table 5. The rural community population in China from 1990 to 1999.

Year	1990	1991	1992	1993	1994	1995	1996	1997	1998	1999
Rural population (×104)	84,138	84,620	84,996	85,344	85,681	85,947	85,085	84,177	83,153	82,038

,

Table 6. The rural community population in China from 2000 to 2009.

Year	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009
Rural population (×104)	80,837	79,563	78,241	76,851	75,705	74,544	73,160	71,496	70,399	68,938

and

Table 7. The rural community population in China from 2010 to 2019.

Year	2010	2011	2012	2013	2014	2015	2016	2017	2018	2019
Rural population (×104)	67,113	65,656	64,222	62,961	61,866	60,346	58,973	57,661	56,401	55,162

show the rural population during 1990−1999, 2000−2009, and 2010−2019, respectively. In addition, when the rural population forecasting was carried out, the population of the first 25 years (1990−2014) was used to train the model, so as to obtain the specific parameter values of the model, and then the fitting effect of the model was tested with the population of the last 5 years (2015−2019). These rural community population data can be found in the “China Statistical Yearbook 2020” compiled by the National Bureau of Statistics of China. In the China Statistical Yearbook’ China’s total population can be divided into urban population and rural population according to urban and rural areas. This study was based on the prediction o’ China’s rural population.

In order to more intuitively observe the changing trend of the rural population, Figure 5 shows the bar chart of the rural population during 1990−2019. It can be seen that during 1990−1995, the rural population showed a slow growth trend, and after 1995, the rural population showed a declining trend year by year.

4.1. China’s Rural Population Forecasting Model

According to the changing trend of the rural community population, this paper used the consistent fractional accumulation nonhomogeneous grey Bernoulli model (called CFANGBM(1, 1, b, c) for short) to fit and forecast the above rural population. The model in this section is based on the consistent fractional accumulation and inverse accumulation operators.

Given the original data sequence:

$$X^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(n)),$$

(22)

then, its consistent r fractional order accumulative sequence is defined as:

$$X^{(r)}=(x^{(r)}(1),x^{(r)}(2),\cdots ,x^{(r)}(n)),$$

(23)

where:

$$x^{(r)}(k)={\displaystyle \sum _{i=1}^k\left(\begin{array}{l}k-i+[r]-1\\ k-i\end{array}\right)\frac{x^{(0)}(i)}{i^{[r]-r}}},r\in {\mathrm{R}}^{+},k=1,2,\cdots ,n,$$

(24)

The consistent $r$ fractional order inverse accumulative sequence of Equation (22) is defined as:

$$X^{(-r)}=(x^{(-r)}(1),x^{(-r)}(2),\cdots ,x^{(-r)}(n)),$$

(25)

where:

$$\left\{\begin{array}{l}{x}^{(-r)}(k)=k^{[r]-r}{\displaystyle \sum _{i=1}^k{(-1)}^{k-i}\left(\begin{array}{l}[r]\\ k-i\end{array}\right)x^{(0)}(i)},k=2,3,\cdots ,n,\\ x^{(-r)}(1)=x^{(0)}(1).\end{array}\right.$$

(26)

Let the original sequence X⁽⁰⁾ and its consistent $r$ fractional order accumulation sequence X^(r) be shown in Equations (22) and (23), respectively. Then, the CFANGBM(1, 1, b, c) is represented as the following differential equation:

$$\frac{\mathrm{d}{x}^{(r)}(t)}{\mathrm{d}t}+a{x}^{(r)}(t)=(bt+c){(x^{(r)}(t))}^{\gamma },\gamma \in \mathrm{R},$$

(27)

which is often called the whitening equation of the CFANGBM(1, 1, b, c).

Multiply both sides of Equation (27) by ${(x^{(r)}(t))}^{-\gamma }$ and obtain Equation (28).

$$\frac{\mathrm{d}{x}^{(r)}(t)}{\mathrm{d}t}{(x^{(r)}(t))}^{-\gamma }+a{(x^{(r)}(t))}^{1-\gamma }=bt+c,$$

(28)

Let $y^{(r)}(t)={(x^{(r)}(t))}^{1-\gamma }$. Then, Equation (28) can be transformed into:

$$\frac{\mathrm{d}{y}^{(r)}(t)}{\mathrm{d}t}+(1-\gamma )a{y}^{(r)}(t)=(1-\gamma )(bt+c),$$

(29)

The constant variation method is used to solve Equation (29), and $y^{(r)}(t)$ can be obtained:

$$y^{(r)}(t)=\left(y^{(r)}(1)+\frac{b}{a^2(1-\gamma )}-\frac{b+c}{a}\right)e^{-a(1-\gamma )(t-1)}+\frac{bt+c}{a}-\frac{b}{a^2(1-\gamma )},$$

(30)

Further, since $y^{(r)}(t)={(x^{(r)}(t))}^{1-\gamma }$, the $x^{(r)}(t)$ is given by:

$$x^{(r)}(t)={\left(\left({(x^{(r)}(1))}^{1-\gamma }+\frac{b}{a^2(1-\gamma )}-\frac{b+c}{a}\right)e^{-a(1-\gamma )(t-1)}+\frac{bt+c}{a}-\frac{b}{a^2(1-\gamma )}\right)}^{\frac{1}{1-\gamma }},$$

(31)

From Equation (24), it can be deduced that $x^{(r)}(1)=x^{(0)}(1)$. Therefore, the time response function of Equation (27) is:

$${\widehat{x}}^{(r)}(k)={\left[\left({(x^{(0)}(1))}^{1-\gamma }+\frac{b}{a^2(1-\gamma )}-\frac{b+c}{a}\right)e^{-a(1-\gamma )(k-1)}+\frac{bk+c}{a}-\frac{b}{a^2(1-\gamma )}\right]}^{\frac{1}{1-\gamma }},$$

(32)

where $k=2,3,\cdots ,n$.

When the fractional order $r$ is given to the model parameters a, b, c, and $\gamma$, Equation (32) can fit (k ≤ n) and forecast (k > n) the original sequence. a, b, and c are obtained by least-square estimation, and the expression is as follows:

$${(a,b,c)}^T={(B^{T}B)}^{-1}{B}^{T}Y,$$

(33)

where, B and Y are as follows, respectively:

$$B=(1-\gamma )\left(\begin{array}{ccc}{z}^{(r)}(2)& -\frac{3}{2}& -1\\ z^{(r)}(3)& -\frac{5}{2}& -1\\ ⋮& ⋮& ⋮\\ z^{(r)}(n)& -\frac{2n-1}{2}& -1\end{array}\right), Y=\left(\begin{array}{c}{y}^{(r)}(1)-y^{(r)}(2)\\ y^{(r)}(2)-y^{(r)}(3)\\ ⋮\\ y^{(r)}(n-1)-y^{(r)}(n)\end{array}\right),$$

(34)

where: $z^{(r)}(k)=\frac{1}{2}(y^{(r)}(k-1)+y^{(r)}(k)),k=2,3,\cdots ,n$.

Substitute $y^{(r)}(t)={(x^{(r)}(t))}^{1-\gamma }$ into Equation (34). Then, B and Y are transformed into:

$$B=(1-\gamma )\left(\begin{array}{ccc}\frac{1}{2}\left({(x^{(r)}(1))}^{1-\gamma }+{(x^{(r)}(2))}^{1-\gamma }\right)& -\frac{3}{2}& -1\\ \frac{1}{2}\left({(x^{(r)}(2))}^{1-\gamma }+{(x^{(r)}(3))}^{1-\gamma }\right)& -\frac{5}{2}& -1\\ ⋮& ⋮& ⋮\\ \frac{1}{2}\left({(x^{(r)}(n-1))}^{1-\gamma }+{(x^{(r)}(n))}^{1-\gamma }\right)& -\frac{2n-1}{2}& -1\end{array}\right),$$

(35)

$$Y=\left(\begin{array}{c}{(x^{(r)}(1))}^{1-\gamma }-{(x^{(r)}(2))}^{1-\gamma }\\ {(x^{(r)}(2))}^{1-\gamma }-{(x^{(r)}(3))}^{1-\gamma }\\ ⋮\\ {(x^{(r)}(n-1))}^{1-\gamma }-{(x^{(r)}(n))}^{1-\gamma }\end{array}\right),$$

(36)

By substituting the obtained a, b, and c into Equation (32), ${\widehat{x}}^{(r)}(k),k=1,2,\cdots ,n$ can be obtained, and the consistent fractional order inverse accumulation is performed on it. The predicted value of the CFANGBM(1, 1, b, c) model can be computed by Equation (37):

$$\left\{\begin{array}{l}{\widehat{x}}^{(0)}(k)=k^{[r]-r}{\displaystyle \sum _{i=1}^k{(-1)}^{k-i}\left(\begin{array}{l}[r]\\ k-i\end{array}\right){\widehat{x}}^{(r)}(i)},k=2,3,\cdots ,n,\\ {\widehat{x}}^{(0)}(1)={\widehat{x}}^{(r)}(1).\end{array}\right.$$

(37)

It can be seen that the key to prediction by using the CFANGBM(1, 1, b, c) model is to find the optimal parameters $r$ and γ, so this problem can be transformed into an optimization problem, and the optimization goal is to minimize the mean absolute percentage error between the predicted value and the real value. The objective function and constraint conditions of the optimization problem are:

$$\mathrm{minimize} \arg (f(r,\gamma ))=\frac{1}{n-1}{\displaystyle \sum _{k=2}^n\left|\frac{{\widehat{x}}^{(0)}(k)-x^{(0)}(k)}{x^{(0)}(k)}\right|}\times 100\%,$$

(38)

subject to:

$$\left\{\begin{array}{l}{(a,b,c)}^T={(B^{T}B)}^{-1}{B}^{T}Y,\\ {\widehat{x}}^{(r)}(k)={\left[\left({(x^{(0)}(1))}^{1-\gamma }+\frac{b}{a^2(1-\gamma )}-\frac{b+c}{a}\right)e^{-a(1-\gamma )(k-1)}+\frac{bk+c}{a}-\frac{b}{a^2(1-\gamma )}\right]}^{\frac{1}{1-\gamma }},\\ {\widehat{x}}^{(0)}(k)=k^{[r]-r}{\displaystyle \sum _{i=1}^k{(-1)}^{k-i}\left(\begin{array}{l}[r]\\ k-i\end{array}\right){\widehat{x}}^{(r)}(i)},k=2,3,\cdots ,n,\\ {\widehat{x}}^{(0)}(1)={\widehat{x}}^{(r)}(1).\end{array}\right.$$

(39)

4.2. The Steps of Improved Aquila Optimizer Solving China’s Rural Population Forecasting Model

As the above model is nonlinear, it is difficult to solve it using conventional methods. Therefore, this study used the IAO algorithm to optimize the parameters of the above China rural population forecasting model, so as to predict China’s population data from 2015 to 2019 and calculate the model error. The detailed steps are described as follows:

Step1 Input the original rural population data $X^{(0)}=(x^{(0)}(1),x^{(0)}(2),\cdots ,x^{(0)}(n))$;

Step2 Initialize the parameters of IAO, and initialize the individuals randomly by Equation (1);

Step3 Calculate quasi-opposition individual $X^q$ of current individual X based on Equation (17); and then take Equation (38) as the objective function, and select the N individuals with better fitness value from $X\cup X^q$ as the current population;

Step4 When t < T, calculate the fitness value of each individual, record the best solution X_best, and update the parameters x, y, QF, G₁, and G₂;

Step5 When $t\le (2/3)\ast T$, perform the exploration operation of the algorithm. If the random number $rand\le 0.5$, update the current solution by Equation (2); if the random number $rand>0.5$, update the current solution by Equation (4);

Step6 When $t>(2/3)\times T$, perform the exploitation operation of the algorithm. If the random number $rand\le 0.5$, update the current solution by Equation (11); if the random number $rand>0.5$, update the current solution by Equation (12);

Step7 When $rand\le 0.5$, execute wavelet mutation by Equation (18);

Step8 Let t = t + 1. If t < T, return Step4; otherwise, output the minimum error value fit(X_best) and the corresponding optimal parameter values r and γ.

Figure 6 shows the steps for the IAO to solve the China’s rural community population forecasting model.

4.3. Experimental Result Analysis of China’s Rural Population Forecasting Model

In this section, the CFANGBM(1, 1, b, c) model is compared with the other four models to verify the effectiveness of the proposed model in forecasting China’s rural population. These four comparison models are GM(1, 1) [11], DGM(1, 1) [13], TRGM [12], and FTDGM [14], respectively. In order to compare the advantages and disadvantages of the five forecasting models, the evaluation criteria for the advantages and disadvantages of the models are firstly given. The specific expressions of the four evaluation criteria are shown in

Table 8. The four evaluation criteria of the models.

Error	Abbreviation	Concrete Expression
Mean Absolute Percentage Error	MAPE	$\frac{1}{n-1}{\displaystyle {\sum }_{k=1}^n\left|\frac{{\widehat{x}}^{(0)}(k)-x^{(0)}(k)}{x^{(0)}(k)}\right|}\times 100\%$
Root Mean Square Percentage Error	RMSPE	$\sqrt{\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{\left(\frac{{\widehat{x}}^{(0)}(k)-x^{(0)}(k)}{x^{(0)}(k)}\right)}^2}}\times 100\%$
Mean Square Error	MSE	$\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{({\widehat{x}}^{(0)}(k)-x^{(0)}(k))}^2}$
Mean Absolute Error	MAE	$\frac{1}{n}{\displaystyle {\sum }_{k=1}^n\left|{\widehat{x}}^{(0)}(k)-x^{(0)}(k)\right|}$

.

Table 9. The fitted China’s rural population based on different models (×10⁴).

Year	Real Data	GM(1, 1)	DGM(1, 1)	TRGM	FTDGM	CFANGBM (1, 1, b, c)
1990	84,138	84,138.0	84,138.0	84,138.0	84,138.0	84,138.0
1991	84,620	89,704.7	89,712.5	90,599.5	85009.1	84,405.9
1992	84,996	88,419.9	88,427.0	89,696.5	85,484.0	85,157.3
1993	85,344	87,153.5	87,159.8	88,090.9	85,628.4	85,609.5
1994	85,681	85,905.2	85,910.8	86,194.6	85,496.6	85,681.0
1995	85,947	84,674.8	84,679.8	85,037.3	85,134.4	85,411.9
1996	85,085	83,462.0	83,466.3	84,206.5	84,579.8	84,863.5
1997	84,177	82,266.6	82,270.2	82,671.9	83,864.6	84,094.0
1998	83,153	81,088.3	81,091.3	80,845.7	83,015.3	83,153.1
1999	82,038	79,926.9	79,929.3	79,757.5	82,054.0	82,081.1
2000	80,837	78,782.1	78,783.9	78,994.7	80,999.3	80,910.3
2001	79,563	77,653.7	77,655.0	77,527.3	79,866.5	79,666.1
2002	78,241	76,541.5	76,542.2	75,767.2	78,668.7	78,368.2
2003	76,851	75,445.2	75,445.3	74,744.1	77,416.4	77,032.0
2004	75,705	74,364.6	74,364.2	74,045.6	76,118.7	75,669.4
2005	74,544	73,299.5	73,298.6	72,641.5	74,783.2	74,289.4
2006	73,160	72,249.6	72,248.2	70,943.8	73,416.0	72,898.8
2007	71,496	71,214.8	71,212.9	69,982.3	72,022.3	71,502.8
2008	70,399	70,194.8	70,192.5	69,344.5	70,606.7	70,105.1
2009	68,938	69,189.4	69,186.6	68,000.1	69,172.6	68,708.6
2010	67,113	68,198.4	68,195.2	66,361.4	67,723.1	67,315.0
2011	65,656	67,221.6	67,217.9	65,458.0	66,260.8	65,925.7
2012	64,222	66,258.8	66,254.7	64,877.3	64,787.8	64,541.6
2013	62,961	65,309.8	65,305.3	63,589.4	63,305.9	63,162.9
2014	61,866	64,374.4	64,369.5	62,006.3	61,816.5	61,790.0

displays the fitted China’s rural community population from 1990 to 2014 based on different models. Compared with the real data of rural population, the fitted error is given in

Table 10. The comparison of fitting errors of five forecasting models.

Error	GM(1, 1)	DGM(1, 1)	TRGM	FTDGM	CFANGBM (1, 1, b, c)
MAE	1.6148E+07	1.3233E+07	1.6143E+07	3.4565E+06	1.6636E+06
MAPE(%)	2.198626	1.738619	2.197725	0.477928	0.231867
MSE	3.7681E+14	2.6103E+14	3.7681E+14	1.6026E+13	4.3480E+12
RMSPE(%)	2.510064	2.001781	2.509359	0.535274	0.277700

to measure the performance of different forecasting methods. As can be seen from Table 10, the results of the CFANGBM(1, 1, b, c) based on the IAO algorithm ranked first on all four evaluation criteria, and its fitted data of rural population in China were closest to the real data, which was better than the fitting effect of other comparison models. Among them, the predicted rural population based on GM(1, 1), DGM(1, 1), and TRGM deviated greatly from the real data.

Figure 7 shows the comparison between the predicted data of the five models and the real data, from which it can be seen that the rural population reached the maximum value in 1995, but the GM(1, 1) and DGM(1, 1) methods had a great error in fitting the data at this stage. That is, the GM(1, 1) and DGM(1, 1) models would face great challenges with data with volatility characteristics, though they are effective for monotonically increasing or decreasing data. By adding triangulation correction, the TRGM method had better performance than GM(1, 1) in the later stage. In Figure 7c, the fitted data based on TRGM are closer to the real data, especially after 2009. However, for the whole stage from 1990 to 2014, it could not obtain better results by triangulation correction. Furthermore, the FTDGM method obtains more reasonable fitted data around 1995. The fitting curve of it was closer to the real curve than GM(1, 1), DGM(1, 1), and TRGM. Finally, in Figure 7e, the fitted data obtained by the CFANGBM(1, 1, b, c) method almost coincide with the real data, which is consistent with the results in Table 10. That is, by introducing the linear correction and searching the suitable parameters by the proposed IAO, the CFANGBM(1, 1, b, c) was able to fit data with different characteristics. Especially around the turning point of the population data, the parameters could be adjusted intelligently by the IAO algorithm to accommodate rapid changes in data. Thus, the approach in this study can achieve better fitting effects for more complex data.

The best parameters could be obtained by regarding the data from 1990 to 2014 as the training data. Then,

Table 11. The predicted China’s rural population based on different models (×10⁴).

Year	Real Data	GM(1, 1)	DGM(1, 1)	TRGM	FTDGM	CFANGBM (1, 1, b, c)
2015	60,346	63,452.4	63,447.1	61,157.7	60,320.9	60,422.7
2016	58,973	62,543.5	62,537.9	60,631.1	58,820.0	59,060.7
2017	57,661	61,647.7	61,641.7	59,396.5	57,314.8	57,703.8
2018	56,401	60,764.8	60,758.4	57,865.9	55,806.0	56,351.4
2019	55,162	59,894.4	59,887.8	57,069.0	54,294.2	55,003.1

shows China’s rural population predicted by the five prediction models from 2014 to 2019, which were compared with the real data. The forecasting errors of the five prediction models are listed in

Table 12. The comparison of forecasting errors of five forecasting models.

Error	GM(1, 1)	DGM(1, 1)	TRGM	FTDGM	CFANGBM (1, 1, b, c)
MAE	3.9520E+07	5.5452E+07	3.9460E+07	3.9743E+06	8.3145E+05
MAPE(%)	8.608105	10.303898	8.595082	0.882413	0.181523
MSE	1.5946E+15	3.2963E+15	1.5898E+15	2.5021E+13	8.6257E+11
RMSPE(%)	6.991762	9.755137	6.981294	0.896381	0.164034

.

As can be seen from Table 11, besides the CFANGBM(1, 1, b, c) method, the other four methods had poorer performance in predicted error than fitted error. The MAPE of GM(1, 1), DGM(1, 1), and TRGM reached about 10%, which means the predicted data obtained by these methods was unreasonable. It is obvious that the CFANGBM(1, 1, b, c) method based on the IAO was effective for China’s rural population forecasting with the smallest predicted error. Next, all the fitted and predicted data were used to plot Figure 8. In the enlarged sub-figure in Figure 8, the predicted data obtained by the CFANGBM(1, 1, b, c) method were closest to the real data, while other methods deviated from the real data.

Finally, the predicted data in the next five years (2020−2024) are given in

Table 13. The predicted data from 2020 to 2024 based on different models (×10⁴).

Year	GM(1, 1)	DGM(1, 1)	TRGM	FTDGM	CFANGBM (1, 1, b, c)
2020	59,036.6	59,029.6	56,532.6	52,779.8	53,658.2
2021	58,191.0	58,183.7	55,611.3	51,263.3	52,316.1
2022	57,357.6	57,349.9	54,705.0	49,745.1	50,976.3
2023	56,536.0	56,528.1	53,813.4	48,225.4	49,638.1
2024	55,726.3	55,718.1	52,936.4	46,704.4	48,300.7

. Figure 9 is the trend line of the predicted data from 2014 to 2024 based on the different models. On the right side of the black vertical line in Figure 9 are the predicted data of China’s rural population in the next five years (2020−2024). It is obvious that the trend of the red curve most closely matched the population data from 2014 to 2019. However, the trend curves based on other models were higher or lower, which is not suitable for the future prediction. By optimizing the parameters in the CFANGBM(1, 1, b, c) model through the proposed IAO, the method can be used to data forecast with higher precision more widely.

5. Conclusions

In this study, the proposed IAO algorithm was used to optimize the parameters of the forecasting model CFANGBM(1, 1, b, c), and then to predict China’s rural population. Firstly, an improved Aquila Optimizer with higher computational accuracy and faster convergence speed was proposed. Based on the original AO algorithm, the following two strategies were introduced to improve it: (1) quasi-opposition learning strategy, which enhances population diversity and improves the quality of solution; (2) wavelet mutation strategy, which enhances the ability of the algorithm to jump out of the local optimum, and improves the calculation accuracy of the algorithm. The advantages of the proposed IAO algorithm were verified by comparison with other intelligent algorithms on CEC2017 test functions.

Secondly, the proposed IAO was applied to the CFANGBM(1, 1, b, c) model which is used to forecast the rural population in China, and the model was compared with GM(1, 1), DGM(1, 1), TRGM and FTDGM. The experimental results showed that the CFANGBM(1, 1, b, c) model had smaller fitting and prediction errors compared with the other four models, and its fitting effect and prediction results were closer to the real value. This illustrates the superiority of the CFANGBM(1, 1, b, c) model in predicting the rural population, and also verifies the effectiveness of using the IAO algorithm to solve the model parameters. In addition, the forecast data of China’s rural population from 2020 to 2024 were given.

In future work, the improved algorithm will be considered to be applied to other various applications, such as image processing, text and data mining, feature selection, task scheduling, and others. Further, the CFANGBM(1, 1, b, c) model is also likely to be used for forecast of economic growth, climate changing, and resource consumption, etc.