An Outlier-Robust Point and Interval Forecasting System for Daily PM2.5 Concentration

Air pollution forecasting, particularly of PM2.5 levels, can be used not only to deliver effective warning information to the public but also to provide support for decisions regarding the control and treatment of air pollution problems. However, there are still some challenging issues in air pollution forecasting that urgently need to be solved, such as how to handle and model outliers, improve forecasting stability, and correct forecasting results. In this context, this study proposes an outlier-robust forecasting system to attempt to tackle the abovementioned issues and bridge the gap in the current research. Specifically, the system developed consists of two parts that deal with point and interval forecasting, respectively. For point forecasting, a data preprocessing module is proposed based on outlier handling and data decomposition to mitigate the negative influences of outliers and noise, which can also help the model capture the main characteristics of the original time series. Meanwhile, an outlier-robust forecasting module is designed for better modeling of the preprocessed data. For the model to further improve its accuracy, a nonlinear correction module based on an error ensemble strategy is developed that can provide more accurate forecasting results. Finally, the interval forecasting part of the system is based on a newly proposed artificial intelligence–based distribution evaluation and the results of the point forecasting part to present the range of future changes. Experimental results and analysis utilizing daily PM2.5 concentration from two provincial capital cities in China are discussed to verify the superiority and effectiveness of the system developed, which can be considered an effective technique for point and interval forecasting of daily PM2.5 concentration.

INTRODUCTION
Urbanization, industrialization, and energy consumption have caused the issue of air pollution to become increasingly serious. The air pollution issue is considered a major concern (Andrade et al., 2015) and regarded as the single largest health risk (Wendel, 2014). It can have adverse effects on human beings and bring great economic losses as well as problems for society, affecting areas such as public health (Li et al., 2017), corporate cash holdings (Li et al., 2021), and the tourism industry (Hao et al., 2021). In this context, to solve air pollution issues and accelerate ecological progress, air pollution forecasting, particularly of PM 2.5 (particulate matter with an aerodynamic diameter of 2.5 μm or less) levels, has been acknowledged as a promising technique for air pollution control and treatment (Liu et al., 2019a). However, the performance of air pollution forecasting suffers from many factors, and how to develop a model that can improve forecasting effectiveness is not only a challenging research topic but also a matter of growing concern. As a result, developing a more effective model for modeling PM 2.5 concentration has become an imperative task that cannot be postponed (Yang et al., 2019a).
Over the past few decades, many studies have been conducted to propose an effective air pollution forecasting model. In general, the current air pollution forecasting models mainly belong to three categories, namely, the chemical transport method (CTM), the traditional statistical method, and the artificial intelligence method. In the first category, Timmermans et al. (2017) employed a kind of CTM method to analyze PM 2.5 in China. Zhang et al. (2019) used the community multiscale air quality (CMAQ) method to analyze PM 2.5 pollution events in Qingdao, China. Similarly, Fan et al. (2015) employed third-generation air quality modeling system Models-3/CMAQ to analyze air pollution episodes in one region of China. However, the forecasting performance of CTM methods is easily influenced by the scale and quality of the emissions data (Feng et al., 2015). Furthermore, Stern et al. (2008) suggested that due to incomplete descriptions of physical and chemical processes and limited knowledge of pollution sources, CTM models may produce poor forecasting results. In the traditional statistical method category, the widely employed statistical model mainly consists of an autoregressive integrated moving average (ARIMA) and regression method. Vlachogianni et al. (2011) employed a regression model for forecasting nitrogen oxides (NOx) and PM 10 (particulate matter with an aerodynamic diameter of 10 μm or less). Zafra et al. (2017) developed an ARIMA model using hourly PM 10 concentration data. However, the traditional statistical model is unable to extract the complex and nonlinear features of pollutant concentration data, which may result in undesirable results (Wang et al., 2020a). Fortunately, with the rapid development of advanced technologies, new methods, including artificial neural networks (ANN), support vector regression (SVR), extreme learning machines (ELM), and other artificial intelligence methods, have been proposed and are being widely employed in different forecasting fields, including air pollution forecasting. These novel methods are acknowledged as promising solutions to air pollution forecasting issues (Yang et al., 2019b).
Although many hybrid models have been proposed for air pollution forecasting, there are still some challenging issues that urgently need to be solved. Specifically, as far as we are aware, most previous studies have employed data decomposition to improve forecasting performance while ignoring the significance of handling and modeling outliers in air pollution data, which may lead to the hybrid model being unable to further enhance the forecasting performance. Moreover, artificial intelligence optimization algorithms are used to search for the optimal parameters of methods in a hybrid model, but the forecasting ability can only be improved to a certain degree. Significantly, however, these studies ignore the timeconsuming issues caused by incorporating artificial intelligence optimization algorithms, which may be unable to completely overcome the model's limitations, i.e., the instability of the final results. Furthermore, most previous studies have emphasized the contribution of advanced data decomposition and optimization algorithms while ignoring the significance of mining the characteristics of the original air pollution time series and correcting forecasting results to further improve the model's forecasting performance, despite the growing importance of air pollution forecasting performance. Moreover, another issue with air pollution forecasting, especially daily air pollution forecasting, is that it is mainly focused on point forecasting and thus can only provide deterministic information that is insufficient for real application and cannot provide uncertainty information. As far as we know, interval forecasting can make up for the defects of point forecasting, but this method has been ignored by the relevant researchers despite it being a novel research area that is especially deserving of attention.
In this study, to develop an effective model that overcomes the abovementioned limitations of most previous studies, a novel outlier-robust point and interval forecasting system is proposed for forecasting daily PM 2.5 concentration, which is composed of two parts. The point forecasting part proposes a novel forecasting model based on a data preprocessing module, an outlier-robust forecasting module, and a nonlinear correction module to obtain future deterministic information, whereas the interval forecasting part is designed based on the newly proposed artificial intelligence-based distribution evaluation and point forecasting results to obtain future uncertainty information. Specifically, first, a data preprocessing module is proposed, which takes into consideration the significance of outlier handling and data decomposition, that can be employed to mitigate the negative effects of outlier and noise information so that the model developed can capture the main features and achieve better performance. Second, an outlier-robust forecasting module is designed, based on outlier-robust ELM (ORELM), to forecast preprocessed data. Next, a nonlinear correction module based on an error ensemble strategy is developed to mine information in the forecasting results and further improve the model's forecasting performance. Finally, an artificial intelligence-based distribution evaluation method is designed in the interval forecasting part, which can be combined with the point forecasting results to provide a range of future changes. Case studies utilizing daily PM 2.5 concentrations from two cities in China are designed to validate the developed system.
The main novelty and contribution of this study to current research can be summarized as follows: (1) The forecasting focus of the system developed is not only on point forecasting but also on interval forecasting. Previous daily PM 2.5 concentration forecasting studies mainly focused on point forecasting, which can provide deterministic information and is insufficient for real application, while ignoring the research about interval forecasting, which can provide uncertainty information. Therefore, an outlier-robust forecasting system that consists of point and interval forecasting is successfully proposed in this study, which can make up for the defects of point forecasting and is validated well in two cities in China.
(2) An improved data preprocessing module is designed to solve the outlier data and noise information issues simultaneously. Most previous studies only employed data decomposition for data preprocessing while ignoring the significance of outlier handling and modeling and failed to further enhance the forecasting performance. As a result, this study develops an improved data preprocessing module based on outlier handling and data decomposition that can effectively overcome the limitations caused by outlier and noise information.
(3) Point forecasting performance is further enhanced by proposing an outlier-robust forecasting module and a nonlinear correction module. In the daily PM 2.5 concentration forecasting field, outlier modeling and forecasting result correction are of great importance but have been ignored by the relevant researchers. Thus, on the one hand, the outlier-robust forecasting module is designed to further solve the outlier modeling issue; on the other hand, a nonlinear correction module is developed based on an error ensemble strategy to mine information in the forecasting results and further enhance the model's forecasting performance.
(4) Convincing experiments are designed to verify the effectiveness and superiority of the system developed. For system evaluation, five evaluation metrics are employed in the evaluation of point forecasting, whereas two typical metrics are used in the evaluation of interval forecasting. For model comparison, four experiments are designed to prove the effectiveness and superiority of outlier handling and modeling, data decomposition, nonlinear correction module, and the system developed in daily PM 2.5 concentration forecasting.
The methods are presented in the Methodology section, next the construction of the outlier-robust point and interval forecasting system are discussed, then the Experimental Analysis section presents the experiments, and the final section draws the conclusions of this study.

METHODOLOGY
This section introduces the detailed methods used in the system developed, which includes a data preprocessing module, an outlier-robust forecasting module, a nonlinear correction module, artificial intelligence-based distribution evaluation, and interval forecasting theory.

Data Preprocessing Module
The data preprocessing module is proposed on the basis of outlier handling and data decomposition, which can overcome the limitations caused by outlier and noise information.

Outlier Handling
Outliers in the original time series data will have a negative influence on the development of a model, which may bring poor results. Therefore, in this study, an outlier handling algorithm, the Hampel filter (HF), is introduced into the data preprocessing module. HF, developed by Hampel (1974), is robust against outliers (Liu et al., 2004). A brief explanation of HF can be defined as follows (Wu et al., 2021).
Given a sequence x 1 , x 2 , . . ., x n , W K i is a set of numbers within a moving window and m i is the median value from the moving window, which are defined as follows: where K is the sliding window's half-width.
The new data obtained after using HF to handle the original sequence can be defined as where t denotes a positive integer and S i denotes the median absolute deviation (MAD), which can be given by

