Construction Project Cost Prediction Method Based on Improved BiLSTM

Abstract

In construction project management, accurate cost forecasting is critical for ensuring informed decision making. In this article, a construction cost prediction method based on an improved bidirectional long- and short-term memory (BiLSTM) network is proposed to address the high interactivity among construction cost data and difficulty in feature extraction. Firstly, the correlation between cost-influencing factors and the unilateral cost is calculated via grey correlation analysis to select the characteristic index. Secondly, a BiLSTM network is used to capture the temporal interactions in the cost data at a deep level, and the hybrid attention mechanism is incorporated to enhance the model’s feature extraction capability to comprehensively capture the interactions among the features in the cost data. Finally, a hyperparameter optimisation method based on the improved particle swarm optimisation algorithm is proposed using the prediction accuracy as the fitness function of the algorithm. The MAE, RMSE, MPE, MAPE, and coefficient of determination of the simulated prediction results of the proposed method on the dataset are 7.487, 8.936, 0.236, 0.393, and 0.996%, respectively, where MPE is a positive coefficient. This avoids the serious consequences of underestimating the cost. Compared with the unimproved BiLSTM, the MAE, RMSE, and MAPE are reduced by 15.271, 18.193, and 0.784%, respectively, which reflects the superiority and effectiveness of the method and can provide technical support for project cost estimation in the construction field.

1. Introduction

The scale and complexity of construction work continues to rise in the context of the global construction boom. Engineering cost is one of the most critical factors affecting the success of engineering projects; in the bidding stage for construction projects, accurate cost prediction is the basis for assessing the feasibility of the project and selecting design options, which directly affects how reasonable the bidding price is and the probability of success in winning the bid. In turn, these factors affect the economic feasibility and overall quality of the project. Consequently, an indispensable need exists to provide an accurate method of predicting construction cost.

In previous studies, linear models such as the autoregressive moving average (ARIMA) [1] model have been utilised due to their simple structure and robustness to data size and noise levels. However, the ARIMA model is unsuitable for capturing the nonlinearities of time series in engineering costs [2], which negatively affects prediction accuracy. To solve this problem, the support vector machine (SVM) [3,4,5], backpropagation (BP) neural network [6,7,8], and other machine learning models have been applied to cost prediction. Although such methods can effectively handle nonlinear problems, the SVM has limitations regarding data correlation processing and slow processing speed, and the BP neural network can quickly lose time series data and fall into local minimal values. With the continuous development of deep learning models, long short-term memory (LSTM) networks have achieved good prediction results in engineering cost prediction. Dong et al. [9] compared an LSTM network model to an SVM model through construction cost prediction experiments. The results revealed that the LSTM model can handle high-dimensional feature vectors and the selective recording of historical information better than an SVM while also having advantages in terms of prediction accuracy and parameter adjustment. Cao et al. [10] used LSTM to predict the construction cost of highway projects. The results showed that LSTM provides more accurate predictions in short, medium, and long-term prediction scenarios compared to the traditional model. However, LSTM still has shortcomings in dealing with the strong interactions between the features of the cost data. A bidirectional LSTM [11,12,13,14] (BiLSTM) network is a variant of the LSTM network that has an additional layer providing an inverse structure that can mine more data. Sima et al. [15] proved that the prediction accuracy of a BiLSTM network is higher than that of an LSTM network for time series problems, and suggested using BiLSTM instead of LSTM for issues related to time series analysis. A single model often fails to achieve optimal performance, and combined models can provide better performances and more accurate predictions by combining the strengths of multiple underlying models [16]. Chen et al. [17] incorporated an attention mechanism (AM) into a BiLSTM model for feature weight adjustment, which further improved the prediction effect of the model. Jalil et al. [18] used a particle swarm optimisation (PSO) algorithm to optimise the hyperparameters of a BiLSTM model and reduce its prediction error.

Implicit temporal patterns exist in construction cost data as factors such as materials, labor, and variable project duration. In addition, the data contain non-linear relationships and multi-level dependencies that allow for complex interactions between different cost characteristic indicators. However, existing studies have yet to consider these effects fully, so the constructed models are deficient in generalisation ability and prediction accuracy. Based on the above problems, this paper introduces the BiLSTM network into the field of engineering cost prediction. Grey correlation feature analysis (GRA) is performed to calculate the correlation between indicators and costs for feature indicator selection. The inter-temporal and inter-feature interactions in the cost data are captured by constructing the model prediction part by fusing the hybrid attention mechanism (HAM) and BiLSTM network. They aim to address the shortcomings of the traditional particle swarm optimisation algorithm. An improved particle swarm optimisation algorithm is proposed to find the optimal hyper-parameter configuration of the model, which avoids the tediousness and uncertainty problems of manual parameter adjustment. Experiments were conducted on 120 completed residential building projects in Jiangsu Province from 2019 to 2022, and the results show that the method has significant advantages in mean absolute error, root-mean-square error, and symmetric mean absolute percentage error.

