Multi-step-ahead forecasting of bike-sharing demand using multilayer perceptron model with additional timestamp features

Data preparation

To assess the effectiveness of the prediction models, we gathered data on bike-sharing reservations in Korea. This reservation data spans from January 2020 to December 2020, encompassing a one-year period, and was obtained from Seoul Public Data (Seoul Bike Sharing, 2020). The original reservation data, as illustrated in Fig. 2, contains various columns including bike ID, starting datetime, departure ID, ending datetime, destination ID, and duration (in minutes). To construct the dataset for hourly bike demand, the reservation data needed to be reformatted, counting the total number of rented bikes for each hour. Subsequently, timestamp information was utilized to extract details such as the hour, day, and month. The resulting dataset for this study includes information on the hour, day, month, and the number of rented bikes. Table 2 provides a breakdown of the dataset description in this research.

The next step is Exploratory Data Analysis (EDA), a crucial phase aimed at gaining insights into the characteristics and patterns present in the bike-sharing demand dataset. EDA involves a comprehensive examination of the data through statistical summaries, visualizations, and graphical representations. This process helps uncover trends, distributions, and patterns within the dataset. Identifying such patterns is instrumental in informing subsequent steps of the modeling process, guiding feature selection, and enhancing the overall interpretability of the model.

Figure 3 illustrates the count of rented bikes, depicting the average count across days, months, and hours throughout the week. The visualizations reveal a discernible temporal pattern in the rented bike numbers, highlighting the significance of these temporal dependencies for more accurate rented bike predictions. Specifically, Fig. 3A showcases the overall rented bike count, while Fig. 3B presents the daily average, revealing notable variability in bike counts at different hours and days. Figure 3C underscores that the average number of rented bikes is considerably lower during the winter months (January, February, and December) in South Korea. Conversely, peak bike rentals occur during the spring months (April, May, and June) and fall months (October and November). Furthermore, Fig. 3D displays the average rented bikes categorized by day of the week. It is evident from Fig. 3D that the average trip duration on weekdays peaks during the morning hours (7 to 9 AM) and evening hours (5 to 7 PM). On weekends, a different trend emerges, with the majority of rented bikes starting in the afternoon, specifically between 1 and 8 PM. This temporal analysis provides valuable insights into the dynamics of bike-sharing demand, laying the foundation for a more nuanced understanding of how time-related factors influence bike rental patterns. These findings are essential for refining predictive models and optimizing resource allocation in bike-sharing systems.

The temporal analysis depicted in Fig. 3 not only uncovers seasonal variations in bike rental patterns but also offers implications for operational considerations in bike-sharing systems. The observed peaks in bike rentals during specific hours and seasons suggest a need for dynamic resource management. For instance, during high-demand periods such as the morning and evening rush hours on weekdays, bike-sharing systems may benefit from strategically redistributing bikes to meet increased demand in specific locations. Moreover, demand notably rises during the spring and summer months, reflecting more favorable weather and increased recreational usage. Furthermore, the comprehensive exploration of temporal dependencies in bike-sharing demand, as illustrated in Fig. 3, not only aids in predicting rented bike counts more accurately but also provides actionable insights for the strategic management and optimization of bike-sharing systems. The integration of these findings into operational practices can enhance user experiences, increase system efficiency, and contribute to the sustainable growth of bike-sharing initiatives.

Feature extraction

In the feature extraction phase, the time series dataset undergoes transformation into a collection of paired inputs (X) and corresponding desired outputs (y). The prediction process utilizes the “direct method,” which involves independently forecasting each horizon without reliance on other predictions (Ben Taieb et al., 2012; Bontempi, Ben Taieb & Le Borgne, 2013). The direct strategy autonomously learns each horizon (denoted as $H$) and models the corresponding function, denoted as $f_h$