Data Decomposition
Data decomposition, which has been acknowledged as a promising data preprocessing technique, has been widely used in forecasting fields to solve complex nonstationary, nonlinear time series forecasting issues. Numerous studies have shown that the data decomposition algorithm has a significant influence on forecasting performance. Thus, a suitable decomposition method should be selected to identify and extract the inner characteristics of the original time series. In previous studies, such as Lin et al. (2017) and Yang et al. (2019a), VMD, proposed by Dragomiretskiy and Zosso (2014), has been shown to be superior to other algorithms, including EMD, EEMD, and CEEMD. Therefore, VMD is introduced into the data preprocessing module to mitigate the negative noise information influence.
The main procedure of VMD is as follows.
Step 1: setting parameters of VMD, while each mode y 1 k , center pulsation w 1 k , and Lagrangian multipliers λ are initialized.
Step 2: y k and w are updated by where n is the iterations number, f(t) and y k are the original time series and the kth component, whereas f(w), y i (w), λ(w), and y n+1 k (w) denote the Fourier transforms of f(t), y i (t), λ(t), and y n+1 k (t), respectively.
Step 3: λ can be updated by Step 4: if k y n+1 k − y n k 2 2 / y n k 2 2 < e, the VMD algorithm is stopped, and a series of band-limited modes is returned; otherwise, return to Step 2 to repeat the iteration.

Outlier-Robust Forecasting Module
The basic forecasting model is the important foundation of a hybrid model, which can make a significant difference in forecasting results. If outliers are present within the dataset, the performance of the model developed will be significantly affected. Considering the significance of outlier modeling and outlier robustness, the ORELM model is acknowledged as a potential contributor for modeling data with outliers. Therefore, the ORELM model is introduced into daily PM 2.5 concentration forecasting to design an outlier-robust forecasting module. The original version of the ORELM model is ELM, developed by Huang et al. (2004), which has many merits, such as its simple structure, better performance, fast computation speed, and the fact that it does not need a large number of samples. Furthermore, previous studies have revealed that ELM methods are superior to some typical ANN methods in solving forecasting issues (Yang et al., 2019a), and it has become one of the most promising approaches.
Given a training dataset with M samples, i.e., (x t , y t ), t 1, /, M, the ELM model for input data x t and output data y t can be presented as where L denotes the number of hidden layer nodes, w i and b i denote the input weight and hidden bias, G is the excitation function, β i represents the connected weight between the ith hidden layer node and the output layer, and y represents the forecasting results. By defining the hidden layer output matrix, i.e., H, The ELM model presented in Eq. 8 can be rewritten as The optimal solution of β can be obtained by solving min where H † represents the Moore-Penrose generalized inverse matrix of H; the corresponding formula is As mentioned above, to enhance the ELM model's robustness when modeling data with outliers, the ORELM model is developed by Zhang and Luo (2015). The core idea is redefining the minimum problem as where e represents training error and k is the regularization parameter.
To solve the newly defined problem, the augmented Lagrange multiplier (ALM) algorithm is adopted, and the corresponding iteration process is defined as where λ represents the Lagrange multiplier vector, μ is the penalty parameter, and β t+1 and e t+1 are defined as Frontiers in Environmental Science | www.frontiersin.org October 2021 | Volume 9 | Article 747101

Nonlinear Correction Module
For a forecasting model, forecasting error is inevitable, but the short-term trend in the variation of the forecasting error can be anticipated by establishing a nonlinear model (Vukicevic, 1991).
To further improve the performance of the system developed, a nonlinear correction module based on an error ensemble strategy is developed to mine information in the forecasting results, which is composed of three steps.
Step 1: generating the error sequence Defining the actual value of the t th datum as A(t) and forecasting the value of the t th datum as F(t), the forecasting error value of the t th datum can be obtained by Step 2: developing the error forecasting model Defining the error value of the t−d th datum as E(t−d), according to the detailed error sequence, the forecasting model can be developed and denoted as f. The error forecasting value of the t th datum, named EF(t), can be obtained by Step 3: correcting the forecasting results To obtain a final result on the basis of the original forecasting results and corresponding error forecasting results, an error ensemble strategy based on ORELM is proposed, which fully exploits the advantages of ORELM and is equipped with outlier robustness. By developing an outlier-robust ensemble model, denoted as En.f, the final forecasting results of the t th datum, i.e., FF(t), are

Artificial Intelligence-Based Distribution Evaluation
Distribution evaluation plays a vital role in many fields, such as wind energy evaluation, time series analysis, and interval forecasting. In recent years, in order to further mine data characteristics, researchers have focused on applying different distribution functions to fit the experimental data and obtain a suitable distribution; then, the interval forecasting results can be obtained according to the interval forecasting theory and point forecasting results. However, to the best of our knowledge, the related research is well validated in many fields, but so far, few studies have involved research on or application to daily PM 2.5 concentration forecasting. In this context, four typical distributions, i.e., Weibull, Gamma, Rayleigh, and Lognormal, are introduced in this study to fit the daily PM 2.5 concentration data. In general, the goodness of fit (0 ≤ R 2 ≤ 1) was employed to measure the fitting performance of one distribution. Traditionally, the maximum likelihood estimation (MLE) method is used to estimate the distribution function's parameters. However, the MLE method may not obtain the optimal distribution parameters. To the best of our knowledge, the larger the R 2 value, the more optimal the distribution. As a result, the optimal distribution determination problem can be converted into solving the maximum value problem. Inspired by Wang et al. (2020b), Schwarz et al. (2020), and Ließ et al. (2021), artificial intelligence optimization can be considered a promising technique for searching for the optimal distribution parameters. Based on this idea and considering the limitations of the traditional method, the artificial intelligence-based distribution evaluation is proposed to obtain the optimal distribution in this study. In order to obtain the optimal distribution, specifically, an advanced optimization algorithm named grey wolf optimizer (GWO) is adopted to search for the optimal parameters of specific distribution by maximizing the values of R 2 . In this study, the minus R 2 is defined as the objective function of GWO-based distribution evaluation. Finally, the distribution with the best R 2 value among all distributions is selected as the optimal distribution of PM 2.5 , which can be combined with interval forecasting theory to achieve interval forecasting.

Interval Forecasting Theory
Given the significance level α, actual value A t , and lower and upper limits (L, U), the probability can be given by For a random variable time series, Eq. 19 can be rewritten as Supposing that the forecasting value has a similar distribution function, the estimated variance can be determined, and then the following conditional probability formula can be obtained as The lower and upper limits can be obtained by

CONSTRUCTION OF THE OUTLIER-ROBUST FORECASTING SYSTEM
The outlier-robust point and interval forecasting system is constructed in this section; the details of the system design and system evaluation are as follows.

System Design
The system design is composed of the point forecasting part and the interval forecasting part, which can provide deterministic information and uncertainty information in the future, respectively.

Point Forecasting
The point forecasting part is developed based on three modules, i.e., a data preprocessing module, an outlier-robust forecasting module, and a nonlinear correction module, and consists of three phases as follows.
◆ Phase I: data preprocessing. Considering the negative influence of outlier data and noise information in original daily PM 2.5 concentration data, a data preprocessing module is designed based on outlier handling and data decomposition that can effectively solve the outlier data and noise information issues simultaneously. Specifically, on the one hand, the HF algorithm is employed to detect and correct the outliers in the original data, which can eliminate the outliers' negative influence on the model's development from the perspective of improving data quality. On the other hand, the advanced data decomposition method named VMD is performed to decompose the processed data into a number of modes, which can eliminate the noise's negative effect on the model's performance from the perspective of signal denoising and helps the model effectively capture the main features of the daily PM 2.5 concentration data. ◆ Phase II: forecasting preprocessed data. In this phase, some modes obtained in Phase I are forecasted, and then, the forecasting results can be obtained. As a result, it is necessary to select a suitable model for preprocessed data forecasting. Most previous studies tended to develop an optimized ANN model for mode forecasting while ignoring the instability parameter setting problem of artificial intelligence optimization algorithms, which not only cannot guarantee the forecasting model's stability but also bring time-consuming issues. In this context, the ORELM model, with its simple structure, high forecasting performance, and fast calculation speed, better forecasting ability than some typical ANN models, and better outlier robustness than the original ELM model, is selected to design the outlier-robust forecasting module to obtain the future value of each mode, and finally, by summing the forecasting results of each mode, the results of daily ahead PM 2.5 concentration forecasting can be achieved. ◆ Phase III: forecasting results correction. Most previous studies have focused on using data decomposition and optimization to enhance forecasting ability while ignoring the significance of forecasting results correction. Therefore, in this phase, a nonlinear correction module based on an error ensemble strategy is designed to mine information in the forecasting results and thereby enhance the model's performance, which is composed of three steps, i.e., generating an error sequence, developing an error forecasting model, and correcting the forecasting results. Following these three steps, the results of the point forecasting can be obtained; these results can provide deterministic information in the future.

