Simulation-Optimization for the Planning of Off-Site Construction Projects: A Comparative Study of Recent Swarm Intelligence Metaheuristics

Off-site construction is a modern construction method that brings many sustainability merits to the built environment. However, the sub-optimal planning decisions (e.g., resource allocation, logistics and overtime planning decisions) of off-site construction projects can easily wipe away their sustainability merits. Therefore, simulation modelling—an efficient tool to consider the complexity and uncertainty of these projects—is integrated with metaheuristics, developing a simulation-optimization model to find the best possible planning decisions. Recent swarm intelligence metaheuristics have been used to solve various complex optimization problems. However, their potential for solving the simulation-optimization problems of construction projects has not been investigated. This research contributes by investigating the status-quo of simulation-optimization models in the construction field and comparing the performance of five recent swarm intelligence metaheuristics to solve the stochastic time–cost trade-off problem with the aid of parallel computing and a variance reduction technique to reduce the computation time. These five metaheuristics include the firefly algorithm, grey wolf optimization, the whale optimization algorithm, the salp swarm algorithm, and one improved version of the well-known bat algorithm. The literature analysis of the simulation-optimization models in the construction field shows that: (1) discrete-event simulation is the most-used simulation method in these models, (2) most studies applied genetic algorithms, and (3) very few studies used computation time reduction techniques, although the simulation-optimization models are computationally expensive. The five selected swarm intelligence metaheuristics were applied to a case study of a bridge deck construction project using the off-site construction method. The results further show that grey wolf optimization and the improved bat algorithm are superior to the firefly, whale optimization, and salp swarm algorithms in terms of the obtained solutions’ quality and convergence behaviour. Finally, the use of parallel computing and a variance reduction technique reduces the average computation time of the simulation-optimization models by about 87.0%. This study is a step towards the optimum planning of off-site construction projects in order to maintain their sustainability advantages.


Introduction
Off-site construction (OSC) is a sustainable construction method in which parts of the structure are produced off-site and then transported to the construction site for erection [1,2]. OSC has many advantages for the improvement of productivity [3], quality [4], and safety [5] of the built environment. Besides this, OSC brings environmental improvements to the construction industry [6]. OSC could reduce greenhouse gas (GHG) emissions by 48% and 43% during the operation and construction phases, respectively [7]. Furthermore, OSC could reduce construction waste compared with traditional construction [8]. Consequently, OSC is adopted in both building and infrastructure projects. For instance, OSC is an efficient way to narrow the gap between the supply and demand of both public and private buildings [9]. In addition, many transportation agencies prefer OSC to deliver infrastructure projects to mitigate the traffic interruption resulting from the traditional cast-in-situ method [10].
Despite the benefits of OSC, its complexity and fragmentation require extensive planning from the project managers [1,11]. The need for extensive planning is one of the key barriers that hinder the wider adoption of OSC [12]. For successful project delivery, the project manager has to make critical decisions while considering the dynamics and uncertainty of OSC projects [13]. Some of these decisions are related to resource planning, such as the number of required crews and the on-site and off-site equipment [14]. The rest are related to logistics aspects, such as the location of the storage yards and their capacities. These decisions should be made to deliver the project on time and within the allowable budget. This problem is called the time-cost trade-off problem (TCTP), which has been extensively addressed in the literature considering its different variants and has been solved using both exact and approximate optimization methods. Approximate methods are preferable because they efficiently handle complex projects with large networks in a reasonable computation time [15]. Metaheuristics were widely used to solve single and multi-objective TCTP. The Genetic Algorithm (GA) [16], Ant Colony Optimization (ACO) [17], and Particle Swarm Optimization (PSO) [18] were applied to solve singleobjective TCTP. For multi-objective TCTP, researchers used dynamic programming [19], multi-objective ACO [20], non-dominated sorting genetic algorithm II (NSGA-II), multiobjective simulated annealing, and multi-objective PSO [21]. Despite the contribution of these studies, their common limitations are: (1) The proposed models addressed the traditional cast-in-situ construction method. Hence, they do not fit with OSC projects characterized by a three-stage supply chain (i.e., production, logistics, and installation stages) [22]. (2) The previous models failed to capture the complex interactions between the resources and activities of OSC projects. These complex interactions form queues inside the system, making the analytical methods impractical for the modelling of such queuing systems [23]. (3) the previous models assumed that the activities' durations and costs are deterministic without considering the uncertainty associated with these types of projects [24].
In order to address these gaps, researchers resorted to Discrete-Event Simulation (DES) to capture the complexity and uncertainty of OSC projects without the need to develop sophisticated mathematical models. However, DES is not an optimization tool on its own; it can solely enable project planners to study the effect of the proposed resource planning and logistics decisions on the project duration, cost and productivity [10]. In order to tackle this limitation, researchers tend to combine simulation and metaheuristics into one approach called "simulation-optimization (SO)". Swisher et al. [25] defined SO as a "structured approach to determine optimal input parameter values, where optimal is measured by a function of output variables-steady-state or transient-associated with a simulation model". Given this definition, the solution obtained from the SO approach is the values of the input parameters associated with the DES model. This solution's performance is evaluated by the DES model. The metaheuristic algorithm uses the current solution and its evaluation to find a new set of input values [26]. Therefore, different types of metaheuristics have been integrated with DES in the SO approach to solve multiple construction problems characterized by complexity and uncertainty.
Evolution-based metaheuristics are a subcategory of metaheuristics inspired by the laws of natural evolution. The advantage of these metaheuristics is that the new generation of solutions is created using the fittest solutions in the current generation, hence improving the solutions over the generations [27]. Genetic Algorithms (GAs), which are the most popular evolution-based metaheuristics, have frequently been used in SO for different construction applications such as sewer pipeline installation [28], earthmoving operations [29], the production planning of precast components [30], the planning of prefabricated building construction [31], and bridge construction planning [32].
Swarm Intelligence (SI) metaheuristics are another subcategory of metaheuristics. These metaheuristics are usually inspired by the social and hunting behaviours of colonies and swarms of animals in nature. According to Mirjalili et al. [33], SI metaheuristics offer multiple advantages over evolution-based metaheuristics. Firstly, they retain the information associated with the search space over iterations, while evolution-based metaheuristics dispose the information of previous generations once a new generation is created. Secondly, they are generally easier to implement because they have fewer operators than the evolution-based metaheuristics (i.e., operators for crossover, mutation, elitism, etc.). Despite these advantages, SI metaheuristics have rarely been used in SO for construction applications. Examples of these applications are a concrete placement operation using PSO [34], and the construction planning of bridge deck construction using ACO [35] and PSO [36]. Because of the introduction of PSO and ACO, SI metaheuristics have witnessed a dramatic development, leading to tens of metaheuristics that have potential over PSO and ACO, while most of them use a lower number of tuning parameters. Examples of such recent SI metaheuristics are the firefly algorithm (FA) [37], grey wolf optimization (GWO) [38], the whale optimization algorithm (WOA) [39] and the salp swarm algorithm (SSA) [40].
Given the above-mentioned preliminary review of the use of metaheuristics in SO in the construction domain, this study contributes to the literature by addressing three research questions that have not yet been addressed: (1) What is the status quo of the use of metaheuristics for SO in the construction research field? (2) What is the potential of the use of recent SI metaheuristics in SO to optimize the planning of OSC projects? (3) Which one among the selected SI metaheuristics has better performance in obtaining more optimal planning decisions for OSC? In order to address these questions, the metaheuristics used in the SO of construction applications are first reviewed using a systematic review. Secondly, a DES model of a bridge construction project using the OSC method is developed. Then, its optimization model is formulated by defining the decision variables, objective functions and constraints. Next, the developed model is integrated with the SI metaheuristics under study. In order to reduce the computation time of the SO runs, parallel computing and a variance reduction technique called common rand numbers (CRN) are used. After that, the developed models are applied to a case study. Finally, a comparative analysis between the SI metaheuristics under study is conducted using qualitative and quantitative measures. Following the recommendations of Sörensen [41] and Juan et al. [26] regarding the need to conduct comparative studies for SO, this study contributes by comparing the performance of recent SI metaheuristics in SO of OSC projects.
The rest of this study is organized as follows: SO applications in the construction field are reviewed in Section 2. Secondly, the problem description and the formulation of the optimization model are discussed in Section 3. Thirdly, the adopted SO approach, the DES modelling of an OSC project, and the five SI metaheuristics under study are illustrated and coded in Section 4, including their parameter tuning. In addition, Section 4 elaborates on the methods used to reduce the computation time, namely, the CRN method and parallel computing. After that, the developed simulation and optimization models are applied to a case study discussed in Section 5. Finally, the analysis and discussion of the performance of the SI metaheuristics are provided in Section 6 in terms of their convergence behaviour and statistical results, followed by conclusions in Section 7.