(1) $$y_{t+h}=f_h\left(y_t,\dots ,y_{t-n+1}\right)+w.$$with $t\in \left\{n,\dots ,N-H\right\}$ and $h\in \left\{1,\dots ,H\right\}$ and provides a multi-step prediction by combining the forecasts for each horizon ( $H$). To implement this direct strategy, the sliding window approach can be employed to segment the time series dataset, as illustrated in Fig. 4. In this scenario, input data is generated by assembling windows with a size of n = 24 for learning, aiming to predict the count of rented bikes for the subsequent h = 6. Given that the count of rented bikes was recorded every hour in our study, this prediction model utilizes the last 24 h of historical data (as features) to forecast the count of rented bikes for the next 6 h. In Fig. 4, the windows (indicated by red lines) are gathered and presented as an input matrix X, while the following 6 h of values (depicted by green lines) serve as an output vector y. The pairs of inputs and outputs are derived by shifting the window one-time step ahead on each iteration. Two parameters need consideration in this context: the window size, n, and the prediction step ahead, h. Subsequently, the collection of paired inputs and desired outputs is divided into a training set and a testing set.

By adhering to this process, a series of paired inputs and corresponding desired outputs can be generated, allowing the supervised learning model to glean insights from these input-output pairs. Here, C denotes the collection of rented bike counts, and $c_i$ signifies a specific count of rented bikes within the set, where $c_i\in$C and i = 1, 2…, N. The parameter N represents the total number of data entries. Ultimately, with the given $n$ previous values (or window size) and the forecasting horizon, $h$, the input matrix $X$ can be formulated by constructing a $\left[\left(N-n-h+1\right)\times n\right]$ matrix of input data

(2) $$X=\left[\begin{array}{cccc}{c}_1& \dots & c_{n-1}& c_n\\ ⋮& \ddots & ⋮& ⋮\\ c_{N-n-h}& \dots & c_{N-h-2}& c_{N-h-1}\\ c_{N-n-h+1}& \dots & c_{N-h-1}& c_{N-h}\end{array}\right]$$and the $\left[\left(N-n-h+1\right)\times 1\right]$ output vector

(3) $$y=\left[\begin{array}{c}{c}_{n+h}\\ ⋮\\ c_{N-1}\\ c_N\end{array}\right].$$

Figure 3 indicates a notable correlation between timestamps information and the count of rented bikes. Additionally, we assess the relationship between timestamp attributes and the count of rented bikes. To investigate this association, we utilized Pearson’s correlation coefficient, which ranges from −1 to +1. A negative or positive value denotes a negative or positive correlation, and a higher absolute value signifies a more robust correlation. Attributes demonstrating a substantial correlation with the output class can be utilized as input features to improve the accuracy of the prediction model. Figure 5 demonstrates that the hour of the day, day of the week, and month of the year exhibit a positive correlation with the output.

Lastly, we introduce timestamp attributes as supplementary features for the n previous values generated through the sliding window approach, encompassing elements such as the hour of the day, day of the week, and month of the year. The n previous values from the hour of the day (t), day of the week (d), and month of the year (m) are integrated with the n previous values of the counted rented bikes (c). Ultimately, given the n previous values or window size, the forecasting horizon (h), and the total data entries (N), the updated input matrix X can be formulated by creating a $\left[\left(N-n-h+1\right)\times \left(n\times 4\right)\right]$ input matrix as can be seen in Formula (4).

(4) $$X=\left[\begin{array}{ccccccccc}{c}_1& t_1& d_1& m_1& \dots & c_n& t_n& d_n& m_n\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ c_{N-n-h}& t_{N-n-h}& d_{N-n-h}& m_{N-n-h}& \dots & c_{N-h-1}& t_{N-h-1}& d_{N-h-1}& m_{N-h-1}\\ c_{N-n-h+1}& t_{N-n-h+1}& d_{N-n-h+1}& m_{N-n-h+1}& \dots & c_{N-h}& t_{N-h}& d_{N-h}& m_{N-h}\end{array}\right].$$

Proposed MLP model