2. Methodology for the Selection of Characteristic Indicators

Considering the specificity and diversity of construction projects, the factors affecting project costs are complex and diverse. However, not all factors have the same weight and importance regarding project costs [19]. Therefore, when selecting indicators, the principle of moderation should be considered to screen out indicators that can effectively describe project characteristics or have a significantly impactful project cost.

GRA is a data analysis method based on grey system theory that aims to study the interrelatedness between multiple indicators [20]. The data associated with each indicator are summed to determine the relative degree of each indicator, and the grey correlation between indicators is calculated to determine the degree of influence of each indicator on a problem. Compared to the traditional correlation analysis method, GRA can more accurately reflect the correlations between indicators and their degree of influence. Additionally, GRA has the advantages of a simple model, a small data volume, and interpretable results, which are suitable for selecting construction cost indicators. The steps involved in the grey correlation calculations are as follows.

• Raw data standardisation: the values of indicators in the raw data are mapped to intervals that ensure the same scale and range of variation among different variables.
• Absolute difference matrix construction: we take the absolute value of the difference between the target indicator and each candidate indicator to form an absolute difference matrix as:

  $${\mathrm{\Delta }}_i\left(k\right)=\left|x_0\left(k\right)-x_i\left(k\right)\right|$$

  (1)

  where $x_0$ and $x_i$ denote the standardised target and candidate indicators, respectively, $k$ denotes the kth sample, and $i$ denotes the ith candidate indicator.
• Correlation matrix construction: a correlation matrix is derived by considering each element in the absolute difference matrix to calculate the minimum and maximum values of the absolute difference matrix as:

  $${\zeta }_i\left(k\right)=\frac{{\mathrm{\Delta }}_{\min }+\rho {\mathrm{\Delta }}_{\max }}{{\mathrm{\Delta }}_i\left(k\right)+\rho {\mathrm{\Delta }}_{\max }}$$

  (2)

  where ${\mathrm{\Delta }}_{\min }$ and ${\mathrm{\Delta }}_{\max }$ denote the minimum and maximum values of the matrix, respectively, and $\rho$ denotes the resolution factor, which was set to 0.5 in this study.
• The grey correlation is calculated as:

  $$r_{0,i}=\frac{1}{n}{\displaystyle \sum }_{k=1}^n{\mathsf{\zeta }}_i\left(k\right)$$

  (3)

In this study, GRA was used to select characteristic indicators. By analysing the frequency of influencing factors and the difficulty of obtaining them in previous literature [21,22,23,24,25,26], we selected 18 indicators that have significant impacts on construction cost as candidate indicators: foundation type, structural type, fortification intensity, management level, façade material, aboveground floor area, underground floor area, number of aboveground floors, number of underground floors, aboveground floor height, underground floor height, interior wall decoration, types of doors and windows, number of elevators, roof type, concrete prices, reinforcing steel prices, and duration. The correlations between the candidate indicators and the unilateral cost were calculated via GRA, and the characteristic indicators were selected according to the resulting correlation value.

3. Predictive Modelling Design

The structure of the prediction model is shown in Figure 1, which consists of four parts: an input layer, a BiLSTM layer, an attention layer, and an output layer.

In this case, the input layer pre-processes the cost data to transform it into a shape that the BiLSTM can process. The BiLSTM layer is used to capture inter-temporal interactions in the costing data, and the attention layer is used to obtain inter-feature interactions in costing data. The output layer contains the dropout layer and the dense layer, where the dropout layer randomly switches off some neurons with probability p during training to set their outputs to zero, thus reducing the risk of overfitting. The weight matrix W of each neuron is multiplied by the probability p to maintain consistency during testing to compensate for the random switching-off operation during training, thus improving the model generalisation and robustness. In the dense layer, the neurons establish connections with each neuron in the previous layer, each connection has a corresponding weight, and each feature in the input data is passed through these connections to the dense layer where a weighted summation is performed in the neuron and nonlinear mapping is performed on the weighted output sequence to produce the output of the neuron.

3.1. Capturing Inter-Temporal Interactions

An LSTM network is a variant structure of a recurrent neural network (RNN) that solves the problem of gradient vanishing that occurs in RNNs when handling long sequential data. The structure of an LSTM network is illustrated in Figure 2.

An LSTM unit contains three gate structures called the forget gate, input gate, and output gate, respectively. The forget gate determines that information $f_t$ is overlooked based on the hidden state $h_{t-1}$ in the previous moment and the input $x_t$ in the current moment. The input gate takes the information update value $i_t$ and candidate cell state ${\widehat{c}}_t$ through $h_{t-1}$ and $x_t$. We multiply $f_t$ times the cell state $c_{t-1}$ in the previous moment and add the data from the cell update to obtain the new cell state $c_t$, which yields the initial output $o_t$ through $h_{t-1}$ and $x_t$ and the current moment output $h_t$ through $c_t$ and $o_t$. The corresponding formulas are defined as follows:

$$f_t=\sigma \left(w_f{h}_{t-1}+w_f{x}_t+b_f\right)$$

(4)

$$i_t=\sigma \left(w_i{h}_{t-1}+w_i{x}_t+b_i\right)$$

(5)

$${\hat{c}}_t=\tanh \left(w_c{h}_{t-1}+w_c{x}_t+b_c\right)$$

(6)

$$c_t=f_t{c}_{t-1}+i_t{\hat{c}}_t$$

(7)

$$o_t=\tanh \left(w_o{h}_{t-1}+w_o{x}_t+b_o\right)$$

(8)

$$h_t=o_t\tanh \left(c_t\right)$$

(9)

where $w_f$, $w_i$, $w_c$, and $w_o$ represent the weight matrices of the different gates; $b_f$, $b_i$, $b_c$, and $b_o$ represent the corresponding bias vectors; and $\sigma$ and $\tanh$ represent the activation functions.

In contrast to traditional LSTM, a BiLSTM network consists of two LSTM layers for the forward and backward directions and can fully account for both historical and future information. The structure of BiLSTM is illustrated in Figure 3.

The BiLSTM layer performs inter-temporal interaction capture by conducting both forward and backward pass steps. During the forward pass, the data in each time step are processed stepwise, starting from the first time step of the time series data. For each time step, the BiLSTM calculates and updates the hidden state $h_i^f$ of the current time step based on the input data $x_i$ of the current time step and hidden state $h_{i-1}^f$ from the previous time step. In contrast to the forward pass, in the backward pass, the BiLSTM begins processing forward from the last time step in the time series data. For each time step, BiLSTM calculates and updates the hidden state $h_i^b$ of the current time step based on the input data $x_i$ of the current time step and the hidden state $h_{i+1}^b$ from the subsequent time step. Because the information transmitted in the past and the information transmitted in the future regarding the cost data have different degrees of importance, a more comprehensive information representation $o_t$ is obtained by connecting $h_i^f$ and $h_i^b$ through the adaptive assignment of weights to capture the inter-temporal interactions at a deep level, which is calculated as follows:

$$h_t^f=\mathrm{LSTM}\left(x_t,h_{t-1}^f\right)$$

(10)

$$h_t^b=\mathrm{LSTM}\left(x_t,h_{t-1}^b\right)$$

(11)

$$o_t=w_f{h}_t^f+w_b{h}_t^b+b$$

(12)

where $w_f$ and $w_b$ represent the weight matrices of the forward and backward LSTM layers, respectively, and $b$ represents the bias vector.

3.2. Capturing Inter-Feature Interactions

The attention mechanism is a powerful tool in deep learning, and its design is inspired by the information processing mechanism of the human brain [27,28]. The attention mechanism makes the deep learning model pay more attention to the critical parts of the input data while ignoring the unimportant parts. The core idea is to assign different weights dynamically to different inputs so that the model can effectively learn and utilise important information related to the task and adapt more flexibly to the complex relationships in the data, thus improving the performance and generalisation of the model.

Given the complex spatial-temporal interactions between features of cost data, this article proposes a hybrid attention mechanism (HAM) to capture the inter-feature interactions, aiming to focus on both the temporal and feature dimensions of the data and to improve the ability of the model to capture multivariate time-series information. HAM contains two modules: channel attention and spatial attention. The structure is shown in Figure 4.

Here, the channel attention module aggregates the information about the feature dimensions through the tie pooling and maximum pooling operations to obtain the deep feature representation of the feature dimensions, which is then mapped into the attention weights on the feature dimensions by the multilayer perceptron. The $sigmoid$ activation function is used to activate the outputs of the summation to obtain the weights MC between zero and one. Unlike channel attention, the spatial attention module targets the temporal dimension by concatenating the pooled outputs and inputting them into a convolutional layer to extract features with high attentional weights in the temporal dimension, the convolution operation is utilised to emphasise the significance of the different time steps, and the $sigmoid$ activation function is used to activate the outputs of the summed outputs in order to obtain weights between zero and one MS. Connecting the channel-attention module and the spatial-attention module in a parallel way makes the modules consider the features of both the channel and spatial dimensions, avoiding the mutual influence of the channel and spatial dimensions and enabling the model to understand the multidimensional time series data better. This hybrid model helps the model capture key information about the data more comprehensively, which in turn improves the performance and generalisation of the model, as calculated in the following formula:

$$Mc\left(F\right)=\sigma \left(MLP\left(AvePool\left(F\right)\right)+MLP\left(MaxPool\left(F\right)\right)\right)$$

(13)

$$Ms\left(F\right)=\sigma \left(Conv\left(Conn\left[AvePool\left(F\right),MaxPool\left(F\right)\right]\right)\right)$$

(14)

$$F^{\prime }=Mc\left(F\right)\otimes F\otimes Ms\left(F\right)$$

(15)