Literature Review of Simulation-Optimization (SO) in the Construction Research Field
This section answers the first research question mentioned in the previous section by reviewing the SO applications in the construction research field. The specific aims are the identification of: (1) studies that adopted SO in the construction domain, (2) which simulation methods (e.g., DES, agent-based simulation (ABS), system dynamics (SD)) have been used in these studies, (3) which metaheuristics have been applied in these studies, (4) the construction applications addressed, and (5) whether these studies applied parallel computing or VRTs. These objectives can be applied using a systematic review method [42,43]. This method consists of two main steps: the extraction of related studies and their analysis [44]. The first step starts by selecting a number of relevant keywords for the database search. The following search code was used in the Scopus database: (TITLE-ABS-KEY ("optimization") AND TITLE-ABS-KEY ("construction") OR TITLE-ABS-KEY ("infrastructure") OR TITLE-ABS-KEY ("building") AND TITLE-ABS-KEY ("simulation") AND TITLE-ABS-KEY ("DES") OR TITLE-ABS-KEY ("discrete") OR TITLE-ABS-KEY ("continuous") OR TITLE-ABS-KEY ("discrete-event") OR TITLE-ABS-KEY ("system dynamic") OR TITLE-ABS-KEY ("system dynamics") OR TITLE-ABS-KEY ("SD") OR TITLE-ABS-KEY ("agent") OR TITLE-ABS-KEY ("multi-agent") OR TITLE-ABS-KEY ("agent-based") OR TITLE-ABS-KEY ("ABS") OR TITLE-ABS-KEY ("ABM")) AND (LIMIT-TO (DOCTYPE, "ar")) AND (LIMIT-TO (LANGUAGE, "English")) AND (EXCLUDE (SUBJAREA, "MATE") OR EXCLUDE (SUBJAREA, "PHYS") OR EXCLUDE (SUBJAREA, "EART") OR EXCLUDE (SUBJAREA, "CHEM") OR EXCLUDE (SUBJAREA, "CENG") OR EXCLUDE (SUBJAREA, "BIOC") OR EXCLUDE (SUBJAREA, "MEDI") OR EXCLUDE (SUBJAREA, "PHAR") OR EXCLUDE (SUBJAREA, "AGRI") OR EXCLUDE (SUBJAREA, "NEUR") OR EXCLUDE (SUBJAREA, "IMMU") OR EXCLUDE (SUBJAREA, "HEAL") OR EXCLUDE (SUBJAREA, "ARTS") OR EXCLUDE (SUBJAREA, "NURS") OR EXCLUDE (SUBJAREA, "PSYC")). As of November 2021, using this code results in the identification of 745 studies, as shown in Figure 1. Figure 1 shows a common protocol used in the systematic review, called the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA), to screen the relevant studies [45]. After checking the title and abstract of the retrieved studies, 637 studies were found to be irrelevant to the construction research field, namely the implementation stage of construction projects. Then, the full texts of the remaining studies (i.e., 745 -637 = 108 studies) were evaluated, and 38 of them were found to be relevant. By checking the reference list and studies that cited each of the 38 studies, eight more relevant studies were detected, raising the number of related studies to 46, as shown in Figure 1.
The second step of the adopted systematic review method is to evaluate the 46 relevant studies to identify the simulation methods and metaheuristics used, their applications, and the use of parallel computing and VRTs. Figure 2 provides an overview of the classification of the 46 studies based on multiple criteria. Figure 2a shows that most of these studies (59%) are related to the traditional construction method, while Figure 2b shows that 65% of these studies addressed single optimization problems. Regarding the metaheuristics used in these studies, Figure 2c shows that the majority of these studies used evolutionary-based metaheuristics (i.e., GA and evolutionary algorithms (EA)), while one-fifth of these studies adopted SI metaheuristics (i.e., PSO and ACO). Table 1 provides details of this information, including the adoption of parallel computing and VRTs to reduce the computation time, as well as the type of simulation methods used in these studies. Most of these studies integrated DES with optimization, while very few studies conducted SO using ABS and SD, as shown in Table 1. Furthermore, Table 1 shows that three studies used parallel computing, while only one study adopted VRTs. Despite the vast improvements in computing power, computation time reduction is inevitable for five reasons: (1) Construction projects are usually characterized by uncertainty. For example, their activities' durations are best represented by probabilistic distributions; hence, multiple simulation replications are required to obtain reliable results, which in turn prolong the computation time. (2) SO with metaheuristics requires the call of the simulation model frequently to evaluate the fitness of each individual solution of the population in each iteration. (3) Finding the optimum values of the controllable parameters for each metaheuristic requires extensive numerical experiments. (4) Construction researchers tend to hybridize DES with other simulation approaches such as system dynamics [46] or integrate it with neural networks [47]. Such research developments would increase the capabilities of DES but at the expense of the computation time. (5) OSC projects are dynamic in nature, meaning that optimized decisions made in the early planning stage might be infeasible in later construction stages due to inevitable uncertainties [32]. Therefore, there is a need to make near-optimum planning decisions in a reasonable computation time.  Site layout planning and resource allocation Minimize the project's duration DES GA [84] Sim. M, simulation method; VRT, variance reduction technique; DES, discrete-event simulation; GA, genetic algorithm; PSO, particle swarm optimization; SA, simulated annealing; ACO, ant colony optimization; NSGA-II, non-dominated sorting genetic algorithm II; fmGA, fast messy genetic algorithm; EA, evolutionary algorithm ABS, Agent-based simulation; DE, differential evolution; ICA, imperialist competitive algorithm; MC, Monte Carlo simulation; SD, system dynamics. * Sequential simulation-based optimization. Given the above analysis of the literature of SO in the construction domain, a number of research gaps were identified. Firstly, most of the studies used traditional SO methods by solely integrating DES and metaheuristics without benefiting from recent approaches such as reinforcement learning [47], machine learning [85], metamodeling [86], and hybrid simulation [87]. Secondly, although solving stochastic SO models takes a long computation time, hindering its wide adoption [32], very few studies integrated SO models with computation-time-reduction techniques such as VRTs and parallel computing. Thirdly, SI metaheuristics have not been widely applied to the SO of construction applications despite their advantages over evolutionary-based metaheuristics, as discussed in the previous section. Moreover, among these SI metaheuristics, only PSO and ACO were adopted. Because of the development of PSO and ACO, many other SI metaheuristics were proposed and demonstrated to perform better in multiple optimization problems while using fewer parameters. Examples of these SI metaheuristics are FA, GWO, WOA and SSA, to name a few. Therefore, it is more worthy to test and compare such recent SI metaheuristics than to introduce another metaheuristic that reuses existing concepts from previous metaheuristics [41]. Therefore, many researchers have conducted comparison studies among multiple metaheuristics in different research fields in order to evaluate their performances in different optimization problems, such as deterministic construction TCTP [88], bridge deck repairs [89], the design of water distribution networks [90], concrete foundation design [91], the domain decomposition of finite element models [92], supply chain management [93], the production scheduling of precast components [61], aircraft routing problems [94], and the design of steel-frame structures [95,96]. Moreover, the SO of OSC represents a new test field to validate and compare the potential of these new SI metaheuristics to solve a new set of real-life optimization problems characterized by uncertainty and complexity [26]. Given these gaps, the contribution of this study is the comparison of recent SI metaheuristics that have not been investigated before solving a stochastic SO problem for OSC. This study could be seen as one of the early attempts to evaluate the performance of recent SI metaheuristics in the solution of this type of optimization problem characterized by uncertainty. Consequently, it could highlight some defects in these metaheuristics for their developers to tackle and could meanwhile suggest the most promising metaheuristics for a specific application in order for researchers and practitioners to further improve or exploit them.

