A Tourist Flow Prediction Model for Scenic Areas Based on Particle Swarm Optimization of Neural Network

In recent years, China has been expanding domestic demand and promoting the service industry. This is a mixed blessing for the further development of tourism [1-4]. Faced with changing market, emerging conflicts, and growing demand for high-quality management, the tourism industry in China needs to undertake several challenging tasks: optimizing the industrial structure, renovating the growth mode, and improving development quality. It is urgent for the industry to transform from extensive operation to precise operation, from expanding the number of tourists to the rational allocation of tourist resources, and from meeting the basic needs of tourists to providing high-quality tourism services [5-8]. To balance the tourist flow and streamline the resource allocation in scenic areas, it is of great necessity to build an effective prediction model for tourist flow in scenic areas. The prediction results are important reference for scenic areas to implement development and planning, perform maintenance and repair, and provide smart tourism services.

In the early days, the long- and mid-term tourist flows of scenic areas are mostly forecasted by classic time series prediction methods, including linear regression models, autoregressive moving average (ARMA) models, and autoregressive integrated moving average (ARIMA) models [9-13]. For instance, Lim et al. [14] constructed a seasonal ARMA model for tourist flow prediction, which captures the obvious seasonal features in tourist flow. Based on GIOWHA operator, Cenamor et al. [15] combined least squares support vector regression (LSSVR) and ARMA into a hybrid model, selected such three indices as expected number of tourists, per-capita consumption of citizens, and the number of overnight tourists as the inputs of the model, and proved the good prediction effect of the hybrid model. To reduce constraints and speed up computation, Kotiloglu et al. [16] integrated ARMA model with the improved grey Markov model to project the number of tourists and foreign exchange income of Mount Taishan. Inspired by multivariate time series analysis, Mir [17] established an ARMA model based on the relevance between tourist flow and Baidu index of tourist spots, and applied the model to predict the tourist flow of a section of the Great Wall in Beijing.

However, many scholars argued that the traditional time series prediction methods are not suitable for predicting tourist flow in scenic areas, because the tourist flow is affected by various external factors, namely, weather and holidays [18-22]. With the rapid development of artificial intelligence (AI), many highly adaptive and self-learning models have been introduced to tourist flow prediction, including artificial neural network (ANN), random forest (RF), and support vector regression (SVR). Del Vecchio et al. [23] optimized the least squares support vector machine (LSSVM) with a self-designed improved fruit fly optimization (FOA) algorithm, and applied the FOA-LSSVM model to forecast the daily tourist flow in West Lake Scenic Area. Wise and Heidari [24] coupled SVR model with particle swarm optimization (PSO) algorithm into a PSO-SVR algorithm based on seasonal adjustment index, and carried out single- and multi-step predictions of monthly tourist flow in Beijing from 2000 to 2013.

For two reasons, the existing tourist flow prediction models for scenic areas often suffer from overfitting or underfitting, and thus fall into the local optimum trap: the tourist flow in scenic areas is affected by various external factors, which carry complex nonlinear features; the choice of kernel function and free parameters need to be further improved. To make accurate prediction of tourist flow, this paper probes deep into and further optimizes the existing AI methods, and develops a tourist flow prediction model for scenic areas based on the PSO of neural network (NN).

The remainder of this paper is organized as follows: Section 2 sets up a system of factors affecting the tourist flow in scenic areas, and eliminates the factors with low relevance through grey correlation analysis (GCA); Section 3 optimizes the long short-term memory (LSTM) NN with adaptive PSO, and then builds up a tourist flow prediction model for scenic areas; Section 4 explains the prediction workflow of the proposed model; Section 5 verifies the effectiveness of our model; Section 6 puts forward the conclusions.

3.1 LSTM NN

LSTM NN adds three gating units to each neuron, namely, input gate, forget gate, and output gate. Figure 2 shows the internal structure of LSTM neurons. The forget gate enables LSTM NN to selectively memorize and store important information, overcoming the long-term dependence of traditional NN. Hence, LSTM NN is suitable for processing data series B mentioned in the above section.