where $\sigma$ represents the $sigmoid$ function, $MLP$ is a multilayer perceptual machine, $AvePool$ and $MaxPool$ denote tie pooling and maximum pooling, $Conv$ stands for convolution, $Conn$ stands for connection, and $Mc$ and $Ms$ represent the outputs of channel attention and spatial attention, respectively.

3.3. Loss Function Selection

During the training of the prediction model, the Adam optimisation algorithm was selected to update the model parameters based on the losses. Adam is a gradient descent-based optimisation algorithm that combines the advantages of momentum gradient descent and the root-mean-squared transfer algorithm [29]. It can effectively improve the training of a model, converge to optimal solutions faster, and iteratively update the weights and biases of a neural network based on training data to optimise the loss function output value. The loss function of the model was the MSE, which is calculated as follows:

$$MSE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-{\overline{y}}_i\right)}^2$$

(16)

where $n$ denotes the number of samples, $y_i$ denotes the actual value, and ${\overline{y}}_i$ denotes the model output value.

4. Hyperparametric Optimisation Methods

In the formulation of the cost prediction task, the accuracy of the HAM-BiLSTM model not only depends on the comprehensiveness of feature extraction, but is also affected by the combination of hyperparameters used by the model. In this study, the improved particle swarm optimisation (IPSO) algorithm was used to search for optimal hyperparameter combinations automatically to avoid the tedium and instability of manual parameter tuning and further improve model performance.

4.1. PSO Algorithm

The PSO algorithm [30,31,32] is an optimisation algorithm based on group intelligence inspired by the collective behaviors of flocks of birds or schools of fish. In the PSO algorithm, individuals are called particles. Each particle represents a potential solution and particles find the optimal solution by adjusting their speed and position. The inertia weight is an important parameter for balancing and adjusting the global and local search ability of the algorithm. The particle trajectory depends on the individual and social experiences of the particle, which are controlled by a learning factor. The velocity and position update formulas for a particle are defined as follows:

$$V_{i,t+1}=w{V}_{i,t}+c_1{r}_1\left(pbes{t}_i-X_{i,t}\right)+c_2{r}_2\left(gbes{t}_t-X_{i,t}\right)$$

(17)

$$X_{i,t+1}=X_{i,t}+\lambda V_{i,t+1}$$

(18)

where $V_{i,t}$ denotes the velocity of particle $i$ after $t$ iterations, $w$ denotes the inertia weight, $c_1$ and $c_2$ denote the individual learning factor and social learning factor, $r_1$ and $r_2$ denote random numbers, $x_{i,t}$ denotes the position of particle $i$ after $t$ iterations, and $\lambda$ denotes the velocity coefficient.

4.2. Phase Adjustment of Inertia Weights and Learning Factors

In traditional PSO algorithms, inertia weights and learning factors are typically set to fixed values, which limits global optimisation capability and convergence speed. Reference [33] dynamically adjusted inertia weights and learning factors according to changes in fitness, [34] dynamically adjusted inertia weights according to the number of iterations, and [35] dynamically adjusted the learning factor according to the number of iterations. These schemes all resulted in significant improvements in algorithm performance. However, considering the different states of particles in different periods, it is difficult to maximise PSO performance based on a single factor. Therefore, we combine the number of iterations and degree of adaptability and propose adjusting the inertia weights and learning factors in stages according to the following formulas:

$$w=w_{min}+\left(w_{max}-w_{min}\right)*{\left(\frac{f_{i,t}-f_{best}}{f_{i,t}}\right)}^m$$

(19)

$$c_1=c_{1,min}+\left(c_{1,max}-c_{1,min}\right)*{\left(\frac{f_{i,t}-f_{best}}{f_{i,t}}\right)}^m$$

(20)

$$c_2=c_{1,min}+\left(c_{1,max}-c_1\right)$$

(21)

$$m=\frac{t^2}{T^2}$$

(22)

where $w_{max}$ and $w_{min}$ are the maximum and minimum values of the inertia weight $w$, respectively; $f_{i,t}$ and $f_{best}$ are the fitness value of particle $i$ and the value of globally optimal particle fitness in iteration $t$, respectively; $c_{1,max}$ and $c_{1,min}$ denote the maximum and minimum values of the individual learning factor $c_1$, respectively; $c_2$ denotes the social learning factor; and $t$ and $T$ denote the current iteration and a maximum number of iterations, respectively. Additionally, $w\in \left[0.4,0.9\right]$ and $c_1\in \left[1,4\right]$.

Stage regulation divides the iterative process of the algorithm into three phases.

• Exploration phase: at the beginning of the iteration, when the value of $m$ is small, the inertia weights and learning factors are relatively less affected by adaptation, which reduces the magnitude of changes in the weight and learning factors, thereby ensuring that the particles have larger inertia weights and stronger individual learning abilities, promoting global exploration and escape from local optima.
• Equilibrium phase: in the middle of the iteration, when the value of $m$ is moderate, the inertia weights and learning factors become smoother through adaptation, which improves the social learning ability of the particles and helps them better utilise collective experiences for refined searching and adjustment.
• Convergence phase: in later iterations, when the value of $m$ is large, the inertia weights and learning factors are more sensitive to the influence of adaptation, which enhances the social cognitive ability of the particles, accelerates the convergence of the algorithm and helps guide particles to converge to the globally optimal solution more quickly.