Problem Statement and Formulation
Precast full-span construction using a launching gantry to construct bridge decks is considered an example of infrastructure OSC projects in this research. The sequence of its construction processes and the associated optimization problem are illustrated here.

Problem Description
OSC projects are characterized by their three-echelon supply chain, starting from the production stage at the casting yard, followed by the logistics and installation stages, as shown in Figure 3 [97]. The production stage starts by seizing a rebar cage mould, a reinforcement crew, and casting yard space to place the steel reinforcement and the stressing ducts of the bottom slab and the webs of the box girder. Then, the inner mould is loaded to the rebar mould by the preparation crew. After that, the reinforcement of the top slap is installed by the reinforcement crew. After finishing the rebar cage, it is lifted by the yard crane to an outer mould, where the casting crew pours the concrete. After the concrete pouring, the girders are cured using either traditional methods or accelerated methods. Next, the preparation crew removes the inner mould so that the pre-stressing crew can perform the first stage of post-tensioning. The second stage of the post-tensioning is conducted after moving the girder to the storage area, the capacity of which depends on the location of the casting yard. At this point, the production process of the girder is finished, and the installation process at the construction site can commence after the transportation of the girders to the bridge construction location by specific trailers. At the access point of the bridge construction site, the girder is unloaded from the trailer and loaded onto a trolley. Then, the trolley moves to the location of the launching gantry while the unloaded trailer returns to the storage area to transport the other girders. After the girder arrives at the launching location, the launching gantry repositions itself to the new span's location. Then, the launching gantry picks up the girder from the trolley and lays it on the bridge piers while the unloaded trolley returns to the access point of the bridge. Finally, the permanent bearings are grouted to carry the load of the girder after transferring its load from the temporary bearings. It is worth mentioning that the trolley-loading process is allowed only after the load transfer of the girders. The sequence of the aforementioned production, transportation and installation processes, along with their associated resources, is represented in Figure 3. The figure shows that the girders, in order to proceed from one process to the next one, need to wait until all of the resources required to accomplish the next process are available. The main objective is to find the optimum/near-optimum number of human and equipment resources required to accomplish the production, logistics and installation processes, as well as the number of daily working hours and weekly working days, and the location of the production facility (i.e., the casting yard) that minimize the project's duration and costs.

Optimization Model Formulation
This section elaborates on the development of the optimization model for the infrastructure OSC project, including the model parameters, decision variables, objective functions, and constraints.

Model Parameters and Decision Variables
The model indices, parameters and variables are summarized in Table 2. The model parameters include: (1) the duration of the production, logistics and installation processes of the OSC project; (2) the hourly cost of the production, logistics and installation equipment, and human resources; (3) the daily site indirect cost; (4) the mobilization and demobilization cost of the production, logistics and installation equipment, and human resources; and (5) the hourly storage cost.

Indicies
Prefabricated box girders; g ∈ G, G is the set of girders. p Production processes; p ∈ P, P is the set of production processes. l Logistics processes; l ∈ L, L is the set of logistics processes. i Installation processes; i ∈ I, I is the set of installation processes. e Equipment resources; e ∈ E, E is the set of equipment. c Human resources (i.e., crews); c ∈ C, C is the set of crews.

T gp
Duration spent by a girder g in a production process p. T gl Duration spent by a girder g in a logistics process l. T gi Duration spent by a girder g in an installation process i. W gp Time waited by a girder g for a production process p. W gl Time waited by a girder g for a logistics process l. W gi Time waited by a girder g for an installation process i. CIC Daily indirect or overhead costs.

MCE pe
Mobilization cost of an equipment e used in a production process p.

MCE le
Mobilization cost of an equipment e used in a logistics process l.

MCE ie
Mobilization cost of an equipment e used in an installation process i.

MCC pc
Mobilization cost of a crew c working on a production process p.

MCC lc
Mobilization cost of a crew c working on a logistics process l.

MCC ic
Mobilization cost of a crew c working on an installation process i.

FC pe
Fixed cost of an equipment e used in a production process p.

RE pe
Hourly cost of an equipment e used in a production process p.

RE le
Hourly cost of an equipment e used in a logistics process l.

RE ie
Hourly cost of an equipment e used in an installation process i.

TE pe
Duration spent by an equipment e working on production processes P.

TE le
Duration spent by an equipment e working on logistics processes L.

TE ie
Duration spent by an equipment e working on installation processes I.

RC pc
Hourly cost of crew c working on a production process p.

RC lc
Hourly cost of crew c working on a logistics process l.

RC ic
Hourly cost of crew c working on an installation process i.

TC pc
Duration spent by crew c working on production processes P.

TC lc
Duration spent by crew c working on logistics processes L.

TC ic
Duration spent by crew c working on installation processes I.
φ pc Overtime cost adjustment factor of crew c working on production processes P.
φ lc Overtime cost adjustment factor of crew c working on logistics processes L.
φ ic Overtime cost adjustment factor of crew c working on installation processes I. CS Hourly storage cost.
ne pe , ne le , and ne ie The maximum number of equipment e available for production p, logistics l, and installation i processes, respectively.
nc pc , nc lc , and nc ic The maximum number of crew c available for production p, logistics l, and installation i processes, respectively.

PYD min and PYD max
The minimum and maximum distances between the construction site and available locations for the production facility.

YS min and YS max
The minimum and maximum storage capacity of the production yard.

ST min and ST max
The minimum and maximum storage time of girders at the production yard.

DW H min and DW H max
The minimum and maximum number of daily working hours.

WWD min and WWD max
The minimum and maximum number of working days per week. Number of equipment resource e used in a production process p. NE le Number of equipment resource e used in a logistics process l. NE ie Number of equipment resource e used in an installation process i.