2.png

Figure 2. The internal structure of a LSTM neuron

The forget gate determines with a certain probability which information to discard from the neuron state. Let O_t-1 be the output of the previous neuron, and I_t be the input of the current neuron in LSTM NN. Then, the forget gate will read O_t-1 and I_t, and output a number between 0 and 1. The output number will be assigned to the state of the previous neuron S_t-1, indicating the probability of discarding the neuron state. If the number is 1, the neuron state will be retained fully; if the number is 0, the neuron state will be completely discarded.

${{f}_{t}}=Sigmoid({{W}_{f}}\cdot [{{O}_{t-1}},{{I}_{t}}]+{{\varepsilon }_{f}})$     (6)

Unlike the forget gate, the sigmoid activation function of the input gate determines which information needs to be updated and added to the neuron. A state vector S_t is generated by the tanh function to replace the state of the previous neuron S_t-1:

$S_{t}=f_{t} \cdot S_{t-1}+u_{t} \cdot S_{t}$      (7)

${{u}_{t}}=Sigmoid({{W}_{i}}\cdot [{{O}_{t-1}},{{I}_{t}}]+{{\varepsilon }_{i}})$     (8)

$S_{t}=\tanh \left(W_{C} \cdot\left[O_{t-1}, I_{t}\right]+\varepsilon_{C}\right)$     (9)

The output gate determines which information needs to be filtered and which needs to be output. The neuron state is processed by the tanh function. The output number, which falls between -1 and 1, is multiplied with the output of the sigmoid activation function. The final output can be expressed as formula (11):

${{h}_{t}}=Sigmoid({{W}_{h}}[{{O}_{t-1}},{{I}_{t}}]+{{\varepsilon }_{h}})$     (10)

${{O}_{t}}={{h}_{t}}\cdot \tanh ({{S}_{t}})$      (11)

Through the above analysis, this paper builds a 4-layer LSTM NN, including an input layer, an LSTM layer, a fully-connected layer, and an output layer. Among them, the input layer has three dimensions: the number of samples, the time step, and features. The LSTM layer contains n neurons with three gated units (forget gate, input gate, and output gate). Both the LSTM layer and the fully-connected layer use random orthogonal matrices for weight initialization. The difference is that the fully-connected layer adopts a linear activation function to synthesize the features extracted from the LSTM layer. The data series B obtained through GRA were taken as new features, and input to the LSTM NN for training and prediction.

The constructed model was trained in three steps: acquire valuable information through forward propagation and compute the prediction result based on reasonably set parameters; calculate the loss value; constantly correct the learning parameters through backpropagation of the loss. The loss function can be described as:

$E=\frac{1}{2}\sum\limits_{k=1}^{K}{{{({{y}_{lk}}-{{{{b}''}}_{lk}})}^{2}}}$      (12)

where, y_lk is the output of the LSTM NN. To improve the distribution width of activation values in each layer, enhance the expressiveness of the network, and overcome vanishing gradients, batch normalization was added to each layer during the training, using the batches in the learning process, that is, the mean and variance of data series B, which were input in p-dimensional batches, were respectively calculated by:

$\frac{1}{p}\sum\limits_{k=1}^{p}{{{{{b}''}}_{lk}}}\to {{\eta }_{B}}$      (13)

$\sqrt{\frac{1}{p}\sum\limits_{k=1}^{p}{{{({{{{b}''}}_{lk}}-{{\eta }_{B}})}^{2}}}}\to {{\sigma }_{B}}$      (14)

Let γ be a small value. Then, the data series B can be normalized with mean of η_B and variance of σ_B:

$\frac{b_{l k}^{\prime \prime}-\eta_{B}}{\sqrt{\sigma_{B}^{2}+\gamma}} \rightarrow \hat{b_{l k}^{\prime \prime}}$     (15)

3.2 PSO