As the value of $m$ varies in the interval $\left[0,1\right]$, it enables constraining the particles while adjusting their behavior, enabling them to strike a balance between exploration and convergence and thereby optimising the performance of the algorithm.

4.3. Negative Selection Evolutionary Strategy

In the iterative process of the algorithm, when the fitness values of the particles gradually become close to each other, the algorithm may be in a state of convergence or stagnation. At this time, the particles will be classified into advantageous and disadvantageous particles according to the average fitness value, and the disadvantageous particles generally converge to advantageous particles as the iteration progresses in the traditional PSO, which prevents the algorithm from finding a better solution in the solution space.

Inspired by biological evolution, each particle is considered as an individual and each dimension of a particle is considered as a gene segment. Under the premise of retaining the superior particles in the population, a gene crossover operation is used to break the gene structures of inferior individuals to promote their evolution. First, the ratio $C$ of the average fitness of each particle to the fitness of the global optimal particle is calculated using the following formulas:

$$\frac{f_{best}}{f_{ave}}=C$$

(23)

$$f_{ave}=\frac{{{\displaystyle \sum }}_{i=1}^n{f}_{i,t}}{n}$$

(24)

where $f_{best}$, $f_{ave}$, and $f_{i,t}$ are the average fitness of each particle, global optimal fitness, and fitness of particle $i$ after $t$ iterations, respectively, and $n$ is the number of particles.

Second, a critical point $K$ is defined by observing the magnitude of the change in $C$. When $C$ is greater than $K$, inferior particles are selected as bi-parental samples (e.g., particles A and B). Eventually, $x$ gene segments in the same position are randomly selected from the bi-parental samples to be interchanged, where the value of $x$ is in $\left[1,n-1\right]$, to produce new particles A1 and B1. This maintains the diversity of the population and prompts the algorithm to continue exploring the solution space. The negative selection evolutionary strategy is illustrated in Figure 5.

4.4. Hyperparameter Optimisation Process

The hyperparameter optimisation process is shown in Figure 6.

The specific steps are as follows:

• The batch size, number of LSTM layer cells, dropout probability, number of dense layer cells, and number of iterations in the HAM-BiLSTM model are selected as hyper-parameters to be optimised, and the respective search ranges are set.
• We defined the population size $N$, maximum number of iterations $T$, inertia weight $w$ search range, learning factor $c_1$ search range, and particle dimension $d$ of the IPSO algorithm and initialised the particle velocity $V$ and particle position $X$.
• The MSE of the true and predicted values was used as the fitness function.
• The performance was evaluated by training the HAM-BiLSTM model, calculating the fitness of each particle, and recording the average fitness $f_{ave}$ and global optimal fitness $f_{best}$ of the current particle population.
• We determined whether $f_{best}$ was less than the predefined value $g_{best}$ or whether the current number of iterations $t$ reached $T$. If so, then we output the optimal solution.
• The individual and population extremes were updated based on particle fitness.
• We updated $w$, $c_1$, and $c_2$ according to Equations (18)–(21).
• We determined whether the ratio $C$ of $f_{best}$ to $f_{ave}$ was greater than the critical value $K$ according to Equations (22) and (23). If so, the negative selection evolution strategy was executed.
• The $V$ and $X$ values of the current particle were updated according to Equations (7) and (8) before returning to step four.

5. Forecasting

The prediction process is shown in Figure 7.

The model prediction steps are as follows.

• Data pre-processing: the sample data were divided into training and testing sets and normalised.
• Model training: the model was trained using the training set, and the IPSO algorithm was utilised to determine the optimal combination of hyperparameters for the model.
• Model prediction: the trained model performed predictions on the testing set by employing an optimal combination of hyperparameters.
• Output results: the predictions of the model were inversely normalised to the original data range and outputted.

6. Experimental Results and Analysis

6.1. Data Preprocessing

This study is based on the data collection of residential cost data provided by Guanglianda Indicator Network and material price data provided by the Guangcai Network. The balance of the sample needs to be considered in the sample collection process, including the fact that different areas will have different land costs and code standards, high-rise buildings will require more complex and robust structural design and engineering than mid-rise buildings, and different forms of delivery represent different levels of renovation, all of which will affect the balance of the sample. Therefore, we restricted the region to Jiangsu Province during the sample collection process, the number of building floors to mid-rise residential buildings (18 floors and below), and the form of residential delivery to be simple. We initially screened 167 pieces of data from 2019 to 2023, excluding invalid samples with missing data, and ultimately obtained 156 valid sample data.

6.1.1. Relevance Analysis

The correlation values of the candidate indicators are listed in

