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

: Rural community population forecasting has important guiding signiﬁcance 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.


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.

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.

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 ≤ (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.

The Process of Initialization
The Aquila Optimizer uses the method of random initialization to generate the initial population. The specific expression is as follows: 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.

Expanded Exploration (X 1 )
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: where X 1 (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.
where D is the dimension of the problem and N is the population number.

Narrowed Exploration (X 2 )
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: where X 2 (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.
where s is a constant value fixed at 0.01, and u and v are random numbers between 0 and 1. σ is calculated by the following formula: where β 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: where: where r 1 takes a value between 1 and 20 to fix the number of search cycles, and U is a fixed value of 0.00565. D 1 is an integer from 1 to the length of the search space (D), and ω is a fixed value of 0.005.

Expanded Exploitation (X 3 )
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).
where X 3 (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.

Narrowed Exploitation (X 4 )
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: where X 4 (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 1 represents the various movements of the Aquila tracking prey, which is generated using Equation (14). G 2 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.
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 nearoptimal 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 pseudocode of the original AO algorithm is given in Algorithm 1 below.

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., δ, ω, etc.) Output: The optimal value fit(X best ) While (t < T) Calculate the fitness values of population, and record the best solution (X best ) Update the parameters such as x, y, QF, G 1 , Update the current solution based on Equation (2) When fit(X 1 (t + 1)) < fit(X(t)), replace X(t) by X 1 (t + 1) else Narrowed exploration (X 2 ): Update the current solution based on Equation (4) When fit(X 2 (t + 1)) < fit(X(t)), replace X(t) by X 2 (t + 1) end if else if rand ≤ 0.5

Quasi-Opposition Learning Strategy
Quasi-opposition learning (QOL) [35] strategy was proposed on the basis of oppositionbased 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: 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: where m = LB+UB 2 is the midpoint of the current search space, and r 1 and r 2 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.

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: where X i (t) (i = 1, 2, · · · , 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. σ is the wavelet mutation coefficient. Its calculation formula is as follows: where ψ(ϕ/α) = e −(ϕ/α) 2 /2 · cos(5ϕ/α) is the Morlet wavelet function, and 99% of its energy is concentrated between −2.5 and 2.5, so ϕ is the random number between −2.5α and 2.5α. The a is the scaling parameter, and its expression is: 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 new i is obtained, and greedy selection is made between the mutant individual X new i and the original individual X i , i.e.,: This process ensures that individuals with better fitness value enter the next iteration, thus improving the optimization ability and convergence speed of the algorithm.

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: Initialize the population X randomly and set parameters such as δ, ω, etc. Calculate quasi-opposition individual X q i of current individual X i based on Equation (17) Select the N individuals with better fitness value from X ∪ X q as the current population Judge whether the individual position is beyond the boundary Calculate the fitness values of population, and update the best solution (X best ) Update the parameters such as QF, G 1 , G 2 , etc.
Step1 Initialize population size N, variable dimension D, maximum iteration T, and parameters δ, ω, etc.; Step2 Initialize the population according to Equation (1); calculate quasi-opposition individual X q i of current individual X i based on Equation (17); and then select the N individuals with better fitness value from X ∪ 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 1 , and G 2 ; Step4 When t ≤ (2/3) × T, perform the exploration operation of the algorithm. If the random number rand ≤ 0.5, perform the expanded exploration and update the current solution by Equation (2); if the fitness value of X 1 (t + 1) is less than the fitness value of X(t), X 1 (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 2 (t + 1) is less than the fitness value of X(t), X 2 (t + 1) is used to update X(t); Step5 When t > (2/3) × T, perform the exploitation operation of the algorithm. If the random number rand ≤ 0.5, perform the expanded exploitation and update the current solution by Equation (11); if the fitness value of X 3 (t + 1) is less than the fitness value of X(t), X 3 (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 4 (t + 1) is less than the fitness value of X(t), X 4 (t + 1) is used to update X(t); Step6 When rand ≤ 0.5, execute wavelet mutation according to Equation (18); 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.

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)).

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 [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 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 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. 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    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, 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, Figures 1-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.

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. Tables 5-7     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.

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: then, its consistent r fractional order accumulative sequence is defined as: where: The consistent r fractional order inverse accumulative sequence of Equation (22) is defined as: where: Let the original sequence X (0) 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: which is often called the whitening equation of the CFANGBM(1, 1, b, c). (27) by (x (r) (t)) −γ and obtain Equation (28).

Multiply both sides of Equation
Let y (r) (t) = (x (r) (t)) 1−γ . Then, Equation (28) can be transformed into: The constant variation method is used to solve Equation (29), and y (r) (t) can be obtained: Further, since y (r) (t) = (x (r) (t)) 1−γ , the x (r) (t) is given by: x (r) (t) = (x (r) (1)) From Equation (24), it can be deduced that x (r) (1) = x (0) (1). Therefore, the time response function of Equation (27) is: where k = 2, 3, · · · , n. When the fractional order r is given to the model parameters a, b, c, and γ, Equation (32) can fit (k ≤ n) and forecast (k > n) the original sequence. a, b, and c are obtained by leastsquare estimation, and the expression is as follows: where, B and Y are as follows, respectively: where: z (r) (k) = 1 2 (y (r) (k − 1) + y (r) (k)), k = 2, 3, · · · , n. Substitute y (r) (t) = (x (r) (t)) 1−γ into Equation (34). Then, B and Y are transformed into: By substituting the obtained a, b, and c into Equation (32),x (r) (k), k = 1, 2, · · · , 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): 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: subject to:

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), · · · , 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 ∪ 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 1 , and G 2 ; Step5 When t ≤ (2/3) * T, perform the exploration operation of the algorithm. If the random number rand ≤ 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) × T, perform the exploitation operation of the algorithm. If the random number rand ≤ 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 ≤ 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.

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.
Mean Absolute Error MAE 1 n ∑ n k=1 x (0) (k) − x (0) (k) Figure 6. The process of the IAO solving the China's rural population forecasting model. Table 9 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 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 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. Table 11. The predicted China's rural population based on different models (×10 4 ).

Year
Real  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. 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.

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.