NC pc
Number of crew c working in a production process p.

NC lc
Number of crew c working in a logistics process l.

NC ic
Number of crew c working in an installation process i.

YS
Storage capacity of the production yard.

ST
Storage time of girders at the production yard.

PYD
The distance between the construction site and available locations to set up the production yard.
The decision variables of the optimization model cover a number of resource, logistics, and overtime planning decisions that control the project's duration and cost. They include: (1) the amount of equipment and human resources required for the production, logistics and installation processes; (2) the number of daily working hours and the number of working days per week; (3) the location of the production yard; and (4) the storage time and capacity. In total, the model includes 11 types of decision variables, as shown in Table 2. All of these variable types are integers, except one binary variable related to a piece of equipment used during the production stage. This variable represents the decision of the use of either a steaming machine to accelerate the concrete curing process or the traditional curing method.

Objective Functions and Constraints
In this study, the objective is to minimize both the project duration and the total costs. Equations (1) and (2) calculate the project duration in working (DWD) and calendar days (DCD), respectively. Equation (1) shows that the project duration includes the sum of the working duration of each production, logistics and installation process (T p , T l , and T i , respectively), plus the times waited by each girder due to the unavailability of resources in each production, logistics and installation process (W p , W l , and W i , respectively). The waiting times are obtained from the simulation model of the OSC project, as illustrated in the next section. On the other hand, the working duration of the processes is represented by probability distributions fitted from the historical data of previous projects, except for the working duration of logistics processes (T l ). T l depends on the distance between the production facility and the construction site (PYD), and how long the girders stay at the storage area (ST).
Regarding the project's cost (TPC), Equation (3) shows that it equals the sum of the direct (DC) and indirect costs (IC). Equation (4) shows that the indirect costs include the overheads (CIC) and the mobilization and demobilization costs (MC) of each piece of equipment and human resource used in the production, logistics and installation processes. The latter is calculated using Equation (5). The second cost component is the direct cost (DC) associated with the production, logistics and installation processes. It consists of three cost components, as shown in Equation (6). The first is the sum of the direct costs of the equipment resources assigned to accomplish the production, logistics and installation processes (CE), as illustrated in Equation (7). The second cost component is the sum of the direct costs of the human resources assigned to conduct the processes of the three supply chain stages (CC), as shown in Equation (8). The third cost component is the storage costs (SC), as shown in Equation (9). The previous cost estimation equations were adapted from Salimi et al. [10] and Marzouk et al. [35].
Because there is an inverse relationship between the two objectives (i.e., the total project's duration and costs), they are integrated into a non-dimensional fitness function using the function transformation method shown in Equation (10) [98].
where f t (x) is the fitness value of solution vector x (i.e., a configuration of the decision variables); DCD(x) and TPC(x)are the project duration and cost corresponding to solution x, respectively; DCD * and TPC * are the minimum values of the project duration and cost, respectively; and W t and W c are the relative weights of the project duration and cost, respectively (note that W t + W c = 1). The values of these weights (i.e., W t and W c ) are project dependent [35]. In some cases, delivering the project in a shorter duration is prioritized over the reduction of the total costs, such as the case of constructing field hospitals to face COVID-19 [99]. In other cases, such as the limited budget, the reduction of the total costs does matter more than shortening the project duration.
Indeed, it is commonly known in the literature that the domain of each decision variable is defined in a single equation as a constraint in the optimization model [100,101]. Following this concept results in the definition of the decision variables of the proposed optimization model in Equations (11)- (22). Usually, in OSC projects, the availability of equipment and human resources at each supply chain stage (i.e., production, logistics and installation) is limited. Besides this, increasing the number of resources at the production and construction stages impacts their productivity and safety. Therefore, the number of resources assigned to each production, logistics and installation process must be less than a maximum number defined in each supply chain stage, as shown in the constraints (11)(12)(13)(14)(15)(16).
Constraints (17) and (18) ensure that the location of the production facility (PYD) is within the allowable range and that its storage capacity (YS) does not exceed the specified limit. Constraint (19) limits the storage time of the girders inside the storage area. Constraints (20) and (21) ensure that the number of daily working hours (DW H) and the working days per week (WWD) are within the allowable ranges. Finally, constraint (22) indicates whether a steaming machine to accelerate the concrete curing is used or not during the production stage.

Solution Methods
This section elaborates on the methods adopted to achieve the research objectives. Firstly, the adopted SO approach, and its simulation and optimization components are explained in Section 4.1. Secondly, the methods used to reduce the computation time of the SO are discussed in Section 4.2.

The Adopted Simulation Optimization Approach
Simulation and optimization methods were integrated over decades using various approaches to solve complex and stochastic problems. Interested readers can refer to studies by Amaran et al. [102], Tekin and Sabuncuoglu [103], Wang and Shi [104], and Juan et al. [26] for more information on the classifications and terminologies established in this domain. This study adopts an SO approach in which the simulation model performs as an evaluation function (EF) in the optimization process [26]. Figure 4 shows the general framework of this SO approach. This approach consists of two parts: the optimization part and the simulation part. The optimization part starts by defining the decision variables of the optimization problem and the range of their allowable values. Then, an initial population of potential solutions of size P is generated randomly. Next, these generated solutions are identified to the simulation model in order to evaluate their fitness values (i.e., the simulation model acts as an EF). Note that each solution represents a vector of values of the input parameters associated with the simulation model. Each solution j is evaluated by a number of replications R to address the uncertainty of the activities' durations in the simulation model. The mean of these replications' results represents the fitness value of each solution. After the estimation of the fitness values of the generated solutions, the optimization part uses these values and the previously generated population to generate the second population of solutions based on the metaheuritic's logic. This cyclic process between the optimization and simulation parts is repeated until the stopping condition of the optimization part is met. Finally, the optimization part returns the best solution found so far and its fitness value.

Simulation Model Development
As explained in Section 3.1, the infrastructure OSC project under study represents a system that consists of a number of sequential processes wherein the precast girders represent passive entities serviced in these processes. Besides this, there is a complex interaction between these processes and different resources. The precast girders compete for these resources in order to move from one process to the next, forming queues inside the system. The representation of this type of complex and dynamic system characterized by queues using analytical models is cumbersome; DES is ideal for the modelling of such systems [23]. Figure 5 shows the DES model of a precast full-span construction project built in the MATLAB environment. The model was built using SimEvents, a generalpurpose piece of simulation software available in the MATLAB Simulink Library [105]. As shown in Figure 5, eleven blocks of SimEvents are used to model the system under study. These blocks are: Entity Generator, Entity Queue, Entity Server, Resource Pool, Resource Acquirer, Resource Releaser, Composite Entity Creator, Composite Entity Splitter, Entity Input Switch, Entity Terminator, and Stop. Initially, the Entity Generator block creates entities that represent the precast girders. Then, the precast girders undergo the different processes/activities modelled by the Entity Server blocks. Each Entity Server block delays the entity by the duration of its respective process. However, for the entities to undergo each process, the resources required to accomplish this process need to be assigned to these entities. The Resource Acquirer blocks are responsible for performing this function by seizing the required equipment and human resources stored in the Resource Pool blocks. Because any entity cannot depart from the Resource Acquirer blocks until it has seized all of the required resources, these blocks consider the waiting times due to the potential unavailability of the required resources. After finishing each process, its associated resources are released to be assigned to other entities. This function is conducted using the Resource Releaser blocks. Note that the system includes other types of entities to represent transport resources, such as trailers and the trolley. Once these resources are created using the Entity Generator blocks, they wait for the incoming precast girders at the Entity Queue blocks. Upon loading the incoming precast girder to the waiting trailer or trolley, both entities (i.e., the precast girder and the trailer or trolley) are combined into one entity using the Composite Entity Creator block. After the girder's unloading, this combined entity is split into its two original entities (i.e., the precast girder and the trailer or trolley) using the Composite Entity Splitter blocks. The unloaded trailer and trolley return to wait for another incoming girder using the Entity Input Switch blocks. After the girder is installed at its final destination, its representative entity is terminated using the Entity Terminator block. After installing all of the girders of the project, the simulation model is stopped using the Stop block.

Swarm Intelligence Metaheuristics
This section illustrates the five SI metaheuristics (FA, GWO, NBA, WOA and SSA) used in the optimization part of the adopted SO approach, as well as the controllable parameters of each metaheuristic. These SI metaheuristics were selected for the investigation of their performance because they are among the recent SI metaheuristics that gave satisfactory results, and they showed robustness in different optimization problems. These parameters balance exploration (globally exploring the search space through a randomized search) and exploitation (locally searching in the promising areas obtained from the exploration process). At the end of this section, the methods used to tune the controllable parameters of the SI metaheuristics are discussed. The five SI metaheuristics were coded in the MATLAB environment.

Firefly Algorithm (FA)
Inspired by the way in which fireflies communicate with each other for mating, Yang [106] proposed a stochastic population-based metaheuristic. The pseudo-code of the FA is shown in Figure 6. The first step is to generate randomly the positions of all of the fireflies in the population. Secondly, we evaluate the brightness (fitness value) of each firefly. Thirdly, we move the fireflies with high fitness values to those with lower fitness values in the case of minimization problems by using Equation (23). It is worthwhile to note that the attractiveness β depends on the distance r between the fireflies, as indicated in Equation (25). After that, the brightness of the new fireflies is evaluated. Finally, these steps are repeated until we meet the defined terminating conditions.
where x i = the position of a firefly attracted to a brighter firefly x j ; r ij = the distance between two fireflies i and j calculated by Equation (24); β 0 = the attractiveness at r = 0, usually equalling 1. γ = the light absorption coefficient; α = a constant number between 0 and 1; and u = a random number between 0 and 1. The controlling parameters of the FA are parameters γ and α. The convergence speed of the FA is greatly affected by the values of γ. However, α controls the exploration of the FA.

Grey Wolf Optimization (GWO) Algorithm
The GWO algorithm, inspired by both the leadership hierarchy and the hunting behaviour of grey wolves, was developed by Mirjalili et al. [33]. The pseudo-code of the GWO algorithm is shown in Figure 7. Firstly, all of the positions of the grey wolf population (i.e., solutions) are randomly generated. Secondly, the simulation model is called upon to evaluate the fitness of each solution. Thirdly, the best three solutions ( where → X α = the position of the first best solution in iteration i;  (29) and (30).
where → a = a parameter decreased linearly from 2 to 0 over the iterations according to Equation (31)  Meng et al. [107] improved the performance of the original Bat Algorithm (BA) by developing the NBA. The pseudo-code of the NBA is presented in Figure 8. Firstly, the positions and velocities of all of the bats in the population are initialized randomly. Secondly, the simulation model is used to evaluate the fitness of each bat. Thirdly, the position of each bat is updated based on its habitat selection. The bats choose randomly between the quantum behaviour (Equation (32)), to forage in different habitats, and the mechanical behaviour (Equations (34)-(37)), to forage in limited habitats. Fourthly, a local search is implemented by generating a solution around the best solution using Equations (39) and (40). Fifthly, the fitness value of the new solutions is calculated, and the parameters of the NBA are updated according to Equations (41)- (44). Finally, the previous steps are iterated until the satisfaction of the stopping criteria.
where x i,j = the position of bat i in a D-dimensional search space and j ∈ {1, 2, , D}; t = the current iteration; g t = the global best position; θ = the contraction-expansion coefficient calculated at each iteration k over the number of iterations K; Maxθ and Minθ = the maximum and minimum values of θ; mean t = the average of the bats' positions; and u and rand = random numbers between 0 and 1.
where f i,j = the frequency of bat i at dimension j; f max and f min = the maximum and minimum values of the frequency; C i = the compensation rate for the Doppler effect in echoes randomly selected between 0 and 1 for each bat i; c = the speed in the air (340 m/s); ε = the smallest constant in the computer; v i,j = the velocity of bat i at dimension j; v g,j = the velocity of the best solution at dimension j; w = the inertia weight given a value between Maxw and Minw; and rand = a random number between 0 and 1.
x t+1 i,j = g t j * 1 + rand n 0, σ 2 (39) where rand n 0, σ 2 = a Gaussian distribution with mean 0 and standard deviation σ 2 ; A t i = the loudness of the bat i; and A t mean = the average loudness of all of the bats.
where f (x i ) = the fitness value of bat i; f (x) = the fitness value of the best solution found so far; A t i = the loudness of bat i at iteration t; r t i = the pulse rate of bat i at iteration t; and α and γ = constants, which are usually equal 0.9. The parameters that affect the performance of the NBA are parameters α, γ, P, G, C, w and θ. α and γ affect the convergence speed of the NBA. P affects the bats' choice to either adopt quantum behaviour or mechanical behaviour. The tuning of G is required to ensure exploration and exploitation. Based on some experiments, Meng et al. [107] suggested that C can be in the interval (0.1,1). w has the same effect of its counterpart in the PSO algorithm, and it is recommended that w is decreased linearly from 0.9 to 0.4. Similarly, θ is recommended to be decreased linearly from 1 to 0.5.

Whale Optimization Algorithm (WOA)
Inspired by the hunting behaviour of humpback whales, Mirjalili and Lewis [27] provided the WOA to solve single-objective optimization problems. The pseudo-code of the WOA is depicted in Figure 9. First of all, the positions of all of the whales in the population are randomly generated. Secondly, the fitness of each whale is calculated by using the simulation model. Thirdly, the position of the fittest whale is identified. Fourthly, each whale chooses to update its position by using either a spiral movement (Equation (46)) or a circular movement (Equations (48) or (53)). This choice is based on a probability of 50%. If a whale chooses the circular movement, the position of this whale is updated using the position of the fittest whale identified so far (Equation (48)), if the absolute value of A is smaller than one; A is a parameter in the WOA updated by Equation (49). However, if the absolute value of A is larger than or equal to one, the whale's position is updated by using the position of a randomly chosen whale (Equation (53)). Finally, the previous steps are repeated until the satisfaction of the predefined terminating conditions.
where → X(i) = the position of the fittest whale found so far in iteration i; → X(i) = the position of a solution in iteration i; → D = the distance between a whale and the food (the fittest whale found so far); b = a positive number used to define the shape of the spiral; l = a random number in the interval of (−1,1); and | f |= the absolute value of f .
where → a = a decreased variable from 2 to 0 over the iterations, updated by Equation (51), as proposed by Zhong and Long [108]; → r = a uniformly distributed random number between 0 and 1.
where µ = the adjustable factor; i = the current iteration; I = the maximum number of iterations.
where → X rand = the position of a randomly chosen whale from its population. The SSA, as proposed by Mirjalili et al. [110], is a nature-inspired metaheuristic optimization algorithm inspired by salps' intelligent foraging behaviour. Salps are aquatic creatures similar to jellyfish. The pseudo-code of the SSA is shown in Figure 10. Firstly, the position of each salp in the population is randomly initialized. Secondly, the simulation model is used to evaluate the fitness of each individual salp. Thirdly, the algorithm identifies the salp position with the best fitness value and assigns it to variable F, which represents the food source that will be pursued by the followers. Fourthly, the controlling parameter C 1 is updated using Equation (54). Finally, the position of the leader is updated by using Equation (55), and the positions of the followers are updated using Equation (56). These steps are repeated until the specified stopping criteria are met.
where C 1 is a controlling parameter, L is the maximum number of iterations, and l is the current iteration.
where x j = the position of the leader in the jth dimension; F j = the position of the food source in the jth dimension; ub j = the upper bound of jth dimension; lb j = the lower bound of jth dimension; and c 2 and c 3 = random numbers generated from a uniform distribution in the interval of (0,1).
where x i j = the position of ith follower salp in the jth dimension. The SSA is controlled by only one parameter-C 1 -besides the population size and the number of iterations. C 1 is scheduled throughout the iterations to ensure both diversification and intensification.

Parameter Selection of the Metaheuristics
Unquestionably, any metaheuristic performance is susceptible to the values of its parameters. In this study, the parameter settings recommended by the developers of the SI metaheuristics are used initially. In order to improve the solution quality of each metaheuristic, their parameters are tuned by conducting many experiments. The final parameter settings of each SI metaheuristic are listed in Table 3. In order to ensure a fair comparison between the SI metaheuristics, each metaheuristic has to use the same number of objective function evaluations [111]. Consequently, the population size of all of the metaheuristics is 50 individuals, and the maximum number of iterations is 50. As shown in the third column of Table 3, offline tuning, the pre-scheduled variation of parameter settings, and the self-adaptation of the parameter settings are three approaches for parameter tuning. The offline tuning is conducted by trial-and-error, where the parameters' values are kept constant throughout the iterations. However, it is time-consuming, and depends on user experience. Furthermore, it ignores the interaction effect among the parameters. Therefore, online tuning methods are developed for the automatic tuning of parameters through iterations with minimal human interference [112]. The online tuning methods used in this study are the pre-scheduled variation and self-adaptation methods. In the pre-scheduled variation method, the parameters' values change through the iterations according to a predefined function. In the self-adaptation tuning method, the parameters are tuned by the metaheuristic itself, similarly to other decision variables.

Reduction of the Computation Time
Two approaches were adopted to provide the project manager with accurate planning solutions in a reasonable time. These approaches are CRN and parallel computing. An explanation of each approach and its implementation is given in the following subsections.

Common Random Numbers (CRN)
CRNs are one of the most popular variance reduction techniques (VRTs) used to reduce the variance of the random output from the simulation model without increasing the number of simulation replications [113]. The CRN is suitable for SO, especially for the comparison of two or more solutions. Its main idea is that the solutions generated at each iteration should be evaluated in the simulation model using the same generated random variates. The CRN is successfully implemented by synchronizing random numbers without affecting the independence between the simulation replications and the stochastic processes. However, there is no guarantee that the CRN could reduce the solutions' variability; it may even backfire. Therefore, Law et al. [113] recommended the implementation of a preliminary analysis before applying CRN. This preliminary analysis tests its efficacy through four performance measures represented by Equations (57)-(60). The first performance measure (Equation (57)) states that the variance of the difference between the fitness values of two candidate solutions (S 2 Z (N)) should be reduced after applying CRN. The second one (Equation (58)) says that this value (S 2 Z (N)) after using CRN should be less than the sum of the variance of the two candidate solutions. The third measure explains that the confidence interval's half-width should be reduced after using CRN (Equation (59)). This reduction in the half-width can reduce the number of replications [114]. The fourth (Equation (60)) states that the fitness values of the two candidate solutions (X 1 , X 2 ) should be positively correlated which results into a positive covariance (Cov(X 1 , X 2 ) > 0); thus, the variance of Z(N) could be reduced according to Equation (61).
Corr(X 1n , X 2n ) CRN > 0 (60) where S 2 Z (N) CRN and S 2 Z (N) no CRN = the sample variance of Z n (Z n is the difference between the fitness values of two candidate solutions (X 1n , X 2n )) with and without using the CRN for N replications, respectively; S 2 1 (N) and S 2 2 (N) = the variance of the fitness values of two candidate solutions (X 1n , X 2n ), respectively, for N replications; h CRN Z (N) and h no_CRN Z (N) = the half-width of the 95% confidence interval of Z n with and without using the CRN, respectively, for N replications; Corr(X 1n , X 2n ) CRN = the correlation between X 1n and X 2n after using CRN; Var(y) = the variance of y; Cov(X 1n , X 2n ) = the covariance of X 1n and X 2n .
Two solutions-called solutions A and B, as shown in Table 4-were arbitrarily chosen to study CRN's feasibility. Table 5 shows the statistical results obtained before and after applying CRN for solutions A and B for 100 replications. Table 5 indicates that the variance of the difference between the fitness value, the project's duration and the project's cost of solutions A and B is reduced by about 48%, 41% and 65%, respectively, after applying CRN. Table 5 shows that the variance of Z n (S 2 Z (N) = 0.0028) after using CRN is less than the sum of the variance of the fitness value of solutions A and B (0.0033 + 0.0031 = 0.0064). Furthermore, using CRN leads to a reduction in the half-width of the fitness value, the project's duration and the project's cost by about 28%, 23% and 41%, respectively. This reduction in the half-width could be used to reduce the required number of replications. For example, only 100 replications are required to reach a half-width of the fitness value equal to 0.0103 when using CRN. However, more than 195 replications are required to reach the same half-width if we do not use CRN. Hence, the computation time can be reduced significantly because one replication takes about 10 s on a workstation (Intel(R) Xeon(R) CPU @ 2.27 GHz, 40.0GB Random Access Memory (RAM)). Finally, using the CRN results in a positive correlation between solutions A and B. Based on these pilot study results, CRN can indeed reduce the variance between the generated solutions.
After the demonstration of the efficacy of CRN, the confidence interval method with a 95% confidence level was used to determine the minimum number of simulation replications required to obtain results with no more than 5% error in their means (i.e., H X < 0.05, where X = the sample mean and H is the half − width). By conducting multiple simulation experiments on different solutions, four replications using CRN were determined to achieve this level of reliability.

Parallel Computing
The basic idea of parallel computing is to divide the total computational load of a given problem into independent parts and allocate them to available processors to be solved simultaneously in order to reduce the computation time [115]. Parallel computing can be implemented on a single computer with multiple processors, a network of connected computers, or a combination of both [116]. In this study, a single computer with multiple processors was selected to implement parallel computing because the other options are rarely available to construction planners.   Four strategies are available to implement parallel computing in metaheuristics. These strategies are low-level, search space decomposition, independent multi-search, and cooperative search strategies [117]. The low-level strategy is usually implemented using the master-slave parallel programming paradigm, as shown in Figure 11. In this paradigm, the master core generates solutions for the initial population and distributes these solutions among the slave cores. Then, the slave cores evaluate the dispatched solutions and send their fitness values to the master core. After that, the master core uses these solutions' fitness values to generate the second generation of solutions according to the metaheuristic logic. This process is repeated until the satisfaction of the terminating conditions of the used metaheuristic. Therefore, there is no communication among the slave cores. For example, suppose the metaheuristic population P consists of 50 solutions, and five slave cores (i.e., N = 5) are available besides the master core. In this case, the process starts at the master core, which randomly generates the initial 50 solutions. The master core then distributes the 50 solutions equally among the slave cores. After that, each slave core evaluates the fitness value of the ten solutions assigned to it using the simulation model. Each solution must be evaluated by a number of replications R to account for the model's uncertainty. Then, the mean fitness value of these replications for each solution is calculated. The five slave cores can evaluate the assigned solutions in parallel. After evaluating the mean fitness values of the ten assigned solutions, the slave cores send the 50 solutions and their mean fitness values to the master core. Next, according to the applied metaheuristic logic, the master core uses the received information to generate 50 new solutions for the new population. This process is iterated until it meets the stopping conditions. As shown in Figure 11, the low-level strategy does not require any modifications to the metaheuristic logic, or the search space, unlike the other parallelization methods. Consequently, the low-level strategy is selected to implement parallel metaheuristics. Figure 11. The master-slave paradigm for the parallelization of SO.

Model Implementation
The case study of the precast full-span bridge deck construction using a launching gantry provided by Mawlana [118] was used in this study. This case study is a bridge which consists of 35 equal spans, with a total length of 875 m. Table 6 shows the probability distributions of each activity duration. As shown in Table 6, the activity durations are represented by triangular distribution. Triangular distribution has been commonly used to represent construction activity durations because it requires only three parameters (i.e., low, mode, and high) that experts can easily estimate if there is no historical data [119]. Furthermore, McCabe [120] indicated that triangular distribution provides a good approximation of the beta distribution, which was demonstrated by AbouRizk and Halpin [119] as an adequate distribution to fit the historical data of activity durations. In this case study, planning for 14 resources, overtime planning, and logistics decision variables considerably affect the project duration and cost. The ranges of the 14 decision variables are listed in Table 7. The chosen overtime policy represented by the number of daily working hours and the number of working days per week affects the workers' productivity and the activities' durations. The detailed information related to each overtime policy, and the hourly cost of the different resources and other cost parameters are available in the Supplementary Materials. In this study, it is assumed that both the project duration and cost have equal weights (i.e., W t = W c = 0.5), meaning that both objectives (i.e., shortening the project duration and saving the total costs) have the same priority [121][122][123]. This information was identified for the developed DES model built using SimEvents (Version 9.6.0 (R2019a)). In this study, SimEvents was selected to model the system under study because of its easier compatibility with parallel computing paradigms than special-purpose construction simulation software [10]. Before the integration of the developed simulation model with the SI metaheuristics, it was validated based on the results obtained by Mawlana and Hammad [32]. They simulated the same project using STROBOSCOPE simulation software. Table 8 compares the results of ten different solutions evaluated by the two simulation models (STROBOSCOPE and SimEvents). An ANOVA F-test was conducted to determine if there is a significant difference between the results of the two simulation models. The obtained p-values for the time and cost are 0.327 and 0.078, respectively. At a significance level of α = 0.05, there is no significant difference between the results of the two simulation models. Table 6. Activities durations of the precast full-span bridge deck construction using a launching gantry [32].    2  76  266  77  271  3  77  247  77  244  4  80  214  80  217  5  81  207  81  208  6  83  205  83  207  7  84  197  84  199  8  89  195  89  202  9  92  193  92  192  10  98  167  97  167 After selecting the low-level strategy, a pilot study is required to determine the optimum number of cores to be used for parallel computing. The number of cores is increased from one to twelve, as shown in Figure 12. Twelve is the maximum number of cores available in the used device. The corresponding time to evaluate the fitness values of a population of 50 individuals is recorded. The results show that the minimum computation time is achieved using ten cores. Using this number of cores reduces the computation time by 87.0% compared with using only one core, which represents the case of not applying the parallel computing. This significant reduction in the computation time matches the results of studies by Mawlana and Hammad [32] and Salimi et al. [10], who indicated a reduction in the computation time by 90.5% and 95.1%, respectively, after using parallel computing. Figure 12 shows that increasing the number of cores does not always reduce the computation time. For example, increasing the number of cores from 11 to 12 increases the computation time. This result is because increasing the number of cores reduces the computation load assigned to each core, whereas the communication overhead increases. Moreover, the even distribution of the computation load to the available cores is an important factor in the reduction of the computation time. For example, the distribution of 50 solutions to ten cores, as shown in Figure 12, lead to the shortest computation time.

Results and Discussion
This section presents a qualitative and quantitative comparison between the five SI metaheuristics after using CRN and parallel computing to reduce the computation time. The stopping condition adopted in this study is reaching the predefined maximum number of iterations (i.e., 50 iterations). Using this stopping condition, each SO run of each SI metaheuristic took on average 11 h on a workstation (Intel(R) Xeon(R) CPU @ 2.27 GHz and 40.0 GB Random Access Memory (RAM)). Note that solving stochastic SO takes a longer computation time than usual because each solution in each population needs to be evaluated by the number of simulation replications to account for the model's uncertainty [10]. In order to address this issue, commercial SO software such as OptQuest combines neural networks with metaheuristics to speed up the SO process [26,124]. The comparisons between the SI metaheuristics are based on five independent runs with different random variates of each SI metaheuristic. Firstly, the analysis of each SI metaheuristic convergence behaviour is discussed in Section 6.1 to provide a qualitative comparison between the SI metaheuristics. In order to measure the competitiveness quantitatively among them, a statistical analysis of each SI metaheuristic is provided in Section 6.2.

Analysis of the Convergence Behaviour
In the literature, the convergence curve is the most common way to assess the performance of an optimization algorithm qualitatively in terms of its ability to achieve exploration and exploitation. Each line graph in Figure 13 shows the convergence curve of each SI metaheuristic.
Generally, Figure 13 shows that the population fitness of the five SI metaheuristics is improved throughout the iterations. As shown in Figure 13, the NBA exhibits steep convergence during the first iterations. However, the search process is stuck in a local optimum point which is better than the solutions obtained by the WOA, FA and SSA. The performance of the NBA demonstrates its ability to detect the promising regions in the search space early and exploit them. On the one hand, this performance might be beneficial because good-quality solutions can be obtained in a reasonable time; however, being stuck in local optima prevents the optimizer from exploring other promising search regions. Despite the fact that GWO shows no improvements for some iterations, it achieves a better balance between exploration and exploitation compared to NBA. This balance between exploration and exploitation prevents GWO from trapping into a local optimum early. Furthermore, GWO shows a relatively rapid convergence in the beginning and reaches a better solution than NBA at the end of iterations. As for FA, its performance is characterized by smooth and steady convergence. Despite its ability to avoid trapping in local optima, it could not obtain better results at the end of its iterations compared with the NBA and GWO. These results indicate that FA might need more iterations to reach more optimum search regions and overcome the NBA and GWO. Furthermore, it could not give good results in the beginning due to its slow convergence. Unlike FA, WOA shows rapid convergence in the beginning, but stagnation in local optima for a large number of iterations made the performance of WOA unappealing. This means that the exploration process in WOA should be improved to avoid such a performance. Finally, SSA neither shows rapid convergence at the beginning of the search process nor reaches better results at the end of the search process.

Statistical Results of the SI Metaheuristics
Although the convergence analysis provides a qualitative comparison among the SI metaheuristics, it cannot assess the extent to which each metaheuristic outperforms its counterparts. For this purpose, some statistical measures were used to measure the central tendency and the variability of the obtained results. Furthermore, the half-width of the 95% confidence interval of the mean was calculated. The Tukey method, a multiple-comparison method in ANOVA, was used to detect which metaheuristic provides results' means which are significantly different from the means of other metaheuristics, and hence to rank the SI metaheuristics with an overall confidence level of 95%. Table 9 summarizes these statistical measures of the project duration obtained by the SI metaheuristics. It can be noticed that the minimum project duration resulting from the optimization algorithms is 79 days. The NBA and WOA find this minimum value. On the other hand, 85 days is the most prolonged project duration recorded from the numerical experiments, and it is obtained by SSA. Figure 14 shows the box and whisker plot of the project durations obtained by each metaheuristic. WOA and SSA provide the lowest mean of the recorded project durations, whereas FA gives the highest mean.
The means of the project durations resulting from GWO and NBA are 81.8 and 82 days, respectively. The values of the percentage error, which is the subdivision of the half-width of the 95% confidence interval of the project duration and the project duration mean (i.e., H X < 0.05, where X = the sample mean and H is the half − width) show that the number of optimization runs results in quite reliable outcomes (below 5%), as shown in the last row of Table 9.   Table 10 lists the statistical results of the project cost for each SI metaheuristic. NBA gives the minimum recorded project cost (i.e., 1,578,000 $), whereas SSA finds the highest project cost. Figure 15 shows that the solutions obtained by GWO have the lowest mean project cost. However, SSA produces solutions with the highest mean project cost. The average costs obtained by NBA and WOA are 1,638,900 $ and 1,650,100 $, respectively. However, higher cost means are reported from FA. The percentage error values in Table 10 show again that the number of optimization runs is sufficient.  For instance, the SI metaheuristics which obtain solutions with a low mean of the project duration have a high mean of the project cost, and vice versa. Therefore, the fitness values of the obtained solutions are used to rank the SI metaheuristics. Figure 16 shows the box and whisker plot of the fitness values of the solutions obtained by the SI metaheuristics.
The results indicate that GWO and NBA produce the best results. For WOA and FA, their generated solutions have higher fitness values. On the other hand, the solutions obtained by SSA have the highest fitness values (the worst results). Besides this, the solutions obtained by GWO have the least variability, followed by the NBA. However, the SSA produces the most scattered solutions. The variability among the other SI metaheuristics solutions is quite similar. In order to make a decisive comparison among the means of the optimization algorithms' fitness values, Tukey's method was selected to conduct the pairwise comparisons [125]; its results are shown in Figure 17. The results show that, with a 95% confidence level, there is no significant difference between the fitness of the solutions generated by GWO, NBA, and WOA. The FA solutions are worse than those obtained by GWO and NBA, while SSA provides the worst solutions among the SI metaheuristics.  Given the above analysis, it can be concluded that GWO, the NBA, and the WOA can find the best solutions, and there is no statistically significant difference between the quality of their solutions. Given the relatively new advent of these three metaheuristics (i.e., GWO, NBA, and WOA), few studies have reported their efficiency in the solution of different construction planning problems such as the optimization of construction duration and schedule robustness using GWO [126], and the construction stage and zone optimization of Rockfill dams using the WOA [127]. On the other hand, the FA generates less optimal solutions than those produced by the others, while the SSA provides the worst results. Regarding the convergence behaviour, the NBA shows rabid convergence behaviour compared with the other optimization algorithms. However, GWO shows the best balance between exploration and exploitation to avoid stagnation into a local optimum, and the FA exhibits smooth and steady convergence behaviour. However, the SSA suffers from becoming stuck in local points during the search process. The performance of the SSA can indeed be enhanced using other equations to control its only control parameter (i.e., C 1 ).

Conclusions
As a new construction method, OSC has the potential to improve the sustainability of the built environment. However, its wider diffusion is hindered by multiple barriers; among them is the need to make extensive planning decisions. This barrier calls for computational methods that help project managers to make the optimum planning decisions. Simulation and optimization are among the key computational methods that can capture the complexity and uncertainty of OSC projects. This study investigated, for the first time, recent SI metaheuristics, namely: the FA, GWO, the NBA, the WOA and the SSA, for the integration of simulation and optimization, and applied them to infrastructure OSC projects to simultaneously minimize projects' duration and cost.
These SI metaheuristics were applied to a case study of a bridge deck construction project using the OSC method. The construction operations were simulated via a DES model. In order to reduce the computation time of the SO process, CRN and parallel computing were integrated into the SO models, reducing the average computation time by 87.0%. Then, a comparative study based on the convergence behaviour and the statistical analysis of the five SI metaheuristics was conducted. The convergence behaviour analysis indicated that the NBA shows rabid convergence behaviour compared with the other SI metaheuristics. However, GWO shows the best balance between exploration and exploitation to avoid stagnation into local optima. The FA exhibits smooth and steady convergence behaviour, whereas WOA and the SSA suffer from becoming stuck in local points. Based on the 95% confidence level, the statistical analysis proved that there is no statistically significant difference between the solution qualities obtained by GWO, the NBA and the WOA. However, the FA and SSA provide less optimal solutions. This comparative analysis proves that the NBA and GWO are very competitive for the SO of OSC projects. Considering the time and effort needed to tune the controllable parameters of the SI metaheuristics, GWO is preferable to the NBA, because GWO has a smaller number of tuning parameters. Besides this, its simplicity makes it more accessible for researchers to modify its logic or hybridize it with other local search methods to further improve its performance. To summarize, this study contributes to the body of knowledge by comparing multiple and recent SI metaheuristics for the stochastic SO of TCTP in infrastructure OSC projects. Furthermore, the study incorporates CRN and parallel computing with the SO models to reduce the computation time. To the best of the authors' knowledge, this study is one of the first to integrate recent SI metaheuristics such as the FA, GWO, NBA, WOA and SSA with DES for the SO of OSC projects.
This study offers a valuable reference to professionals and academics who are interested in the optimization of the planning of OSC projects. The integration of CRN and parallel computing with the most recent and powerful SI metaheuristics such as GWO, the NBA and the WOA can provide project managers with near-optimum resource planning and logistics decisions to reap the full sustainability merits of OSC. This study could be seen as one of the early attempts to evaluate the performance of multiple SI metaheuristics in the solution of a set of optimization problems characterized by uncertainty in the construction industry [26]. Hence, it can guide researchers to test the performance of other metaheuristics for the SO of different construction projects.
Nonetheless, this study's conclusions should be interpreted in the light of some limitations. According to the no-free-lunch theorem, the reported results herein cannot be generalized to other types of construction projects. Furthermore, the developed simulation model ignores some non-physical aspects of OSC projects, such as organizational policies, rework, and the laborers' skill level. These aspects could be considered using SD, forming a hybrid DES-SD model. Such hybrid models could provide another challenging testbed for SI metaheuristics to optimize real-world systems characterized by complexity and uncertainty. Given the potential and simplicity of GWO in SO for OSC, its performance could be further improved by hybridizing it with local search methods such as the Hooke-Jeeves method and the Nelder-Mead simplex method. Furthermore, different kinds of OSC projects could be used to test SI metaheuristics.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/su132413551/s1, Table S1: The cost input parameters of the proposed model. Table S2