Table 1. Candidate indicator correlation values.

Indicator Name	Relatedness	Indicator Name	Relatedness
Foundation type	0.9852	Aboveground floor height	0.9732
Structure type	0.9897	Underground floor height	0.9436
Fortification intensity	0.9603	Interior wall decoration	0.8763
Management level	0.9718	Types of doors and windows	0.8525
Facade material	0.8956	Number of elevators	0.8643
Aboveground floor area	0.9496	Roof type	0.8726
Underground floor area	0.8831	Concrete prices	0.9558
Number of aboveground floors	0.9764	Reinforcing steel prices	0.9623
Number of underground floors	0.9668	Duration	0.9311

, and the candidate indicators with correlation values greater than 0.9 are designated as feature indicators. The foundation type ($x_1$), structure type ($x_2$), fortification intensity ($x_3$), management level ($x_4$), aboveground floor area ($x_5$), number of aboveground floors ($x_6$), number of underground floors ($x_7$), aboveground floor height ($x_8$), basement floor height ($x_9$), concrete prices ($x_{10}$), reinforcing steel prices ($x_{11}$), and duration ($x_{12}$) indicators were selected as the model inputs, and the unit cost ($y$) was defined as the model output.

6.1.2. Data Quantification

For qualitative indicators to be used as inputs for the model, they must be quantified. These indicators include foundation type, structure type, fortification intensity, and management level. The detailed quantification process is defined in

Table 2. Quantitative table of qualitative indicators.

Foundation Type	Value	Structure Type	Value	Fortification Intensity	Value	Management Level	Value
Raft slab foundation	1	Shear wall structure	1	6 degrees	1	Excellent	1
Pile foundation	2	Framework structure	2	7 degrees	2	Good	2
Mantenna foundation	3	Frame shear construction	3	8 degrees	3	Bad	3
				9 degrees	4

.

6.1.3. Data Normalisation

The sample data were divided at a ratio of 8:2 to obtain training and test sets. The input quantities of the model were the 12 feature indicators selected above, and the output quantities were the unilateral costs. Considering the magnitude differences between different feature indicators, to eliminate the impacts of such differences on the model training, reduce the training time, and accelerate the convergence speed, we adopted the max-min method to normalise all sample data. The max-min method eliminates the impact caused by the large magnitude differences by mapping original data to the interval $\left[0,1\right]$ such that different feature indicators have the same scale range. This process is formulated as follows:

$$y_i=\frac{x_i-x_{min}}{x_{max}-x_{min}}\left(i=1,2,\cdots ,n\right)$$

(25)

where $x_{min}$ denotes the minimum value of indicator $x$ in the data and $x_{max}$ denotes the maximum value of indicator $x$ in the data.

Some of the data obtained after normalising the sample data are presented in

Table 3. Normalised data.

Number	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	y
1	1.00	0.00	1.00	1.00	0.54	0.91	1.00	0.00	0.78	0.34	0.47	0.98	0.83
2	1.00	0.00	1.00	1.00	1.00	0.86	1.00	0.00	0.69	0.34	0.47	0.37	1.00
3	1.00	0.50	0.50	1.00	0.20	0.96	0.67	0.00	0.89	0.36	0.47	0.22	0.63
4	0.00	0.00	1.00	0.50	0.18	0.91	0.67	1.00	0.78	0.36	0.47	0.61	0.18
⁝
155	1.00	0.50	0.50	0.50	0.02	1.00	0.18	0.33	0.52	0.94	0.36	0.27	0.23
156	1.00	1.00	0.50	0.00	0.06	1.00	0.27	0.33	0.71	0.04	0.71	0.86	0.59

.

6.2. IPSO Algorithm Performance Test

To verify the superiority of the IPSO algorithm, the PSO and IPSO algorithms were executed on eight benchmark test functions. The population size for the PSO and IPSO algorithms was set to 50, the number of iterations was set to 100, the dimensionality was set to 30, and the number of independent runs was set to 20. As shown in

Table 4. Comparison of PSO and IPSO algorithm test results.

Function	Arithmetic	Min	Max	Average	STD
Sphere	PSO	2.23 × 10−1	1.12	5.53 × 10−1	2.24 × 10−1
	IPSO	3.62 × 10−30	3.32 × 10−26	1.46 × 10−26	1.70 × 10−27
Schwefel	PSO	1.22	4.05	1.97	6.30 × 10−1
	IPSO	1.28 × 10−15	3.61 × 10−13	3.64 × 10−14	9.27 × 10−13
Sum squares	PSO	5.67	2.20 × 101	1.16 × 101	5.06
	IPSO	6.48 × 10−30	5.78 × 10−26	2.89 × 10−26	3.97 × 10−27
Rosenbrock	PSO	6.18 × 101	2.68 × 102	1.26 × 102	5.16 × 101
	IPSO	2.87 × 101	2.90 × 101	2.89 × 101	3.76 × 10−2