Interval Forecasting
The interval forecasting part is developed based on the previous point forecasting results, artificial intelligence-based distribution evaluation, and interval forecasting theory and consists of two phases as follows.
◆ Phase I: distribution evaluation. Distribution evaluation is the crucial basis of the interval forecasting part, whereas the traditional MLE method may not fit the optimal distribution for a specific PM 2.5 concentration dataset. In order to solve this issue, an artificial intelligence-based distribution evaluation is proposed to obtain the optimal distribution. Specifically, the advanced optimization algorithm named GWO is selected to search the parameters of four typical distributions, i.e., Weibull, Gamma, Rayleigh, and Lognormal. Finally, the distribution function with the best R 2 is considered the optimal distribution. ◆ Phase II: obtaining interval forecasting results. According to the point forecasting results, the interval forecasting results are estimated using the optimal distribution determined in Phase I and interval forecasting theory, which can provide uncertainty information in the future.

System Evaluation
This section is designed to provide system evaluation metrics, including point forecasting evaluation and interval forecasting evaluation.

Point Forecasting Evaluation
To evaluate the forecasting performance of the system developed for daily PM 2.5 concentration point forecasting, the three typical metrics listed in Table 1 are selected in this study. In addition to evaluating the forecasting accuracy using these three typical evaluation metrics, to measure the similarity of the forecasting value curve and actual value curve, grey relational analysis (GRA)  is introduced into point forecasting evaluation, which provides a new metric named grey relational degree (GRD) for further analysis. Moreover, forecasting stability is another important metric in practical application. Thus, a metric named variance ratio (VR)  is introduced into the point forecasting evaluation, which can measure the forecasting stability of different models by considering the variance between the actual and the forecast values.

Interval Forecasting Evaluation
To evaluate the interval forecasting performance, two widely used metrics named forecasting interval coverage probability (FICP) and forecasting interval normalized average width (FINAW) are introduced into the interval forecasting evaluation. Specifically, the FICP indicator is selected to measure the probability that the actual observation value falls into the forecasting interval, which can reflect the reliability of the interval forecasting results. The larger the FICP value, the higher the interval forecasting accuracy. Furthermore, to avoid the excessive pursuit of the reliability of interval forecasting while ignoring the effective measure of interval forecasting width for interval forecasting results, FINAW is employed to express the width of interval forecasting results. The FICP and FINAW can be calculated by where c i 1, A i ∈ L i ,U i 0, A i ∉ L i ,U i and N is the length of the time series, where U i and L i denote the forecasting interval's upper and lower limits, respectively, at time i, and R is the range of actual values.

EXPERIMENTAL ANALYSIS
The experiments and analysis are presented in this section, which is mainly composed of data description, experiment design, and four detailed experiments.

Data Description
To validate the ability of the outlier-robust system developed to perform point forecasting and interval forecasting of daily PM 2.5 pollution, Jinan and Zhengzhou are considered as the study areas; two datasets collected from these two study areas are used as illustrative empirical studies in this study. Jinan, the capital city of Shandong Province, is located in the middle of China. Zhengzhou, the capital city of Henan Province, is located in the middle part of the Yellow River. Specifically, two daily PM 2.5 concentration datasets, covering 1 yr from July 1, 2017, to June 30, 2018, are employed in this study. In the experiment, the data, from July 1, 2017, to May 31, 2018, are employed as training data for the development of the proposed system, whereas the data from June 1, 2018, to June 30, 2018, are considered as testing data to test the forecasting performance of the system developed.

Experiment Design
As mentioned above, in this study, an outlier-robust point and interval forecasting system is developed, which is composed of a data preprocessing module, an outlier-robust forecasting module, a nonlinear correction module, artificial intelligence-based distribution evaluation, and interval forecasting theory to obtain future deterministic information and uncertainty information about daily PM 2.5 pollution. To verify the forecasting superiority of the system developed, sufficient empirical research should be carried out. In addition to comparing the performance of the system developed with that of the other types of forecasting models, the contribution of each component proposed or employed in the system developed should also be proved by designing appropriate comparative studies. For this purpose, this study designs four experiments to conduct a convincing evaluation of the system developed. Specifically, in Experiment I, the effectiveness of outlier handling and modeling in the system developed is verified from the perspectives of data preprocessing and model selection. In Experiment II, the effectiveness of data decomposition in the system proposed is compared with other decomposition algorithm-based models and a model without a decomposition preprocess. In Experiment III, a nonlinear correction module is developed to correct the forecasting results, which is designed to compare the proposed system with the model without correcting the process and the model with a simple error-addition strategy.
It should be noted that the experiments for each model in Experiments I-III are carried out 100 times in this study, and the average values of the forecasting results are considered the final forecasting results for practical application and model comparison, which can ensure that the system developed is more reliable, accurate, and independent of random factors to some extent. In Experiment IV, different distributions of daily PM 2.5 concentration are compared to obtain the optimal distribution, and the interval forecasting results based on point forecasting are obtained and evaluated by two typical metrics.

Experiment I: The Effectiveness of Outlier Handling and Modeling
To evaluate the effectiveness of outlier handling and modeling, eight models, i.e., ELM, regularized ELM (RELM), weighted RELM (WRELM), ORELM, HF-ELM, HF-RELM, HF-WRELM, and HF-ORELM, are proposed and tested. The MAE, RMSE, MAPE, VR, and GRD values of these eight models are shown in Table 2. Meanwhile, the results of the different models in the two cities are depicted in Figures 1, 2, which indicate that the ORELM model is superior to ELM, RELM, and WRELM, whereas the HF-ORELM model is superior to the seven other models. As shown in  Summary: by taking Zhengzhou as an example, the improvement percentage values of MAPE between the different models are employed to summarize the contribution and effectiveness of outlier handling and modeling in this study. The detailed results are 8.2222% (HF-ELM vs ELM), 7.4577% (HF-RELM vs RELM), 10.6916% (HF-WRELM vs WRELM), 8.8642% (HF-ORELM vs ORELM), 14.3166% (ORELM vs ELM), 13.7965% (ORELM vs RELM), and 12.0041% (ORELM vs WRELM). It can be concluded that the HF algorithm and ORELM model are suitable for outlier handling and modeling, which make a great contribution to the success of the system developed in this study.

Experiment II: The Effectiveness of Data Decomposition
To verify the contribution of data decomposition in the proposed data preprocessing module and the superiority of the forecasting results of the outlier-robust forecasting module developed, four models, i.e., HF-ORELM, HF-EMD-ORELM-S, HF-EEMD-ORELM-S, and HF-VMD-ORELM-S, are developed and compared in Jinan and Zhengzhou. In detail, the HF-EMD-ORELM-S, HF-EEMD-ORELM-S, and HF-VMD-ORELM-S employ different data decomposition algorithms to decompose the data after outlier handling into some modes, and the simple addition way is used to add all modes' forecasting results to obtain the daily PM 2.5 concentration forecasting results. The MAE, RMSE, MAPE, VR, and GRD values of HF-ORELM, HF-EMD-ORELM-S, HF-EEMD-ORELM-S, and HF-VMD-ORELM-S are shown in Table 3. Moreover, the forecasting results of these four models in the two cities are shown in Figure 3, which indicates that the HF-VMD-ORELM-S model is superior to the original HF-ORELM model and the EMD-or EEMD-based HF-ORELM model. In this experiment, two comparisons can be designed as follows: 1) Comparison I is proposed to validate the superiority of the data decomposition algorithm in the system developed by comparing the HF-VMD-ORELM-S with other decomposition method-based forecasting models, i.e., HF-EMD-ORELM-S and HF-EEMD-ORELM-S. It can be observed that the HF-EMD-ORELM-S model obtains worse forecasting performance compared with the EEMD-and VMD-based models, whereas the VMD-based model achieves better forecasting performance compared with the EMD-and EEMD-based models. For example, for daily PM 2.5 concentration forecasting in Zhengzhou, the MAE, RMSE, MAPE, VR, and GRD values of HF-VMD-ORELM-S are 1.1259, 1.5228, 3.8169%, 0.9523, and 0.9222, respectively, whereas the metric values of HF-EMD-ORELM-S are 4.2089, 5.9923, 14.1646%, 0.9024, and 0.7705, and the values of HF-EEMD-ORELM-S are 2.1140, 2.7875, 7.0637%, 0.8161, and 0.8652. It is obvious that there are significant differences in the forecasting power of these three models, which further demonstrates the significance of selecting a suitable data decomposition algorithm for the data preprocessing module and the system developed. Therefore, in this study, the VMD algorithm is combined with outlier handling to design the data preprocessing module, which also makes great contributions to the success of the system developed. 2) Comparison II is designed to verify the superiority of the outlier-robust forecasting module developed in daily PM 2.5 concentration forecasting. In Experiment I, the superiority of HF-ORELM over ORELM has been proven, which means the outlier handling is an indispensable part of the forecasting system developed. Against this background, there is no need to compare the proposed outlier-robust forecasting module's results with those of the individual ORELM model. As a result, It can be concluded that the VMD algorithm is superior to the EMD and EEMD algorithms and is a promising technique for daily PM 2.5 concentration decomposition, which can also make a great contribution to the success of the system developed.