In this study, a MLP model was employed for the prediction of bike sharing demand. The MLP belongs to the category of feedforward Artificial Neural Networks (ANNs), characterized by one input layer, one or more hidden layers, and one output layer. To train the MLP, a backpropagation algorithm was utilized (Han, Kamber & Pei, 2012; Rumelhart, Hinton & Williams, 1986). Our network is fully connected, meaning each unit receives connections from all units in the preceding layer. Consequently, each unit possesses its own bias, and there exists a weight for every pair of units in two consecutive layers. The net input was computed by multiplying each input by its corresponding weight and then summing the results. Each unit in the hidden layer received a net input, which was subsequently subjected to an activation function. Therefore, with input units are denoted as $x_j$, output unit as y, units in the lth hidden layer represented as $h_i^{\left(l\right)}$, and $\phi$ as the activation function, the network computation can be expressed as follows.

(5) $$h_i^{\left(1\right)}={\phi }^{\left(1\right)}\left(\sum _j{w}_{ij}^{\left(1\right)}{x}_j+b_i^{\left(1\right)}\right)$$

(6) $$h_i^{\left(2\right)}={\phi }^{\left(2\right)}\left(\sum _j{w}_{ij}^{\left(2\right)}{h}_j^{\left(1\right)}+b_i^{\left(2\right)}\right)$$

(7) $$y_i={\phi }^{\left(3\right)}\left(\sum _j{w}_{ij}^{\left(3\right)}{h}_j^{\left(2\right)}+b_i^{\left(3\right)}\right).$$

The backpropagation technique involves assessing the predicted outcome against the target value and adjusting the weights for each training tuple to minimize the mean squared error between the predicted and target values. This iterative process is repeated multiple times to achieve optimal weights, thereby enabling optimal predictions for the test data. Figure 6 depicts the suggested MLP model for forecasting the count of rented bikes, using inputs comprising the rented bikes’ history over the last 24 h, along with attributes such as the time of day, day of the week, and month of the year.

The selection of the number of lags or window size in a forecasting model is crucial as it directly impacts the prediction error. The appropriate choice of lags or window size ensures capturing relevant historical information, enhancing the model’s ability to make accurate predictions. Therefore, in this investigation, we explore the most effective window size for the proposed MLP model. Initially, we partitioned the training data (0.75 of the dataset) into an additional training set and a validation set (0.25). The MLP model underwent evaluation using various window sizes, including 12, 18, 24, 30, and 36, alongside different prediction horizons (PH) such as 1, 3, 6, 12, and 24. The outcomes of the experiment, as depicted in Fig. 7, revealed that a window size of 24 was optimal in the validation set, demonstrating lower errors for PH = 1, PH = 3, and PH = 6. Interestingly, our findings indicated that a higher window size, specifically n = 36, exhibited the lowest Root Mean Square Error (RMSE) for PH = 12 and PH = 24.

Choosing the appropriate window size or number of lags in forecasting is a critical decision that significantly influences the accuracy and effectiveness of predictive models. A window size that is too small may result in overlooking important patterns and trends in the historical data, leading to an oversimplified model that fails to capture the complexity of the underlying patterns. On the other hand, an excessively large window size may introduce unnecessary noise and irrelevant information, hindering the model’s ability to discern relevant signals for prediction. Striking the right balance is essential, as an optimal window size ensures that the model leverages sufficient historical context without being overwhelmed by excessive data, ultimately contributing to more accurate and reliable forecasting results. Hence, a window size of 24 was chosen in our study as the input features for our proposed machine learning models.

The model evaluation employed the holdout method, dividing the dataset into two segments: the initial 75% for training and the remainder for testing. The forecasting models were implemented using Python V3.7.3, Scipy V1.3.2, and Scikit-learn V0.22.1 (Pedregosa et al., 2012), with optimized hyperparameters as can be seen in Table 3. Performance metrics, including RMSE, Mean Absolute Error (MAE), and Coefficient of Determination ( ${\mathrm{R}}^2$), were utilized to assess the forecasting models. RMSE measures the average magnitude of differences between predicted and actual values (number of rented bikes), offering a comprehensive evaluation of prediction accuracy. MAE quantifies the average absolute differences between the predicted and actual number of bikes, providing a straightforward representation of prediction errors. Lastly, ${\mathrm{R}}^2$ measures the proportion of variance in the dependent variable predictable from the independent variable, indicating the goodness-of-fit of the regression model.