Rastigrin	PSO	4.73 × 101	1.01 × 102	7.36 × 101	1.45 × 101
	IPSO	1.82 × 10−1	1.35 × 101	6.46	3.79
Ackley	PSO	7.91 × 10−1	3.18	2.28	4.75 × 10−1
	IPSO	8.88 × 10−16	8.88 × 10−16	8.88 × 10−16	0
Griewank	PSO	2.43 × 10−2	1.04 × 10−1	4.86 × 10−2	1.68 × 10−2
	IPSO	2.84 × 10−4	7.60 × 10−3	8.78 × 10−4	1.60 × 10−3
Penalised	PSO	1.71 × 10−2	2.45 × 10−1	9.51 × 10−2	6.30 × 10−2
	IPSO	2.79 × 10−5	2.62 × 10−4	1.53 × 10−4	5.61 × 10−5

, the minima, maxima, averages, and standard deviations obtained by the IPSO algorithm were better than those obtained by the PSO algorithm and were closer to the global optimum. Figure 8 compares the convergence curves of the PSO and IPSO algorithms, and one can draw the following conclusions: the IPSO algorithm performs better in terms of convergence accuracy, and it enables the particles to gather more stably in the vicinity of the global optimum, therefore improving the ability of the PSO algorithm to overcome local extreme values.

6.3. Hyperparametric Optimisation Results

The hyperparameters to be optimised in the model are the package batch size $N_{batch}$, number of LSTM layer cells $N_{LSTM}$, dropout probability $N_{Dropout}$, number of dense layer cells $N_{Dense}$, and number of iterations $N_{epoch}$. In this paper, the hyperparameter search ranges were set as shown in

Table 5. Hyperparameter search ranges.

Hyperparameter	Search Scope
$N_{batch}$	1–30
$N_{LSTM}$	1–100
$N_{Dropout}$	0.1–0.5
$N_{Dense}$	1–100
$N_{epoch}$	200

.

The results of IPSO and traditional PSO optimisation searches were compared. The number of iterations of the optimisation algorithm was set to 50, and the population was set to 30. The variation of fitness values during the iteration process of the PSO and IPSO algorithms is shown in Figure 9.

One can see that the IPSO algorithm has a better ability to find the optimum compared to the PSO algorithm. The optimisation results are summarised in

Table 6. Optimisation results.

Algorithm	${\mathit{N}}_{\mathit{b}\mathit{a}\mathit{t}\mathit{c}\mathit{h}}$	${\mathit{N}}_{\mathit{L}\mathit{S}\mathit{T}\mathit{M}}$	${\mathit{N}}_{\mathit{D}\mathit{r}\mathit{o}\mathit{p}\mathit{o}\mathit{u}\mathit{t}}$	${\mathit{N}}_{\mathit{D}\mathit{e}\mathit{n}\mathit{s}\mathit{e}}$	${\mathit{N}}_{\mathit{e}\mathit{p}\mathit{o}\mathit{c}\mathit{h}}$
PSO	5	26	0.2	22	200
IPSO	3	32	0.1	18	200

.

6.4. Evaluation Indicators

In this study, we used five evaluation indices to evaluate the model: the mean absolute error (MAE), root-mean-square error (RMSE), mean percentage error (MPE), mean absolute percentage error (MAPE), and coefficient of determination R². The smaller the values of MAE, RMSE, and SMAPE, the better the performance of the model, and the larger the value of R², the better the model fitting ability. In assessing the feasibility of a project, the consequences of underestimating costs are more serious than those of overestimating them. The MPE was used to measure the directionality of the forecasting index. The specific formulas are as follows:

$$MAE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|y_i-{\widehat{y}}_i\right|$$

(26)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(27)

$$MPE=\frac{100\%}{n}{\displaystyle \sum }_{i=1}^n\frac{{\widehat{y}}_i-y_i}{y_i}$$

(28)

$$MAPE=\frac{100\%}{n}{\displaystyle \sum }_{i=1}^n\left|\frac{y_i-{\widehat{y}}_i}{y_i}\right|$$

(29)