Experiment III: The Effectiveness of the Nonlinear Correction Module
As mentioned above, the third module, named the nonlinear correction module, is proposed to correct the results of the outlier-robust forecasting module to further improve the daily PM 2.5 concentration forecasting performance. To prove the superiority and effectiveness of the proposed nonlinear correction module and the system developed for point forecasting, the performance of the point forecasting part developed, i.e., HF-VMD-ORELM+EnError, is compared with HF-VMD-ORELM-S and HF-VMD-ORELM+Error in this section. In detail, the HF-VMD-ORELM-S model without a correcting process is the best model in Experiment II, which can provide the results of the devised forecasting module, whereas the HF-VMD-ORELM+Error model is a model with a simple error-addition strategy. The MAE, RMSE, MAPE, VR, and GRD values of the system developed, HF-VMD-ORELM-S, and HF-VMD-ORELM+Error are listed in Table 4; meanwhile, the results of these three models are shown in Figure 4. Based on Experiment III, the following conclusions can be obtained: 1) The HF-VMD-ORELM+Error model performs better than HF-VMD-ORELM-S model in Jinan but performs worse than HF-VMD-ORELM+Error in Zhengzhou, which indicates that the simple error-addition strategy cannot guarantee the effectiveness of error correction. Therefore, how to correct the forecasting results is a challenging issue in forecasting fields. In other words, the method of correcting forecasting results plays a vital role in the success of the system developed.
In this context, a nonlinear correction module based on an error ensemble strategy is presented to further improve the model's forecasting performance. 2) By comparing the HF-VMD-ORELM+EnError model and the HF-VMD-ORELM+Error model, it can be found that the main difference between these two models is the forecasting results correcting method. As shown in Table 4, the forecasting results obtained  In previous experiments, the superiority of the HF-VMD-ORELM-S model over other benchmark models has been proven. Therefore, considering the system developed performs better than HF-VMD-ORELM-S, we can conclude that the system designed in this study is superior to other models and can be widely employed in daily PM 2.5 concentration forecasting.
Summary: by taking the MAPE metric as an example, the improvement percentage values between different models are employed to summarize the contribution and effectiveness of  It can be concluded that the proposed nonlinear correction module is not only effective for improving the final forecasting results but also is superior to the HF-VMD-ORELM+Error model. Furthermore, the HF-VMD-ORELM+Error model may perform worse than the HF-VMD-ORELM-S model. In other words, the proposed nonlinear correction module is suitable for correcting forecasting results, which can contribute to improving the performance of the system developed.

Experiment IV: Interval Forecasting
In the system developed, the interval forecasting can be achieved by the proposed interval forecasting part according to the results of the point forecasting part, artificial intelligence-based distribution evaluation, and interval forecasting theory. In order to perform the interval forecasting, a distribution evaluation of daily PM 2.5 concentration data is conducted. As mentioned above, the traditional MLE method may not obtain the optimal distribution function for a specific PM 2.5 concentration dataset, whereas artificial intelligence optimization is a powerful technique for determining the optimal distribution. Therefore, in this study, the advanced optimization algorithm named GWO is selected to fit four typical distributions, i.e., Weibull, Gamma, Rayleigh, and Lognormal. In order to prove the superiority of GWO, detailed distribution is also determined by MLE, and the parameters and R 2 values provided by MLE and GWO are presented in Table 5. Furthermore, the comparison is also depicted in Figure 5. On the basis of Table 5 and Figure 5, we find that the GWObased distribution evaluation can obtain the best R 2 values for each distribution function, which indicates that the GWO-based distribution evaluation is superior to the MLE method and is suitable for fitting the detailed distribution. As a result, the results of the GWO-based distribution evaluation can be compared by R 2 . As shown in Table 5, the R 2 values of Weibull, Gamma,Rayleigh,and Lognormal are (0.9879,0.9936,0.9873,and 0.9976) and (0.9779,0.9818,0.9666,and 0.9879) in Jinan and Zhengzhou, respectively. It can be observed that the Lognormal distribution achieves the largest R 2 , which means that the Lognormal distribution can effectively fit the daily PM 2.5 concentration data in Jinan and Zhengzhou. Thus, the optimal Lognormal distribution obtained can be combined with the point forecasting results and interval forecasting theory to achieve the final interval forecasting.
The interval forecasting results under different significance levels are depicted in Figure 6, and the corresponding evaluation metric values are listed in Table 6. From Table 6, we can find that the interval forecasting performances for Jinan and Zhengzhou are different at the same significance level. For example, when α 0.30, the FINAW and FICP values for Jinan and Zhengzhou are (0.6268, 100.0000%) and (0.4695, 96.6667%), respectively. The main reasons for this phenomenon are that the interval forecasting performance largely depends on the point forecasting performance. As the system developed has achieved excellent point forecasting performance, it has also achieved ideal interval forecasting results. Moreover, for the same dataset, the FINAW and FICP values under five significance levels are different. For instance, for the Zhengzhou dataset, the FINAW and FICP values for α 0.20 and α 0.25 are (0.7711, 100.0000%) and (0.6097, 96.6667%), respectively. Furthermore, the solid lines represent the actual values, and the shaded areas represent the forecasting intervals in Figure 6. Obviously, as most of the observations fall into the shaded area, the interval forecasting ability of the system established can be considered effective and good. According to the abovementioned analysis and discussion, we can reasonably conclude that the system developed can be a promising tool for daily PM 2.5 concentration interval forecasting.