In the traditional PSO algorithm, N particles are randomly initialized. After t iterations, the i-th particle will move at the velocity of {v_i1(t), v_i2(t), …, v_iD(t)} in the D-dimensional search space. At this time, the particle position is {x_i1(t), x_i2(t), …, x_iD(t)}. Starting with the current velocity and position, the i-th particle will search for the best-known individual solution p_best and the best-known global solution g_best. Let v_id(t) be the velocity of the i-th particle at the t-th iteration in the d-th dimension subspace; r₁, and r₂ be any numbers from 0 to 1. Then, the velocity and position of particles in the swarm can be respectively updated by:

$\begin{align}& {{v}_{id}}(t+1)=\omega {{v}_{id}}(t)+{{c}_{1}}{{r}_{1}}\left[ {{p}_{best}}(t)-{{x}_{id}}(t) \right] \\& +{{c}_{2}}{{r}_{1}}\left[ {{g}_{best}}(t)-{{x}_{id}}(t) \right] \\\end{align}$      (16)

 ${{x}_{id}}(t+1)={{x}_{id}}(t)+{{v}_{id}}(t+1)$     (17)

where, c₁ and c₂ are acceleration factors; ω is the inertia weight; p_best is the best-know position of the i-th particle at the t-th iteration; g_best is the best-known position of the swarm at the t-th iteration.

In the traditional PSO algorithm, the particles may converge prematurely with the growing number of iterations t:

${{\sigma }_{PSO}}=\sqrt{{{\sum\limits_{i=1}^{N}{\left( \frac{{{\lambda }_{i}}-\bar{\lambda }}{\lambda } \right)}}^{2}}}\le \beta $     (18)

where, β is a preset threshold; σ_PSO is the variance of swarm fitness. To limit the variance, a normalization factor called fitness λ is introduced:

$\lambda=\left\{\begin{array}{l}

\max \left|\lambda_{i}-\bar{\lambda}\right|, \max \left|\lambda_{i}-\bar{\lambda}\right|>1 \\

1, \max \left|\lambda_{i}-\bar{\lambda}\right|<1

\end{array}\right.$     (19)

The particles need to be disturbed more intensely to prevent premature convergence. The exploration and search abilities of particles can be controlled by adjusting the inertia weight ω, such that the particles could get out of the local optimum trap. Let t_max be the maximum number of iterations for the swarm. Then, ω can be expressed as:

$\omega ={{\omega }_{\max }}-\frac{t\cdot ({{\omega }_{\max }}-{{\omega }_{\min }})}{{{t}_{\max }}}$    (20)

where, ω_min and ω_max are the minimum and maximum of inertia weight, respectively. Traditionally, ω is adjusted nonlinearly with the Gaussian function:

$\omega ={{\omega }_{\min }}+({{\omega }_{\max }}-{{\omega }_{\min }}\times \exp \left[ \frac{-{{t}^{2}}}{0.2\times {{t}_{\max }}{{)}^{2}}} \right]$     (21)

However, this gradual adjustment approach will weaken the optimization performance of the swarm. In the early iterations, the particles have relatively weak local search ability, and converge rather slowly. It is difficult for a particle to make detailed searches in the global best-known area. In the later iterations, the particles have relatively weak global search ability, failing to jump out of the local optimum trap. To solve the problems, the adjustment method for ω was optimized as follows, considering the difference between the variances of swarm fitness in adjacent iterations:

$\omega =\left\{ \begin{align}  & {{\omega }_{\max }}-({{\omega }_{\max }}-{{\omega }_{\min }})\times \frac{\ln \left[ \bar{\lambda }(t+1) \right]-\ln \left[ {{\lambda }_{i}}(t+1) \right]}{\ln \left[ \bar{\lambda }(t+1) \right]},\Delta {{\sigma }_{PSO}}>{{\lambda }_{i}}(t+1)\ge \bar{\lambda } \\ & {{\omega }_{\min }},{{\lambda }_{i}}(t+1)<\min (\begin{matrix}   {\bar{\lambda }} & \Delta {{\sigma }_{PSO}}  \\\end{matrix}) \\ & {{\omega }_{\max }},{{\lambda }_{i}}(t+1)>\max (\begin{matrix}   {\bar{\lambda }} & \Delta {{\sigma }_{PSO}}  \\\end{matrix}) \\ & {{\omega }_{\min }}-({{\omega }_{\max }}-{{\omega }_{\min }})\times \frac{\ln \left[ {{\lambda }_{i}}(t+1) \right]-n\left[ \bar{\lambda }(t+1) \right]}{\ln \left[ {{\lambda }_{i}}(t+1) \right]},\Delta {{\sigma }_{PSO}}<{{\lambda }_{i}}(t+1)<\bar{\lambda } \\\end{align} \right.$     (22)

where, λ_i(t+1) is the fitness of the i-th particle at the t+1-th iteration; λ is the mean fitness of the swarm at the t-th iteration. It can be seen that, if Δσ_PSO is relatively large, the particle distribution is dispersed, and ω needs to be reduced to accelerate the convergence; if λ_i is greater than`λ, the particle is outside the best-known area, and ω needs to be increased to relocate the optimal solution; if λ_i is smaller than λ, the particle has good learning ability, and ω needs to be further reduced to find the optimal solution quickly; if Δσ_PSO is extremely small, the particle swarm will converge prematurely, and the particles need to be disturbed, making the swarm more diverse.

4a.png

(a) 30 iterations

4b.png

(b) 150 iterations

Figure 4. The particle distribution in the search process

The particle distributions of the traditional PSO algorithm and the proposed adaptive PSO algorithm were recorded, aiming to demonstrate the ability of our algorithm to optimize NN. The particle distributions after 30 and 150 iterations are presented in Figures 4(a) and 4(b), respectively. Compared with the traditional PSO algorithm, the proposed algorithm expands the search range of particles in the early and later iterations, during the search for the optimal solution. In this way, our algorithm effectively prevents particles from falling into local optimum trap and avoids premature convergence. Moreover, the particles of our algorithm clustered around the optimal solution in the early phase of the search, indicating that our algorithm ensures fast and effective convergence.

To verify the performance of the proposed LSTM NN, this paper carries out simulations on Matlab2010, with the dataset of daily tourist flow in a scenic area from 2018 to 2019 as the training set, and the dataset of daily tourist flow in that scenic area from January to September 2020 as the test set. The original data series exist as a 3×1,578 matrix.

The convergence of loss function in the LSTM NN is displayed in Figure 5. Obviously, the loss values on test set and training set both tended to be stable at around 30 iterations, indicating that the network has basically converged. Thus, the maximum number of iterations was set to 40. The loss values on test set and training set were about 0.0342 and 0.0328, respectively.

Figure 6 provides the prediction results of our NN on the tourist flow in the scenic area. It can be seen that the predicted values basically overlapped with the actual values. The prediction accuracy was still desirable, even if the tourist flow surged up on holidays in May and October.

Furthermore, our model was compared with several time series prediction methods, including RF, SVM, recurrent neural network (RNN), the LSTM optimized by traditional PSO (PSO-LSTM). The prediction results and errors of these methods are given in Tables 2 and 3, respectively. As shown in Table 2, the peak season of the scenic area lasts from April to October, according to the prediction results of all contrastive methods. Hence, the predictions are basically consistent with the actual situation.

5.png

Figure 5. The convergence of loss function in the LSTM NN

6.png

Figure 6. The prediction results of our model

Table 2. The comparison of prediction results

Table 3. The prediction errors of different models

Note: RMSE, MAE, MAPE, and R² are short for root mean square error, mean absolute error, mean absolute percentage error, and coefficient of determination, respectively.

As shown in Table 3, the RF, which is based on random sampling, had the greatest prediction error. This means RF is not a good choice for the prediction of tourist flow with strong nonlinear features. The PSO-LSTM had lower errors than SVM and RNN, indicating that LSTM NN can effectively overcome the long-term dependence in time series data. Our model boasted the smallest errors, because the redundant information is removed through GRA of influencing factors. This dimensionality reduction treatment ensures the extraction of important information from the data series, while improving the training effect of the NN. The ultralow MAPE (<6.98) fully demonstrates the effectiveness of our model.