Performance comparison of prediction models

The proposed MLP model was tested on the Seoul bike-sharing dataset and exhibited favorable outcomes in reducing prediction errors compared to alternative forecasting models. This research involved a comparative analysis between the proposed model and various machine learning-based prediction algorithms, including Support Vector Regression (SVR), K-Nearest Neighbour (KNN), Decision Tree (DT), Adaptive Boosting (AdaBoost), Random Forest (RF), and Linear Regression (LR). The effectiveness of each model was assessed based on lower RMSE and MAE, along with higher ${\mathrm{R}}^2$ values in the testing sets. The predictive input features comprised the past 24 h of rented bikes to predict the next 1, 3, 6, 12, and 24 h of bike rentals. The study demonstrated that utilizing the preceding 24 h of rented bikes as input features could minimize prediction errors on the validation sets. Additionally, the proposed MLP model incorporated timestamp features (hour, day, month) as supplementary attributes, distinguishing it from other models that solely relied on the last 24 h of rented bikes for input.

Table 4 displays the average values of RMSE, MAE, and ${\mathrm{R}}^2$ from various forecasting models aimed at predicting the demand for bike sharing (rented bikes) across different prediction horizons (PHs). The proposed MLP model demonstrated superior performance, exhibiting the lowest RMSE and MAE, as well as the highest ${\mathrm{R}}^2$ in the testing sets compared to alternative prediction models. Specifically, the RMSE values for the proposed MLP model were 344.020, 719.812, 889.602, 920.725, and 959.997 for prediction horizons of 1, 3, 6, 12, and 24 h, respectively. Furthermore, the proposed model achieved the highest ${\mathrm{R}}^2$ (coefficient of determination) in comparison to other models, with values of 0.973, 0.882, 0.820, 0.807, and 0.790 for PHs of 1, 3, 6, 12, and 24 h, respectively. Conversely, both Decision Tree (DT) and AdaBoost exhibited the poorest performance across all prediction horizons compared to other models. Ultimately, the findings indicated that as the prediction horizon (PH) increased, the majority of prediction models produced higher RMSE and MAE but with lower ${\mathrm{R}}^2$ values in the testing sets.

To assess the effectiveness of the suggested MLP-based prediction model, graphical representations of the prediction outputs and original values from the testing sets can be examined. As illustrated in Fig. 8, the predicted values from the proposed MLP model are visually compared with the original hourly rented bikes in the testing sets. For this graphical comparison, we opted to analyze the initial 100 points from both the testing sets and the prediction output. The findings indicate that the prediction output from the proposed MLP model closely aligns with the testing sets or original data, particularly evident with a prediction horizon (PH) of 1 h (Fig. 8A). The model’s prediction output for a PH of 1 h exhibits the lowest error when contrasted with the models for PHs of 3, 6, and 24 h. Additionally, as the prediction horizon increases, the forecasting model experiences higher prediction errors. Specifically, for a PH of 24 h (Fig. 8E), the model generates a higher error (i.e., the disparity between predicted and original values) compared to the models for PHs of 12 (Fig. 8D), 6 h (Fig. 8C) and 3 h (Fig. 8B). Based on this preliminary experiment, it can be inferred that the proposed MLP model shows great promise and could be effectively applied for predicting hourly demand in real-world bike-sharing scenarios.

To further evaluate the predictive accuracy of the proposed model for PH = 1, a comprehensive residual analysis was conducted. The residual plot, presented in Fig. 9A, illustrates the relationship between residuals and predicted values, showing no discernible patterns or signs of heteroscedasticity. This indicates that the errors are randomly distributed, which is a desirable property in a well-performing predictive model. Additionally, the histogram of residuals in Fig. 9B reveals an approximately symmetric distribution with a slight positive skew. The descriptive statistics of the residuals are as follows: the mean residual is 11.990, indicating a negligible overall prediction bias; the standard deviation is 343.810, reflecting a moderate level of dispersion in the errors; the skewness is 0.400, suggesting mild asymmetry toward higher residual values; and the kurtosis is 10.100, indicating heavier tails than a normal distribution, which may reflect the presence of a small number of relatively large prediction errors. These results support the robustness of the model’s predictive performance and offer valuable insights into the distribution and behavior of its residuals.