CONCLUSION
Forecasting air pollution is not only a challenging research topic but also a growing concern. To model and forecast the complex PM 2.5 concentration time series, in this study, a novel outlierrobust point and interval forecasting system is developed, which attempts to mitigate or solve some of the challenges in current studies. In the system developed, the point forecasting part is designed to provide future deterministic information on daily PM 2.5 concentration, whereas the interval forecasting part is devised to present future uncertainty information. More specifically, three modules, named the data preprocessing module, the outlier-robust forecasting module, and the nonlinear correction module, are proposed in the point forecasting part. The data preprocessing module, considering the negative influence of outliers and noise on the development of the model, is designed on the basis of outlier handling and data decomposition for the purpose of solving the outlier data and noise information issues simultaneously. Moreover, in order to obtain a forecasting model with outlier robustness, the ORELM model with superior performance is selected to design the outlierrobust forecasting module to forecast each mode. Furthermore, the nonlinear correction module is developed based on an error ensemble strategy, which can mine information in the forecasting results and further improve the model's forecasting performance. Afterward, the interval forecasting part is developed based on artificial intelligence-based distribution evaluation and interval forecasting theory, which can be incorporated with the point forecasting results to obtain the range of future changes. The experimental results illustrate that the system developed can not only perform better than other compared models in point forecasting but also provide uncertainty information in the future. Moreover, outside the field of daily PM 2.5 concentration forecasting, the system developed can also be employed to solve other challenging issues, including energy forecasting, economic forecasting, and financial forecasting.
Although the system developed shows better performance in daily PM 2.5 concentration forecasting, there are still some limitations that must be considered in future research. For example, other features neglected in this study may make the forecasting more reliable and practical and can further improve the forecasting performance, which can be a promising subject for future studies. Moreover, only 1-day-ahead forecasting is conducted by the system developed; how to achieve multi-day-ahead forecasting is of great importance and worth studying but is still a challenging task. As a result, further studies about multi-day-ahead point and interval forecasting for PM 2.5 concentration can be considered an important research direction in future studies.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: China Air Quality Online Monitoring and Analysis platform (https://www.aqistudy.cn).

AUTHOR CONTRIBUTIONS
ZY and XF conducted the research.


extfieldsofstudy:
Computer Science
Engineering

text:
RACH in Self-Powered NB-IoT Networks: Energy Availability and Performance Evaluation

NarrowBand-Internet of Things (NB-IoT) is a new 3GPP radio access technology designed to provide better coverage for a massive number of low-throughput low-cost devices in delay-tolerant applications with low power consumption. To provide reliable connections with extended coverage, a repetition transmission scheme is introduced to NB-IoT during both Random Access CHannel (RACH) procedure and data transmission procedure. To avoid the difficulty in replacing the battery for IoT devices, the energy harvesting is considered as a promising solution to support energy sustainability in the NB-IoT network. In this work, we analyze RACH success probability in a self-powered NB-IoT network taking into account the repeated preamble transmissions and collisions, where each IoT device with data is active when its battery energy is sufficient to support the transmission. We model the temporal dynamics of the energy level as a birth-death process, derive the energy availability of each IoT device, and examine its dependence on the energy storage capacity and the repetition value. We show that in certain scenarios, the energy availability remains unchanged despite randomness in the energy harvesting. We also derive the exact expression for the RACH success probability of a {randomly chosen} IoT device under the derived energy availability, which is validated under different repetition values via simulations. We show that the repetition scheme can efficiently improve the RACH success probability in a light traffic scenario, but only slightly improves that performance with very inefficient channel resource utilization in a heavy traffic scenario.

I. INTRODUCTION
The Internet of Things (IoT) is a novel paradigm that is rapidly gaining interest in modern wireless telecommunications to support connections of billions of miscellaneous innovative devices. Third Generation Partnership Project (3GPP)

A. NarrowBand-Internet of Things
NB-IoT is a new 3GPP radio-access technology developed from existing LTE functionalities, whereas some features of its specification deemed unnecessary for LPWA IoT needs have been stripped out [9]. Because of this, NB-IoT can provide unique advantages for various IoT services over other technologies like 2G, 3G or LTE. LPWA networks mainly require deep/wide coverage, low power consumption, massive connections, and lower cost. The inherent characteristics of NB-IoT make it a good fit for LPWA deployment as shown in Fig. 1.  Extending battery lifetime is one of the features of NB-IoT. There are several ways to reduce power consumption and achieve lifespans of IoT devices for more than 10 years. In the most simple way, we can switch the device to sleep mode when it does not work, so it doesn't waste power while waiting. Power management is fundamentally the tradeoff between message frequency, device sleep cycles, and business case needs. Deep coverage is another feature of NB-IoT, where it is designed to improve indoor coverage by 20 dB compared to conventional GSM/GPRS. This is achieved by a higher power density, as radio transmissions are concentrated on a narrower carrier bandwidth of just 180 kHz. The Coverage Enhancement (CE) feature additionally offers the capability to repeat the transmission of a message when there exist poor coverage conditions but at the expense of a lower data rate.
NB-IoT reuses the LTE design extensively, such that the time required to develop full specifications and products is significantly reduced. In the downlink of NB-IoT, OFDMA (Orthogonal Frequency Division Multiplexing) technology is adopted with the sub-carrier spacing of 15 kHz and 12 subcarriers make up the 180 kHz channel [10]. In the uplink, SC-FDMA (Single Carrier Frequency Division Multiple Access) technology including single sub-carrier and multiple subcarrier are adopted. A single sub-carrier technology with 3.75 kHz and 15kHz is adopted with carrier spanning over 48 and 12 respectively. In this paper, we focus on NB-IoT Physical Random Access CHannel (NPRACH) with single sub-carrier spacing of 3.75 kHz.
NPRACH resources consist of the assignment of time and frequency resources and occur periodically. In the NPRACH, a random access preamble is transmitted which is the first step of random access procedure that enables the device to establish a data connection with its associated Evolved Node B (eNB) [10]. To improve the quality of service and reduce the power consumption of IoT devices, efficient RACH procedures need to be proposed and analyzed [9]. In [11]- [13], mathematical models of contention-based RACH focusing on the Signalto-Interference-plus-Noise Ratio (SINR) outage or collision have been studied. However, to the best of our knowledge, most works have focused either on studying the SINR outage without considering collision or studying the collision from the single cell point of view and most results are based on the uplink power control due to its analysis simplicity.
Our previous work [14] has provided the preamble transmission model without considering the collision, [15] has considered the collision model without considering the repetition scheme, and [16] has considered the repetition scheme based on the uplink power control for simplicity. However, the IoT device adopts the cell specific maximum transmission power without power control when transmitting more than two preamble repetitions, which is motivated by the statement in 3GPP standard [5]. In this scenario, it is unknown 1) to what extent the transmission power affects the RACH success; 2) to what extent the repetition transmission scheme affects the RACH success; 3) how to choose the repetition value in different traffic scenarios to balance the RACH success probability and data transmission channel resources. To solve these problems, we present a novel mathematical framework to analyze and evaluate the RACH success probability, taking into account the SINR outage events as well as the collision events at the eNB, where the IoT devices adopt fixed transmission power.
Generally speaking, in the NB-IoT network, the physical layer parameters and network topology can strongly affect the RACH performance, due to that the received SINR at the eNB can be severely degraded by the mutual interference generated from massive IoT devices. In this scenario, the random positions of the transmitters make accurate modeling and analysis of this interference even more complicated. It is worth noting that stochastic geometry has been regarded as a powerful tool to model and analyze mutual interference between transceivers in the wireless networks [17], such as conventional cellular networks [18], wireless sensor networks [19], cognitive radio networks [20], and heterogenous cellular networks [21].
However, conventional stochastic geometry works [17]- [21] focused on analyzing normal uplink or downlink scheduled data transmission channel, where the intra-cell interference is not considered, due to the ideal assumption that each orthogonal sub-channel is not reused in a cell. The work [22] still focused on the scheduled data transmission, though it has considered the intra-cell interference over non-orthogonal subchannels. All these stochastic geometry works are different from the RACH analysis in this paper, where massive IoT devices in a cell may randomly choose and transmit the same preamble to their associated eNB to request for channel resources for data transmission, and different preambles represent orthogonal sub-channels. Thus only IoT devices choosing the same preamble (i.e., sub-channel) have correlations. Note that IoT devices associated to the same eNB may choose the same preamble and collisions occur when two or more IoT devices transmit the same preamble simultaneously, such that the intra-cell interference and collisions are considered.
In practical NB-IoT networks, mutual interference among transmissions is much more intricate than the conventional LTE network systems. In this paper, we develop a novel mathematical framework for NB-IoT networks using stochastic geometry for RACH analysis with some challenges: 1) the distribution of the random transmission distances is considered, due to that the IoT device transmits with fixed transmission power; 2) the intra-cell interference is considered, due to that the massive IoT devices in a cell may randomly choose and transmit the same preamble using the same sub-channel; 3) the temporally correlated mutual interference is considered, due to that each IoT device remains spatially static during each repetition; 4) both the SINR outage and the collision events are considered, due to that collisions occur when two or more IoT devices transmit the same preamble simultaneously; 5) the network point of view analysis is considered instead of the single cell point of view, which is difficult in capturing both interference and collision generated from IoT devices transmitting with different transmit distances, due to the many concurrent transmissions and the interference experienced in each cell is different.

B. Energy Harvesting
Human-operated cellular devices, such as smart phones, can be charged at will. But IoT devices are often located at remote and hard to reach locations, such as underground or in tunnels, without access power supply, which may be inconvenient, dangerous, expensive, or even impossible to change the battery. Hence, the battery energy storage highly determines the lifetime of the whole device, which provides a strong motivation for powering IoT devices by harvesting energy, such that networks consisting of energy harvesting or rechargeable batteries can survive perpetually [23]. Practically, energy can be harvested from renewable environmental sources including thermal, solar, wind, etc. [24] and radio signals of different frequencies such as radio broadcasting [25]. In these cases, the network performances are often tied closely to the efficiency in energy harvesting and utilization.
It is difficult to predict the time and the amount when the energy is available, as the energy arrival process is also random and dynamic. Energy buffer, i.e., battery storage, which collects harvested energy for signal processing and communication, is introduced in order to mitigate the unpredictability of the energy. Moreover, energy harvesting rates achievable today still fall short of typical power consumption levels. The harvested energy need to be accumulated in storage modules (e.g., capacitors or batteries) to a sufficient energy level to operate the IoT device. Nevertheless, the energy buffer is limited. It remains challenging to assign the amount of energy in terms of uncertain energy sources. To have an indepth study on the harvesting and utilization of energy, we need an analytically simple yet practically accurate model of the harvested energy. Several models on the harvested energy have been used in literature, assuming a time-slotted system. The arriving of harvested energy is known as a deterministic way in [26] while unknown in [27]. In [28] [29], the stationary Markovian models of the harvested energy were studied, which are analytically simple and are thus useful to provide insights for solving some key theoretical problems. However, the validity of these models has not been formally justified with empirical measurements, and hence it is not known if the insights are useful in practice. In [30], a more general analytical model for the harvested energy was provided with support from empirical measurements. In their empirical measurements, the harvested energy may be a Markovian process. In [31], sleep/wake-up strategies for various factors were studied, which determines the optimal parameters of the solar energy harvest based strategy using a bargaining game model. However, to the best of our knowledge, there has been no work studying energy harvesting in NB-IoT networks. In this work, we consider NB-IoT networks with the IoT devices harvesting energy from nature, e.g., solar cells, microbial fuel cells, and water mills, etc.

C. Contributions and Outcomes
The contributions of this paper can be summarized as follows: 1) Using stochastic geometry, we present a tractable analytical framework for the self-powered NB-IoT network via energy harvesting from natural resources, in which the locations of the eNBs and the IoT devices are modeled as two independent Poisson Point Processes (PPPs) in the spatial domain.
2) We model the arrival of harvesting energy as independent Poisson arrival processes. Using tools from Markov stochastic process, we characterize the fraction of time each IoT device kept in ON state as the energy availability. We first model the temporal dynamics of the energy level as a birth-death process, and then derive the expression of energy availability of each IoT device using hitting time analysis.
3) Based on the derived energy availability of the IoT device, we drive the expression for the RACH transmission success probability of a randomly chosen IoT device with fixed transmission power under both SINR outage and collision conditions in the NB-IoT network. Furthermore, we develop a realistic simulation framework to capture the randomness locations, the preamble transmission as well as the RACH collision, and verify our derived RACH success probability. 4) Our results show that the energy availability of the IoT device increases at first and then remain unchanged, the RACH success probability increases, and the repetition efficiency decreases when increasing the repetition value. As such, the repetition value needs to be optimized. It is also noticed that in certain scenarios, the energy availability remains unchanged despite randomness in the energy harvesting. The results also show that there is an upper limit on transmission power and too large transmission power will waste energy.
The rest of the paper is organized as follows. Section II presents the network model. Section III derives energy availability of IoT device and analysis the actual result conditioning on some specific energy utilization strategies. Section IV derives the RACH transmission success probabilities of a randomly chosen IoT device. Section V presents the simulation framework. Finally, Section VI concludes the paper.

II. SYSTEM MODEL A. Network Description
We consider an uplink stochastic geometry model for NB-IoT system, consisting of a single class of eNBs and IoT devices, that are spatially distributed in the Euclidean plane R 2 following two independent homogeneous Poisson Point Processes (PPPs) Φ B and Φ D with intensities λ B and λ D , respectively. Without loss of generality, we assume each IoT device remains spatially static during time slots [16]. Same as [18], we assume that each IoT device associates with its geographically nearest eNB, where a Voronoi tesselation is formed. A standard power-law path-loss model is considered, where the path-loss attenuation is defined as r −α , with the propagation distance r and the path-loss exponent α. In addition, we consider a Rayleigh fading channel, where the channel power gain h is assumed to be exponentially distributed random variable with unit mean, i.e., h ∼ Exp(1). All channel gains are assumed to be independent and identically distributed (i.i.d.) in space and time.