$$R^2=1-{\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}/{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(30)

where $y_i$ is the true value, ${\widehat{y}}_i$ is the predicted value, ${\overline{y}}_i$ is the mean of the true values, and $n$ is the number of data used.

6.5. Comparative Experiments and Analyses

Figure 10 compares the predicted and real value curves of the proposed model on the test set for the unilateral cost of construction projects. The predicted and real value curves fit well together, which shows that the model has good prediction performance for the unilateral cost.

To verify the superiority of BiLSTM in capturing the interactions between data tensors, BiLSTM is compared with the traditional benchmark models BP, SVM, and LSTM. The prediction results and errors are shown in Figure 11. The prediction results demonstrate that BiLSTM has the best fit to the true values, and the overall error variation of BiLSTM is relatively smooth, as seen from the error curves.

Table 7. Comparison of single model performance evaluation.

Model	MAE	RMSE	MPE/%	MAPE/%	R2
BP	33.980	39.174	−0.124	1.787	0.921
SVM	29.254	36.068	−0.154	1.517	0.933
LSTM	25.294	30.236	0.472	1.329	0.953
BiLSTM	22.758	27.129	0.054	1.187	0.962

shows the results of the comparison of each single model on the five evaluation indicators. For the construction cost data, the average errors of LSTM and BiLSTM with time series capturing capability are relatively small, and the MAE, RMSE, and MAPE of BiLSTM are reduced by 2.536, 3.107, and 0.142%, respectively, compared with those of LSTM, which verifies the superiority of the BiLSTM bi-directional processing mechanism in capturing the inter-temporal interactions of the data.

To validate the effectiveness of the model components further, the model was compared with BiLSTM, AM-BiLSTM, HAM-BiLSTM, and PSO-HAM-BiLSTM. The prediction results and error comparisons are shown in Figure 12.

As can be seen from the figure, although the prediction results of each model have trends similar to that of the true value curve, the error curves of each model have a large difference. Comparing the error curves of HAM-BiLSTM and AM-BiLSTM reveals that the number of errors greater than 20 decreases from 13 to 8, reflecting that HAM enhances the robustness of the model in terms of its ability to capture inter-feature interactions due to traditional AM. The HAM-BiLSTM with hyperparameters optimised by IPSO has a smoother trend of error curves and a smaller range of error fluctuations relative to the HAM-BiLSTM with hyperparameters optimised by the conventional PSO, reflecting the fact that IPSO can further improve the accuracy of the model.

Table 8. Comparison of performance evaluation of combined models.

Model	MAE	RMSE	MPE/%	MAPE/%	R2
BiLSTM	22.758	27.129	0.054	1.181	0.962
AM-BiLSTM	20.299	23.042	0.838	1.093	0.973
HAM-BiLSTM	16.047	19.372	0.647	0.856	0.981
PSO-HAM-BiLSTM	13.723	16.594	0.161	0.717	0.986
IPSO-HAM-BiLSTM	7.487	8.936	0.236	0.393	0.996

shows the comparison results of the models on the five evaluation metrics, the mean percentage error of the five models is positive, avoiding the serious consequences of underestimating the cost, and all the other error metrics of HAM-BiLSTM after optimising the hyperparameters by IPSO are significantly reduced. In particular, the average absolute error of HAM-BiLSTM is reduced by 6.711 and 4.252%, the RMSE is reduced by 7.757 and 3.670%, and the average absolute percentage error is reduced by 0.325 and 0.237%, respectively. This further demonstrates the effectiveness of HAM in capturing the interactions between the features and enhancing the feature extraction capability of the model. The average absolute error of HAM-BiLSTM after IPSO optimisation of hyperparameters is reduced by 8.560 and 6.236%, the RMSE is reduced by 10.436 and 7.658%, and the average absolute percentage error is reduced by 0.463 and 0.324%, in comparison with those of HAM-BiLSTM and AM-BiLSTM after PSO of the hyperparameters, indicating that the improved PSO algorithm has a better search capability, which can better overcome the problem of falling into local optimal solutions and then find a more comprehensive search of the hyperparameter space to find a better combination. Thus, the improved particle swarm optimisation algorithm has stronger search capability and can better overcome the problem of falling into local optima, thus searching the hyperparameter space more comprehensively and finding better hyperparameter combinations.

7. Conclusions

In this study, a new construction cost prediction method was developed to address the problem of strong interaction among construction cost data. Based on the exploration and improvement of the BiLSTM network in the field of deep learning, the temporal interactions in the data are captured at a deep level, and the HAM module is added to the network, which enhances the ability of the network to capture the interactions among the features and then constructs an effective construction project cost prediction model. In addition, to select the feature indicators that have greater impacts on the cost as model inputs, GRA was used to assess the importance of the factors related to influencing the cost, and the most valuable feature indicators were screened according to the principle of proportionality. To avoid the tediousness and uncertainty of manual parameter tuning, an improved particle swarm optimisation algorithm was designed to find the optimal hyperparameter combination for the model automatically. Through comparative experiments and analyses, the performance of the method in construction cost prediction work was verified.

Even though the results obtained with the improved model are very positive compared to the BiLSTM network, it must be recognised that it has limitations. Firstly, since HAM contains two attention modules, it embodies a powerful feature extraction capability, which will pay more computational resources accordingly. Secondly, IPSO is committed to improving the algorithm’s ability to find the global optimal solution, and we do not reduce the number of iterations and running time of the algorithm relative to the original algorithm.

In future research, we will devote ourselves to overcoming these limitations by focusing on an improved method to simplify the model structure in the feature extraction problem, and a hyperparameter optimisation method that combines computing accuracy and computing speed in the hyperparameter optimisation problem. In addition, for the construction cost prediction problem, not only is an accurate prediction model needed, but also the selection of characteristic indexes is equally important. We will analyse the influencing factors of the cost in depth to contribute more value to the development of this field.