To determine whether the proposed MLP model provides a statistically significant improvement in performance over the baseline models, a paired t-test was carried out based on the RMSE values. As shown in Table 5, the results indicate that the proposed MLP model significantly outperforms all competing models. Specifically, the differences between the MLP and each of the baseline models—including SVR (p = 0.0153), KNN (p = 0.0092), DT (p = 0.0007), AdaBoost (p = 0.00001), RF (p = 0.0036), and LR (p = 0.0090)—are statistically significant at the 0.050 level. These findings confirm that the MLP model consistently delivers superior predictive performance, providing more accurate forecasts across all time horizons compared to the other machine learning models evaluated in this study.

Furthermore, the proposed MLP model was evaluated against two classical statistical forecasting methods—ARIMA and Exponential Smoothing—over a 24-h prediction horizon, as can be seen in Table 6. The results demonstrate that the MLP model significantly outperforms both traditional approaches across all evaluated metrics. Specifically, the MLP achieved a RMSE of 959.997, a MAE of 633.722, and a coefficient of determination (R²) of 0.790. In comparison, the ARIMA model recorded an RMSE of 4,760.520, MAE of 3,908.360, and R² of 0.520, while the Exponential Smoothing model achieved an RMSE of 1,571.560, MAE of 1,318.520, and R² of 0.740.

These findings indicate that the MLP model provides more accurate and reliable forecasts than the ARIMA and Exponential Smoothing models. The superior performance of the MLP can be attributed to its ability to capture complex nonlinear patterns in the data, which traditional linear models may fail to address effectively.

The impact of additional timestamp features on prediction error

Our proposed MLP-based model incorporates the previous 24 h of extracted features from timestamps as supplementary input, distinguishing it from other prediction models that rely solely on the past 24 h of rented bikes as attributes. Consequently, it is essential to investigate the impact of these additional timestamp features on the prediction error of other models. Table 7 displays the RMSE of the prediction models in the testing sets with various feature types for different prediction horizons (PHs). The outcomes reveal that the combination of the last 24 h of rented bikes and timestamp features as attributes led to lower RMSE for all PHs compared to prediction models using only previous values of rented bikes as features—except for DT (6, 12, and 24 h PHs), AdaBoost (1 and 24 h PHs), and RF (6, 12, and 24 h PHs). In the context of short-term forecasting (i.e., 1 h), the combination of rented bikes and timestamp features produced lower RMSE compared to models with rented bike attributes, except for AdaBoost. Notably, in our experimental findings, the inclusion of additional timestamp features diminished the performance of DT, AdaBoost, and RF, resulting in higher RMSE for longer-term forecasting (i.e., 24-h PH) when compared to models without timestamp features.

As depicted in Table 7, the results indicate a notable enhancement in performance by incorporating additional timestamp features, particularly for SVR and Linear Regression. The average reduction in RMSE for these models is 92.738 and 61.859, respectively, when compared to models without timestamp features. However, the impact of additional timestamp features is less favorable for the DT-based algorithms, resulting in an average RMSE increase of 108.349 compared to the non-timestamp DT model. Lastly, the proposed MLP, when augmented with extra timestamp features, has significantly lowered the RMSE by 3.655, 55.198, 49.036, 23.364, and 4.779 for prediction horizons of 1, 3, 6, 12, and 24 h, respectively, in contrast to the MLP model without timestamp features.

Comparison with previous study and other deep learning models

A prior investigation successfully demonstrated the prediction of bike-sharing demand using various models and feature types. To enhance the accuracy of hourly bike-sharing demand predictions, a Gradient Boosting Machine model was proposed in a previous study, incorporating weather and date information as input features (Sathishkumar, Park & Cho, 2020a). This model utilized the Seoul bike-sharing 2018 dataset, obtained from a Korean bike-sharing company’s reservation data spanning December 2017 to November 2018. Features in the dataset included the number of rented bikes, weather conditions, and date information for each hour.