B. Random Access Procedure
In the uplink of NB-IoT, data can only be transmitted via the dedicated uplink data transmission channel, Narrowband Physical Uplink Shared CHannels (NPUSCH), which is scheduled by the associated eNB. Before resource scheduling, the IoT device needs to execute a RACH to request uplink channel resources with the associated eNB. The RACH procedure of NB-IoT has the simplified message flow as for LTE, however, with different parameters [10]. There are two types of RACH procedures: contention-free and contention-based. The former is used to perform handover, whereas the latter is used otherwise [32]. The four-step contention-based RACH used by NB-IoT is described as follows.
In step 1, the device first transmits a preamble (msg1) on the NPRACH during the first Random Access Opportunity (RAO). The eNB periodically informs the devices about a set of up to 64 orthogonal preamble sequences from which the IoT device can make a choice. Collisions occur when two or more devices transmit the same preamble sequence simultaneously. In step 2, the IoT device sets a Random Access Response (RAR) window and waits for the eNB to transmit a RAR (msg2) with an uplink grant for the transmission of a message in the following step. If the device that sent a preamble sequence does not receive a msg2 from the eNB within a certain period of time, it enters a backoff period, trying to access the network once this period has expired. In step 3, the IoT device that successfully receives its RAR transmits a Radio Resource Control (RRC) Connection Request (msg3) with identity information to eNB. If two or more devices have chosen the same preamble sequence in a RAO, they will receive the same grant in the RAR message, and thus, their msg3 transmissions will collide. In step 4, the eNB transmits a RRC Connection Setup (msg4) to the IoT device when it successfully receives msg3. More details on the RACH can be found in [8] [14].
Massive connections in NB-IoT make the simultaneous RACH requests under a limited number of available preambles one of the main challenges, thus we focus on the contention of preamble in step 1 of contention-based RACH, with the assumption that steps 2, 3, and 4 of RACH are always successful whenever step 1 is successful. If step 1 in RACH fails, i.e. the associated RAR message was not received, the IoT device needs to transmit another preamble in the next available RAO. That is to say, a RACH procedure is always successful if the IoT device successfully transmits the preamble to its associated eNB. In this case, the failure of this preamble transmission can result from the following two reasons: 1) the eNB cannot decode the preamble due to the low received SINR; 2) the eNB successfully decoded the same preamble from two or more IoT devices in the same time and causes the collision.
It is known that the collision event in step 1 of RACH can be detected by the eNB, when the collided IoT devices are separable in terms of the power delay profile [8]. Our model follows the assumption of collision handling in [15], where collision events are detected by eNB after it decodes the preambles in step 1 of RACH, and then no response will be fed back from the eNB to the IoT devices, such that it can not proceed to the next step of RACH [33].

C. Narrowband Physical Random Access CHannel
NB-IoT technology occupies a frequency band of 180 kHz bandwidth [34], which corresponds to one resource block in LTE transmission. In the NPRACH, a preamble is transmitted based on symbol groups on a single subcarrier. A preamble consists of four symbol groups transmitted without gaps and can be repeated several times using the same transmission power. Frequency hopping is applied to symbol group granularity, i.e. each symbol group is transmitted on a different subcarrier. To serve UEs in different coverage classes, the NB-IoT network can configure up to 3 NPRACH resource configurations in a cell. In each configuration, a repetition value from the set {1, 2, 4, 8, 16, 32, 64, 128} is specified for repeating a basic preamble [35]. In this model, we consider a single repetition value N T as [16], where the channel resources assignment of NPRACHs only takes place at the beginning of each transmission time interval (TTI, i.e., a time slot) as shown in [16].

D. Energy Harvesting Model
We assume that each IoT device is supported solely by the energy harvested from the surrounding environment (e.g., solar, kinetic, wind) and is equipped with a rechargeable battery with the finite capacity to buffer energy. Note that the assumption of a finite energy buffer is realistic since it is not possible to have infinite energy within an IoT device with limited physical dimensions. We model the energy arrival process of a randomly chosen IoT device in one TTI as an independent Poisson process with intensity µ 0 . This assumption is based on the fact that most energy harvesting modules contain small sub-modules harvesting energy independently, e.g., small solar cells harvesting energy in a solar panel, where the net energy harvested can be argued to be a Binomial process, which approaches to the Poisson process in the limit when the number of sub-modules grows large. This assumption is not uncommon, e.g., see [36] [37].
Since the energy arrivals are random and the energy storage capacities are finite, there is some uncertainty associated with whether the IoT device has enough energy to serve itself at a particular time. Under such a constraint, the IoT device needs to be kept OFF, and be allowed to recharge before it has sufficient energy to serve itself in a given time slot. Thus, at any given time, an IoT device can be in either of the two operational states: ON or OFF. In this paper, the decision to toggle the operational state of an IoT device, i.e., turn ON or OFF, is taken by the device independently, and it is not influenced by the operational states of other IoT devices. For example, one IoT device may decide to turn OFF if its current energy level reaches below a certain predefined level, and to turn ON after harvesting enough energy over the other threshold.
We define the fraction of time that a randomly chosen IoT device remains ON as the energy availability η 0 . In order to obtain η 0 , we first need to characterize how the energy available at the IoT device changes over time. We model the available energy of an IoT device as a finite-state continuous-time Markov process (CTMC). In particular, let {M (t) : t 0} be a stationary, homogeneous, and irreducible Markov process with state space M = {m 1 , m 2 , · · · , m i } that specifies the energy state at time t. For this setup, we define the energy required of a randomly chosen IoT device in each repetition E 0 as the unit energy [38]. Thus, the real energy state at a randomly chosen IoT device is directly proportional to the E 0 and the proportionality implies the maximum repetition value the IoT device can support. Without loss of generality, discretizing the real energy state of the IoT device by dividing by E 0 , we get our state space of a randomly chosen IoT device as M = {0, 1, · · · , M 0 }, in which M 0 = E/E 0 and E is the real energy storage capacity of the randomly chosen IoT device. Thus, the energy state m ∈ M of the IoT device implies the maximum repetition value the IoT device can support with the energy. As such, the temporal dynamics of the battery energy levels can be modeled by the CTMC, in particular, the birth-death process illustrated in Fig. 2. When the IoT device is ON , the energy increases according to the energy harvesting rate µ 0 units energy per second and decreases at a depletion rate of ν 0 units energy per second. Note that the data transmission after a successful RACH can be extended following the analysis of RACH success probability. Since the main focus of this paper is analyzing the contention-based RACH in the NB-IoT network, we assume that the actual intended packet transmission is always successful (i.e., the data transmission success probability is one) if the corresponding RACH succeeds and that the RACH either always fails as shown in Fig. 3(a) or succeeds as shown in Fig. 3(b) to obtain the lower bound and upper bound 1 of the energy availability η 0 . Thus, we have is the energy required for the RACH in each repetition; E DA 0 is the energy required for the data transmission in each repetition. Thus, we have incorporated the energy depletion rate for failure RACH and successful RACH cases separately as follows: where P is the transmission power of each IoT device and A a is the non-empty probability (i.e., IoT device data buffer is non-empty) described in Section IV. The main notations of the proposed protocol are summarized in TABLE I. 1 In our model, we take into account the RACH in one TTI with a single packet sequence transmission as shown rather than multiple TTI transmissions. We need to know the energy availability in the first TTI. To calculate energy availability for simplicity, we assume that the RACH either always fails as shown in Fig .3(a) or succeeds as shown in Fig. 3(b) to obtain the lower bound and upper bound of the energy availability η 0 . The results of our work are the foundation of further multiple TTIs analysis.  The Rayleigh fading channel power gain r The distance between an IoT device and its associated BS α The path-loss exponent N T The RACH repetition value µ 0 The energy arrival rate E 0 The energy required for each repetition E The real energy storage capacity ν 0 The energy depletion rate The energy required for the RACH in each repetition The energy required for the data transmission in each repetition Aa The non-empty probability P The The noise power I 0 The aggregate interference of the typical IoT device Z D The set of interfering IoT devices γ th The SINR threshold L The number of available preambles are reserved for the contention-based RACH λ Da The density of active IoT devices choosing the same preamble c c = 3.575 is a constant η 0 The energy availability

III. ENERGY AVAILABILITY OF IOT DEVICES
We assume that the energy harvesting processes are independent among the IoT devices, which ensure the independence of the current operational state (ON or OFF) of the IoT devices. An IoT device toggles its operational state solely on its current energy level. Essentially, we focus on a general strategy S{M min We note that, it is strictly sub-optimal when the IoT device toggles to OFF at M min  Table II. In our work, the IoT device will be  For the brevity of exposition, we denote this strategy by S{N T }. In practical, the probability that the IoT device is available may be different for each device due to the differences in the capability M 0 of the energy harvesting modules and the repetition value N T . However, we limit ourselves to the IoT device with the same repetition number N T and storage capability M 0 with the same energy availability probability.
For strategy S{N T }, the time for which an IoT device remains in the ON state after it toggles from the OFF state is given by T on 0 {N T }, and the time for which it remains in the OFF state after toggling from the ON state is given by For notational simplicity, the cutoff value {N T } will be dropped wherever appropriate. It is noticed that both T on 0 and T off 0 are general random variables. We first obtain T on 0 as [39] T on where ξ 0 (τ ) denotes the current energy level of a randomly chosen IoT device at time τ . For this setup, the energy availability of a randomly chosen IoT device depends only on the means of T on 0 and T off 0 , which is shown as and E[T on 0 ] is the mean time the randomly chosen IoT device remains in the ON state, and E[T off 0 ] is the mean time it remains in the OFF state.
Proof. Let {T on 0 (k)} and {T off 0 (k)} be the sequences of the kth cycle of ON and OFF times, respectively. The availability can now be described by the fraction of time the randomly chosen IoT device remains in the ON state as  [40], we have E[T off 0 ] = N T /µ 0 (i.e., the time required to harvest N T units of energy), which is the sum of N T exponentially distributed random variables, each with mean 1/µ 0 . Substituting into (3), we obtain To derive the mean ON time E[T on 0 ], we first need to obtain the transmission matrix Q 0 for the birth-death process corresponding to the randomly chosen IoT device. According to the Kolmogorov differential equations [41], where the first column corresponds to the energy state 0 and the states are in the ascending order. Then we obtain the following Lemma 1. from [39]. where Proof. See Appendix A.
Substituting (8) into (7) and then plugging into m = N T , we have the mean ON time for strategy S{N T } with energy storage capacity M 0 as Submitting (9) into (5), we obtain the energy availability in the following Theorem 1.. Theorem 1. (Energy Availability). The energy availability of a randomly chosen IoT device is given by In (10), it can be shown that the energy availability increases with increasing µ 0 and M 0 . For illustration, the relationships between η 0 and N T , M 0 , µ 0 are analyzed in Section V.

IV. RACH TRANSMISSION SUCCESS PROBABILITY
In the NB-IoT repetition scheme, an active IoT device will repeat the same preamble N T times (i.e., the dedicated repetition value). In each repetition, a preamble is composed of four symbol groups transmitted without gaps, where the first preamble symbol group is transmitted via a sub-carrier determined by pseudo-random hopping (i.e., the hopping depends on the current repetition time and the Narrowband physical Cell ID, a.k.a NCellID [10]), and the following three preamble symbol groups are transmitted via sub-carriers determined by the fixed size frequency hopping [42]. This frequency hopping algorithm is designed in a way that different selections of the first subcarrier lead to hopping schemes which never overlap. Specifically, if two or more IoT devices chose the same first sub-carrier in a single RAO, the following sub-carriers (i.e., in the same RAO) would be same, due to that these two hopping algorithms lead to one-to-one correspondences between the first sub-carrier and the following sub-carriers (i.e., these IoT devices either collide on the full set or not collide at all in a single RAO). In this setup, the RACH success refers to the preamble being successfully transmitted to the associated eNB (i.e., received SINR is greater than the SINR threshold) and no collision occurs (i.e., no other IoT devices successfully transmit a same preamble to the typical eNB simultaneously).
For an IoT device to be able to initiate an uplink transmission, the energy harvested by this device should be sufficient to perform RACH and data transmission. We have defined the energy availability of a randomly chosen IoT device after harvesting enough energy as η 0 in Theorem 1. In this section, we first formulate the SINR outage condition to drive the preamble transmission success probability and then facilitate the analysis of the RACH success probability.

A. SINR Definition
Recall that each IoT device transmits a randomly chosen preamble to its associated eNB to request for channel resources, where different preambles represent orthogonal subchannels, and thus only IoT devices choosing the same preamble have correlations. The RACH analysis in this work needs to take into account both the inter-and intra-cell interference 2 . The received power at a typical eNB from a randomly chosen IoT device of interest is therefore P 0 = P h 0 r −α 0 , where h 0 and r 0 are the channel power gain and the distance from the typical IoT device to its associated eNB respectively. Using the received power over the link of interest and the interference power, the SINR received at the typical eNB at the origin can be written as where σ 2 is noise power, and I 0 is aggregate interference of the typical IoT device with In (12), h j and r j are channel power gain and the distance from the interfering IoT devices to the typical eNB, P is the transmission power of the IoT devices, and Z D is the set of interfering IoT devices for the typical IoT device. We note that only the active IoT devices choosing the same preamble will generate interference. The density of active IoT devices choosing the same preamble is obtained as follows.
Remark 1. (The density of active IoT devices choosing the same preamble). Note that inactive IoT devices (those without enough energy or data packets in buffer) do not attempt RACH, such that they do not generate interference. According to the repetition scheme mentioned earlier, each active IoT device will contend on all L = 48 sub-carriers due to the single repetition value configuration, and thus each preamble has an equal probability (1/L) to be chosen. As only active IoT devices will try to request uplink channel resources, we define the non-empty data packets probability of each IoT device A a ∈ [0, 1] following a Bernoulli process. We also have the energy availability η 0 of each IoT device, then according to the thinning process [44], the density of active IoT devices choosing the same preamble can be expressed as 2 In LTE, the PRACH root sequence planning is used to mitigate intercell interference among neighboring BSs (i.e., neighboring BSs could use different roots to generate preambles) [8]. However, as in [14] [43], we focus on providing a general analytical framework of cellular networks considering both the inter-and intra-interference without using PRACH root sequence planning. That is to say, we consider intra-cell interference due to the fact that the IoT devices in the same cell may choose the same preamble and we consider the inter-cell interference due to the fact that the IoT devices in different cells share the preamble sequence pool among eNBs. We also obtain the approximation results of the networks with perfect PRACH root sequence planning (e.g., no inter-cell interference from any cells) in Lemma. 3. But the extension taking into account PRACH root sequence planning will be treated in future works.

B. RACH Success Probability
We formulate the RACH success probability under both SINR outage and collision conditions. We perform the analysis on an eNB associating with a randomly chosen active IoT device in terms of the RACH success probability. The RACH success refers to the preamble being successfully transmitted to the associated eNB (i.e., received SINR is greater than the SINR threshold) and no collision occurs (i.e., no other IoT devices successfully transmits a same preamble to the typical eNB simultaneously). First, we formulate the SINR outage condition. The typical IoT device transmits a preamble successfully if any repetition successes, and in a single repetition, a preamble is successfully received at the associated eNB if its all four received SINRs are above the SINR threshold γ th . Thus, the preamble transmission success probability of a randomly chosen IoT device under N T repetitions is expressed as I is the probability that all four (a preamble consists of four preamble symbol groups) time-correlated preamble symbol groups in the n T th repetition are successfully transmitted, II is the probability that all N T repetitions of a preamble transmission are failed, and In (15), γ th is the SINR threshold, and SINR n T ,1 (r 0 ), SINR n T ,2 (r 0 ), SINR n T ,3 (r 0 ) and SINR n T ,4 (r 0 ) are the received SINRs of the four symbol groups in the n T th repetition of the typical IoT device. Based on the Binomial theorem, the preamble transmission success probability in (14) can be rewritten as where is the binomial coefficient, and P 0 [θ 1 (r 0 ), θ 2 (r 0 ), · · · , θ n T (r 0 )] is the probability that all of 4 × n T (a preamble consists of four preamble symbol groups) time-correlated preamble symbol groups are successfully transmitted. We note that the preamble transmission success probability in (16) depends on the transmission distance r 0 . Taking into account that each IoT device associates to its geographically nearest eNB, r 0 is the minimum distance between the eNB and the typical IoT device. The PDF of the shortest distance between any point BS and the IoT device with radius r 0 is [45] where ε = 1 when λ Da λ B and ε = 1.25 when λ Da λ B . For ease of presentation, we set l = 4 × n T , and the probability that all of 4 × n T preamble symbol groups are successfully transmitted is presented in the following Lemma 2.
Lemma 2. The probability that all of 4×n T received SINRs at the eNB from a randomly chosen IoT device exceed a certain threshold γ th is expressed as Proof. See Appendix B.
Next, we formulate the RACH success probability taking into account both the SINR outage and the collision. The RACH success probability is represented in the following Theorem 2.
Theorem 2. In the energy harvesting NB-IoT network, the RACH success probability of a randomly chosen IoT device is derived as where and Part I the Probability Mass Function (PMF) of the number of intra-cell interfering 3 IoT devices for a typical BS N = n derived following [47,Eq.(3)], where c = 3.575 is a constant related to the approximate PMF of the PPP Voronoi cell and Γ(·) is the gamma function. Part II is the preamble transmission success probability of a randomly chosen IoT device obtained by substituting (18) into (16). Part III is the preamble transmission failure probability that the transmissions from other n intra-cell interfering IoT devices are not successfully received by the BS, i.e., the non-collision probability of the typical IoT device conditioning on n.
Remark 2. It is evident from (21) that the transmission success probability (II in (19)) of the typical IoT device increases, whereas the non-collision probability (III in (19)) decreases with increasing the repetition value N T and decreasing the the received SINR threshold γ th . Therefore, there exists a tradeoff between transmission success probability and non-collision probability. For illustration, the relationship among RACH access success probability (P 0 ), the transmission success probability (P 0 with III=1), and the non-collision probability (P 0 with II=1) versus repetition value N T and the received SINR threshold γ th is shown in Fig. 4. Note that the RACH success probability of a randomly chosen IoT device in networks with perfect PRACH root sequence planning could be obtained by only considering the intra-cell interference. Considering that the practical Voronoi cells do not have a constant radius, we use the average radius D = 1/ √ πλ B [18] to approximate it 4 . Thus, the RACH success probability is given in the following Lemma 3. Lemma 3. The RACH success probability of a randomly chosen IoT device in NB-IoT networks with perfect PRACH 4 Note that for the networks with perfect PRACH root sequence planning, it is difficult to characterize the radius of the practical cell. We use D = 1/ √ πλ B to approximate the radius of the practical cell and the results depend on the practical cell shape. Our results have a good match when choosing a proper deployment area and λ B as shown in Fig. 10 and Fig. 12. root sequence planning is derived as where P[N = n] is given in (20) and

V. SIMULATION AND DISCUSSION
In this section, we verify our analytical results by comparing the theoretical RACH success probabilities with the results from Monte-Carlo simulations. For numerical verification, we compute the RACH success probability from Monte-Carlo simulations as follows. We simulate the spatial model described in Section II in MATLAB. The eNBs and IoT devices are deployed via independent HPPPs in a 2 × 10 4 km 2 circle area. Each IoT device associated with its nearest eNB. The IoT devices and the eNBs remain spatially static during a TTI. The channel fading gains between the IoT devices and eNBs are modeled by exponentially distributed random variables. Unless otherwise stated, we set P = 0.02 w, T r = 6 ms, T g = 31 ms, µ i = 0.05, λ B = 0.1 eNBs/km 2 , λ D = 10 2 IoT devices/km 2 , γ th = 20 dB, α = 4, A a = 0.001, the bandwidth of a subcarrier is BW= 3.75 kHz, and thus the noise is σ 2 = −174 + 10log10(BW) = −138.3 dBm. For each realization of this setup, the uplink communication has declared a success if 1) the calculated SINR exceeds a predetermined threshold γ th and 2) other uplink communications using the same preamble do not exceed γ th . In all figures of this section, "Analytical" and "Simulation" are abbreviated as "Ana." and "Sim.", respectively.  (10). We assume that the RACH in each TTI succeeds or fails on the full set (i.e., the RACH always succeeds or all the RACH always fails), and then we derive the upper bound and the lower bound of the energy availability. As such, we could give the availability region for various values of energy availability.