For our experiment, we additionally employed the publicly available Seoul bike-sharing 2018 dataset from the aforementioned study and conducted machine learning-based forecasting. The input features for our proposed MLP model consisted of the last 24 h of rented bikes and timestamp information (hour, day, month). In contrast, other prediction models utilized only the last 24 h of rented bikes as input. Table 8 illustrates the averages of RMSE, MAE, and ${\mathrm{R}}^2$ from forecasting models, focusing on predicting bike-sharing demand (rented bikes) for a prediction horizon of 1 h on the Seoul bike-sharing 2018 dataset. Notably, our proposed MLP model exhibited superior performance, yielding the lowest RMSE and MAE, as well as the highest ${\mathrm{R}}^2$ in the testing sets when compared to alternative prediction models, including SVR, KNN, DT, AdaBoost, RF, and Linear Regression.

Figure 10 depicts the forecasted results produced by the proposed MLP model for a prediction horizon (PH) of 1 h. The graphical comparison includes the initial 100 points from both the testing sets and the prediction output. In terms of short-term prediction, the findings indicate that the proposed MLP model, leveraging the information from the past 24 h of rented bikes and timestamp data, can accurately predict the number of rented bikes for the following hour. The corresponding performance metrics for the proposed model are an RMSE of 135.453, an MAE of 91.522, and an ${\mathrm{R}}^2$ value of 0.952.

We conducted a comparative analysis with a prior study on a bike-sharing demand prediction model utilizing a Gradient Boosting model (Sathishkumar, Park & Cho, 2020a). The Gradient Boosting model yielded RMSE, MAE, and ${\mathrm{R}}^2$ values of 174.680, 109.890, and 0.920 in the testing sets, respectively. The results suggested that incorporating weather information enhanced the prediction accuracy of the models, with temperature and hour being identified as the most influential variables for hourly rented bike predictions. Unlike their approach of using current weather and date information as input features and rented bikes in each hour as the output value, our strategy involved using the previous values of rented bikes and timestamp information as input features to predict the next hour(s) of rented bikes. The RMSE, MAE, and ${\mathrm{R}}^2$ values generated by our proposed MLP model in the testing sets are 135.453, 91.522, and 0.952, respectively. Our experiment demonstrated that by leveraging the last 24 h of rented bikes and timestamp information, the prediction of the next hour’s rented bikes can be more accurate compared to the previous study.

In our study, we exclusively utilized the most recent data on rented bikes and timestamp information as input features for machine learning models. In a real-world scenario, incorporating other inputs such as current weather conditions might be challenging, as this information must be sourced from external databases. Therefore, the objective of our work is to achieve accurate prediction outputs from models using readily available features (from bike-sharing reservation data), such as the last known information of rented bikes and timestamp information, aiming to enhance the rebalancing system for bike-sharing companies.

In addition, we compared the proposed model with approaches employed in previous studies (Boonjubut & Hasegawa, 2022; Gammelli et al., 2022; Mehdizadeh Dastjerdi & Morency, 2022; Yang et al., 2020) that utilized RNN, LSTM, and GRU. Table 9 summarizes the performance across different prediction horizons (1, 3, 6, 12, and 24 h), where the number of rented bikes in the previous hour was used as the primary input for future demand forecasting. The results clearly indicate that the proposed MLP consistently achieves lower RMSE values and higher R² scores compared to the deep learning models.

The superior performance of the MLP can be attributed to the integration of temporal attributes in addition to lagged demand values. By incorporating features such as hour of the day, day of the week, and month of the year through a sliding window approach, the model effectively captures daily commuting cycles, weekly usage patterns, and seasonal variations. These findings demonstrate that with well-designed feature engineering, a relatively simple MLP can outperform more complex recurrent architectures. While LSTM and GRU are effective in modeling sequential dependencies, their performance may be limited when temporal context is not explicitly represented. By combining lag sequences with timestamp-based features, the proposed MLP achieves a strong balance between simplicity and predictive accuracy, highlighting its practical applicability for both short-term operational decisions and long-term strategic planning in bike-sharing systems.