A. Analysis of the Energy Availability
We first observe that the energy availability of the IoT device increases with increasing the preamble repetition value N T at first and then remains unchanged. This is due to the fact that for the same storage capacity, increasing the repetition value increases the operation time of the IoT device, i.e., the IoT device spends more time in ON state. In addition, since M c 0 = N T , increasing the repetition value N T , i.e., increasing the cutoff value, results in more time needed for the IoT device to harvest sufficient energy, i.e, stay in OFF state for more time before the transmission. As such, to obtain a higher energy availability, a larger repetition value N T is needed, but if the repetition value is overestimated, the IoT device will waste the potential resource for data transmission and lead to lower resource efficiency. That is to say, the repetition value needs to be optimized. Interestingly, we also observe that for energy storage capacity much larger than cutoff value, the energy availability approaches a specific value for different preamble repetition values, e.g., η 0 = 0.3 for the lower bound when and η 0 = 0.92 for the upper bound when M 0 − M c 0 ≥ 160, which reveals that this setup is surprisingly reliable if it is designed properly, despite the randomness in the energy harvesting.  (10). We first observe that the energy availability of the IoT device increases with increasing the energy harvesting rate µ 0 when the storage capacity M 0 is not large enough, e.g.,M 0 = N T . This is due to the fact that for the same storage capacity and the cutoff value, increasing the energy harvesting rate results in less time needed for the IoT device to harvest sufficient energy, i.e, stay in OFF state for less time before the transmission. Interestingly, we also observe that for energy storage much larger than cutoff value, e.g., M 0 − M c 0 ≥ 10, the energy availability approaches a specific value for different preamble repetition values. That is to say, if it is designed properly, the energy availability is independent of N T .

B. Validation of the RACH success probability
In this section, we simulate an NB-IoT network and evaluate the RACH success probability based on the energy availability analyzed above. Fig. 7 plots the RACH success probability of a randomly chosen IoT device versus the energy availability η for various preamble repetition values N T . We observe that the RACH success probability deteriorates as the energy availability of IoT devices increases. This is due to the fact that increasing energy availability increases the number of active devices, which leads to lower received SINR and a higher probability of collision. In order to obtain fundamental insights on the RACH success probability due to the preamble repetition value N T , in the following analysis, we use unchanged  Fig. 8 plots the RACH success probability of a randomly chosen IoT device versus the transmission power P for various preamble repetition values N T . We first observe a good match between the analysis and the simulation results, which validates the accuracy of the developed mathematical framework. As expected, we observe that the RACH success probability increases as the transmission power of IoT devices increase. This can be explained by the reason that whilst increasing the transmission power P leads to higher interference power, it also leads to increased received signal power, thereby improves the overall SIR and hence the RACH success probability. It is worth noting that the RACH success probability increases faster at first (e.g., when P ≤ 0.1) and then gradually becomes steady, which reveals that there is a limit value of the cell maximum transmission power. Interestingly, we observe that the RACH success probabilities with a higher repetition value, e.g., N T = 8 become stable earlier than those with a lower repetition value, e.g., N T = 1, due to that the higher chance of the RACH to succeeds in higher repetition value case.  Fig. 9 plots the RACH success probability of a randomly chosen IoT device versus the SINR threshold γ th for various preamble repetition values N T . As expected, the RACH success probability degrades with an increase in the SINR threshold. According to Fig. 4, increasing SINR threshold γ th leads to lower preamble transmission success probability but higher non-collision probability, thereby decreases the overall RACH success probability. There is a tradeoff between preamble transmission success probability and non-collision probability.  Fig. 10(a) and Fig. 10(b) plot the RACH success probability of a randomly chosen IoT device versus the density ratio λ D /λ B in the network with and without inter-cell interference, respectively. We first observe that the RACH success probability decreases with the increase of the density ratio between IoT devices and BSs (λ D /λ B ), which is due to the following two reasons: 1) increasing the number of IoT devices generating interference leads to lower received SINR at the eNB; 2) increasing the number of IoT devices leads to higher probability of collision. In addition, increasing repetition value increases the RACH success probability due to that it offers more opportunities to re-transmit a preamble with the time and frequency diversity.  Fig. 11(a) and Fig. 11(b) plot the RACH success probability of a randomly chosen IoT device versus the density ratio λ D /λ B in the network with light traffic (A a = 0.001, λ Da λ B , and ε = 1) and heavy traffic (A a = 0.015, λ Da λ B , and ε = 1.25), respectively. In the light traffic scenario, RACH success probability increase when increasing the repetition value. However, in the heavy traffic scenario, the RACH success probability cannot improve much when the repetition value increased. That is to say, the repetition scheme can efficiently improve the RACH success probability in a light traffic scenario, but only slightly improves that performance with very inefficient channel resource utilization in a heavy traffic scenario. Fig. 12: RACH success probability of a randomly chosen IoT device versus the densities ratio λ D /λ B and SINR threshold γ th . Fig. 13: Repetition efficiency of a randomly chosen IoT device versus the densities ratio λ D /λ B and SINR threshold γ th . Fig. 12 and Fig. 13 plot the RACH success probability and the Repetition efficiency (i.e., ζ = P 0 /N T ) of a randomly chosen IoT device versus the density ratio λ D /λ B and SINR thresholds γ th , respectively. Obviously, increasing the repetition value increases the success opportunities of RACH, but decreases the repetition efficiency, which reveals that if the repetition value is overestimated, the IoT device will waste the potential resource for data transmission and lead to lower resource efficiency. In addition, we observe that the repetition efficiencies decrease seriously when increasing the repetition value from N T = 1 to N T = 2, especially for the lower density ratio. This is due to the fact that when there are fewer active IoT devices (e.g., λ D /λ B = 10 4 ) to contention for the same resources, it is likely to succeed with small repetition values, which contributes to a relatively high channel resource utilization.

VI. CONCLUSION
In this paper, we have presented a comprehensive RACH success probability analysis for the NB-IoT network with energy harvesting from natural resources. We have derived a closed-form expression of the energy availability and analyzed the results for some realistic strategies. We found that energy availability remains unchanged if the network is properly designed despite randomness in energy harvesting. We further analyzed the RACH under the repetition transmission scheme in the NB-IoT system based on energy availability. We derived the exact expression for the RACH success probability under time-correlated interference. Different from existing works, we considered both SINR outage and collision from the network point of view in the NB-IoT network, where each IoT device adopts fixed transmission power. Our results have shown that increasing the repetition value increases the RACH success probability but decreases the repetition efficiency. In the light traffic scenario, RACH success probability can meet the requirement with small repetition values. However, in the heavy traffic scenario, very high repetition value leads to a low channel resources utilization. Our results also have shown that there is an upper limit on transmission power and too large transmission power will waste energy. where A r and E r denote the rth column of A 0 −1 and I respectively, r = 1, 2, · · · , M 0 . Then (A.1) can be rewritten in the form In order to calculate A 0 −1 , we begin by considering the elements of the first column A 1 = (a 11 , a 21 , · · · , a M01 ) T of A 0 −1 . That means
The Laplace transform of the aggregate interference is characterized according to the definition of L I0 (s) = E I0 (e −sI0 ) as