Joint Scheduling of Yard Crane, Yard Truck, and Quay Crane for Container Terminal Considering Vessel Stowage Plan: An Integrated Simulation-Based Optimization Approach

: The joint scheduling of quay cranes (QCs), yard cranes (YCs), and yard trucks (YTs) is critical to achieving good overall performance for a container terminal. However, there are only a few such integrated studies. Especially, those who have taken the vessel stowage plan (VSP) into consideration are very rare. The VSP is a plan assigning each container a stowage position in a vessel. It affects the QC operations directly and considerably. Neglecting this plan will cause problems when loading/unloading containers into/from a ship or even congest the YT and YC operations in the upstream. In this research, a framework of simulation-based optimization methods have been proposed ﬁrstly. Then, four kinds of heuristics/metaheuristics has been employed in this framework, such as sort-by-bay (SBB), genetic algorithm (GA), particle swarm optimization (PSO), and multiple groups particle swarm optimization (MGPSO), to deal with the yard crane scheduling problem (YCSP), yard truck scheduling problem (YTSP), and quay crane scheduling problem (QCSP) simultaneously for export containers, taking operational constraints into consideration. The objective aims to minimize makespan. Each of the simulation-based optimization methods includes three components, load-balancing heuristic, sequencing method, and simulation model. Experiments have been conducted to investigate the effectiveness of different simulation-based optimization methods. The results show that the MGPSO outperforms the others.


Introduction
Maritime transport accounts for 80% of global trade [1], which indicates the importance of this kind of transport. One advantage of maritime transport is lower unit transportation cost, which drives the sustainability of maritime transport. Specifically, maritime transport consists of various kinds of transports, such as container transport, bulk transport, and tanker transport, among which container transport is the most popular one due to the generality and operational efficiency. Further, it is noted that about 60% of the world's deep-sea general cargos are transported through container transport, and the routes between some countries are even containerized up to 90% [2]. Thus, maritime transport has attracted our attention.
Maritime transport helps the development of international trading. This transport is especially important for ocean countries. In the world maritime transport system, port used, though they are considered suitable for solving complicated system problems [9]. Thus, the simulation-based optimization methods have attracted our attention and are to be adopted in this research.
This research focuses on solving the simultaneous YCSP, YTSP, and QCSP by using simulation-based optimization approaches. A joint schedule for moving export containers from their current YSP positions to their VSP positions is developed with the YSP and VSP constraints taken into consideration. The objective is to minimize makespan. In the simulation-based approaches, heuristics/metaheuristics, including Method1 (sort-by-bay (SBB)), Method2 (genetic algorithm (GA)), Method3 (particle swarm optimization (PSO)), and Method4 (multiple groups particle swarm optimization (MGPSO)) have been used as sequencing methods. A small-sized experiment is used to demonstrate the solution found by Method4 (MGPSO). In addition, extensive experiments have been conducted to investigate the effectiveness of these developed approaches. Finally, a statistical t-test has been conducted to validate the robustness of the experimental results. The results showed that Method4 (MGPSO) outperforms the others.
The rest of this paper is organized as follows. Section 2 introduces background knowledge and relevant studies. Section 3 defines the problems that are formulated into a MILP. Section 4 proposes a framework for developing simulation-based optimization methods. Section 5 includes some experiments and analyses of the experimental results. Finally, Section 6 includes a conclusion and suggestions for future research directions. Figure 1 shows operations in the three main areas of a container terminal. The seaside area includes facilities such as berths for accommodating ships and QCs for loading and unloading containers. The yard side area stores containers before their transportation to the seaside or yard side areas for further processing. It is a buffer between the two areas. In this area, YCs are the main MHE for storing/retrieving containers into/from a block. YTs are the main transportation means used between the seaside and yard side areas. The landside area connects to the inland and controls the transportation of containers into and out of the container terminal. The main MHE used in the landside area includes inland trucks (ITs) and trains.

Overall View of Operations in a Container Terminal
Mathematics 2021, 9, x FOR PEER REVIEW 3 of 28 swarm optimization (PSO), have been widely used. According to [9], the simple heuristics (RULEs) and GAs top the list. However, it is found that the simulation methods have been seldom used, though they are considered suitable for solving complicated system problems [9]. Thus, the simulation-based optimization methods have attracted our attention and are to be adopted in this research. This research focuses on solving the simultaneous YCSP, YTSP, and QCSP by using simulation-based optimization approaches. A joint schedule for moving export containers from their current YSP positions to their VSP positions is developed with the YSP and VSP constraints taken into consideration. The objective is to minimize makespan. In the simulation-based approaches, heuristics/metaheuristics, including Method1 (sort-by-bay (SBB)), Method2 (genetic algorithm (GA)), Method3 (particle swarm optimization (PSO)), and Method4 (multiple groups particle swarm optimization (MGPSO)) have been used as sequencing methods. A small-sized experiment is used to demonstrate the solution found by Method4 (MGPSO). In addition, extensive experiments have been conducted to investigate the effectiveness of these developed approaches. Finally, a statistical t-test has been conducted to validate the robustness of the experimental results. The results showed that Method4 (MGPSO) outperforms the others.
The rest of this paper is organized as follows. Section 2 introduces background knowledge and relevant studies. Section 3 defines the problems that are formulated into a MILP. Section 4 proposes a framework for developing simulation-based optimization methods. Section 5 includes some experiments and analyses of the experimental results. Finally, Section 6 includes a conclusion and suggestions for future research directions. Figure 1 shows operations in the three main areas of a container terminal. The seaside area includes facilities such as berths for accommodating ships and QCs for loading and unloading containers. The yard side area stores containers before their transportation to the seaside or yard side areas for further processing. It is a buffer between the two areas. In this area, YCs are the main MHE for storing/retrieving containers into/from a block. YTs are the main transportation means used between the seaside and yard side areas. The landside area connects to the inland and controls the transportation of containers into and out of the container terminal. The main MHE used in the landside area includes inland trucks (ITs) and trains. A container terminal usually has to deal with three kinds of containers: import, export, and transshipment. These containers are moved among the seaside, yard, and landside areas. Planning storage positions for containers in different areas is an essential task for container terminal planners, and the changing of the positions of containers requires collaboration and cooperation among QCs, YTs, and YCs. A different operational procedure is used for different kinds of containers. For import containers, they are unloaded by QCs from ship firstly, transported by YTs subsequently, and stored in blocks by YCs. The A container terminal usually has to deal with three kinds of containers: import, export, and transshipment. These containers are moved among the seaside, yard, and landside areas. Planning storage positions for containers in different areas is an essential task for container terminal planners, and the changing of the positions of containers requires collaboration and cooperation among QCs, YTs, and YCs. A different operational procedure is used for different kinds of containers. For import containers, they are unloaded by QCs from ship firstly, transported by YTs subsequently, and stored in blocks by YCs. The export containers use the reverse procedure. Finally, these import containers are transported by ITs or trains in the landside area to their consignees. The transshipment containers use half of the import procedure and half of the export procedure. The YCs, YTs, and QCs are essential MHE in a container terminal, and they can affect the overall performance of the container terminal considerably. To best utilize these MHEs, a joint schedule for these MHEs is necessary.

Yard Storage Plan
A yard storage plan (YSP) is a plan arranging the storage positions for containers in a block. A container yard usually consists of multiple blocks; each being configured in the Asia or Europe type [8]. Figure 2 shows a storage block of Asia type with one YC equipped. Containers are stacked along the X, Y, and Z dimensions in the block, and each container's storage position is denoted as (x, y, z), where the x, y, and z are bay, row, and tier numbers, respectively. The x, y, and z are subject to the constraint Equation (1).
(1) export containers use the reverse procedure. Finally, these import containers are transported by ITs or trains in the landside area to their consignees. The transshipment containers use half of the import procedure and half of the export procedure. The YCs, YTs, and QCs are essential MHE in a container terminal, and they can affect the overall performance of the container terminal considerably. To best utilize these MHEs, a joint schedule for these MHEs is necessary.

Yard Storage Plan
A yard storage plan (YSP) is a plan arranging the storage positions for containers in a block. A container yard usually consists of multiple blocks; each being configured in the Asia or Europe type [8]. Figure 2 shows a storage block of Asia type with one YC equipped. Containers are stacked along the X, Y, and Z dimensions in the block, and each container's storage position is denoted as (x, y, z), where the x, y, and z are bay, row, and tier numbers, respectively. The x, y, and z are subject to the constraint Equation (1).
To access an export container, the YC and its hoist are used. The YC moves along the ±X direction to reach the bay position, and then the hoist moves along the ±Y direction to reach the row position of the target container. Both the YC and its hoist can move simultaneously. After arriving above the position of the target container, the hoist then lowers in the -Z direction to pick up the container along the +Z direction. Finally, the hoist moves the container to the YC waiting in the truck lane.

Vessel Stowage Plan
A vessel stowage plan (VSP) is a plan arranging the vessel stowage positions for export containers in a vessel. Figure 3 shows the storage positions of some bays in the vessel. The assignment of stowage positions for a container is based on the container's attributes, such as size, weight, and destination. In the vessel, each bay for storing 20-foot containers is assigned with an odd bay number, with the number increasing from the bow to the stern. In addition, two 20-foot bays can form one 40-foot bay assigned with an even bay number. For instance, together Bays 09 and 10 form Bay 11. In the VSP, each container i has a position denoted as ( , , ) [10], where the indicates a bay number; the indicates a row number; and the indicates a tier number. The bay structure and VSP affect the QC operations, which requires the following constraints to be considered [10]. To access an export container, the YC and its hoist are used. The YC moves along the ±X direction to reach the bay position, and then the hoist moves along the ±Y direction to reach the row position of the target container. Both the YC and its hoist can move simultaneously. After arriving above the position of the target container, the hoist then lowers in the −Z direction to pick up the container along the +Z direction. Finally, the hoist moves the container to the YC waiting in the truck lane.

Vessel Stowage Plan
A vessel stowage plan (VSP) is a plan arranging the vessel stowage positions for export containers in a vessel. Figure 3 shows the storage positions of some bays in the vessel. The assignment of stowage positions for a container is based on the container's attributes, such as size, weight, and destination. In the vessel, each bay for storing 20-foot containers is assigned with an odd bay number, with the number increasing from the bow to the stern. In addition, two 20-foot bays can form one 40-foot bay assigned with an even bay number. For instance, together Bays 09 and 10 form Bay 11. In the VSP, each container i has a position denoted as (b i , r i , t i ) [10], where the b i indicates a bay number; the r i indicates a row number; and the t i indicates a tier number.
The > 0 is set in case container i precedes container j, which violates the constraint; Otherwise = 0. The penalty costs are considered in the objective function to rule out the solutions violating the VSP constraints seriously.

Relevant Studies
Many studies have been devoted to dealing with various container terminal problems, including QCSP [12][13][14][15], YTSP [16][17][18][19], and YCSP [20][21][22][23][24][25][26][27][28][29][30]. However, individual studies tend to achieve local optima due to the lack of considering the interrelationship of relevant problems. Table 1 lists some integrated studies of some container terminal problems. Chen et al. [31] studied the QC, YT, and YC simultaneously. The authors treated the three problems as a hybrid flow shop scheduling problem (HFSP). Good coordination of these MHE was considered a key factor for the minimization of service times for ships. A mathematical model was formulated for solving these problems. However, the mathematical model becomes computationally intractable due to NP-hard. As an alternative, a Tabu search algorithm was another proposed to deal with a problem of big size. Zheng et al. [11] solved the QCSP together with the YCSP simultaneously, considering the YSP and VSP. The YSP assigns the storage locations for import containers, while the VSP assigns the The bay structure and VSP affect the QC operations, which requires the following constraints to be considered [10].

•
For some export containers to be stored at the same bay and same row, the container at a lower tier should be loaded first (i.e., a sequence constraint).

•
Containers of the same bay should be assigned to the same QC (i.e., an equipment constraint).
The first one is considered as a sequence constraint, while the second one is considered as an equipment constraint. Each violation of the sequence constraints is punished with a cost defined in the cost matrix C ij , as shown in Equation (2).
The c ij > 0 is set in case container i precedes container j, which violates the constraint; Otherwise c ij = 0. The penalty costs are considered in the objective function to rule out the solutions violating the VSP constraints seriously.

Relevant Studies
Many studies have been devoted to dealing with various container terminal problems, including QCSP [12][13][14][15], YTSP [16][17][18][19], and YCSP [20][21][22][23][24][25][26][27][28][29][30]. However, individual studies tend to achieve local optima due to the lack of considering the interrelationship of relevant problems. Table 1 lists some integrated studies of some container terminal problems. Chen et al. [31] studied the QC, YT, and YC simultaneously. The authors treated the three problems as a hybrid flow shop scheduling problem (HFSP). Good coordination of these MHE was considered a key factor for the minimization of service times for ships. A mathematical model was formulated for solving these problems. However, the math-ematical model becomes computationally intractable due to NP-hard. As an alternative, a Tabu search algorithm was another proposed to deal with a problem of big size. Zheng et al. [11] solved the QCSP together with the YCSP simultaneously, considering the YSP and VSP. The YSP assigns the storage locations for import containers, while the VSP assigns the storage position for export containers. However, this study did not consider YT. Cao et al. [32] dealt with the YCSP and YTSP simultaneously. The authors formulated a mixed-integer programming (MIP) model. However, due to computational intractability, they used Benders' decomposition method to find a solution effectively, and experiments have been conducted to investigate the effectiveness. Chen et al. [33] dealt with the QCSP, YCSP, and YTSP simultaneously by using a constraint programming (CP) model. This research considered multiple vessels with YTs shared among these vessels. The objective was to minimize the empty travel times of the QC, YC, and YT. However, this study did not take non-crossing constraints and safety margins into consideration. In addition, the CP model was found hard to solve even in a small-sized problem; thus, the authors proposed a three-stage approach, based on heuristics and disjunctive graphs, as an alternative means. This approach was found to be able to deal with a problem of up to 500 containers. Wu et al. [34] proposed a model for the integrated problem for YC and AGV in a container terminal with YSP and export container taking into consideration. The authors regarded loading as a short-term planning problem while YSP as a long-term planning problem. Furthermore, the processing sequences of containers on QCs are assumed known. A MIP model was formulated to minimize berthing times for ships. A GA was proposed to solve a big instance. Xue et al. [35] studied the integration of yard storage allocation, QCSP, and YTSP. The authors formulated a MIP for this integrated problem, which includes two weighted objectives. The first objective aims to minimize makespan, and the second objective aims to minimize the total YT travel distance. Though this study considered sequence constraints for QCs, it did not consider the non-crossing and safety distance constraints. The yard storage assignment considered block assignments instead of container slots. A two-stage heuristic algorithm was proposed to counter the computational complexity. The first stage uses an ant colony optimization (ACO) to allocate yard storage, while the second stage uses a greedy algorithm and a local search algorithm to deal with the integrated YCSP and QCSP. The authors claimed that exact approaches are impractical in solving the integrated problem. However, the generation of a lower bound remains possible. He et al. [36] studied the QCSP, YTSP, and YCSP simultaneously. The authors formulated a MIP for these problems, intending to minimize the total departure delay for all ships as well as the overall energy consumption for all tasks. An integrated simulationbased optimization approach, which includes GA and PSO, was proposed, where the GA performs the global search while the PSO performs the local search. This study illustrated the optimal trade-off between time-saving and energy-saving. Luo et al. [37] considered the integrated problem of YTSP and container storage problem. The import containers at an automated container terminal are considered. The authors formulated a MIP model for solving problems of a small size. In addition, a GA was used to deal with a big instance. However, that study did not consider QC interference and export containers. Azenodo et al. [10] solved the QCSP and the 3D stowage planning (3D SP) problem together by proposing a framework. The two problems were considered being interrelated and combinatorial. A hybrid approach combining GA with simulation was proposed as the solution means. The numbers of 30 ports, 2 QCs, and 1500 TEUs were used to test the effectiveness of this approach. According to these studies, the addition of the QCSP results in an average increase of 45.82% in load/unload time for the 3D SP problem solution. This could help the charterer avoid having to pay the ship owner for unplanned extra usage of the vessel. Kizilay et al. [38] proposed a MIP and a CP model to optimize the assignment and scheduling of QC and YC. This study treated containers as groups, called shipments, and a shipment belongs to the same port destination and the same customer. QC and YC are assigned to handle groups of containers (shipment) consecutively to simplify the problem complexity. Containers of the same shipment are stored in the same vessel bay and storage block. It was found that the CP model was more efficient than the MIP model. Yang et al. [39] studied integrated scheduling of QCs, AGVs, and YCs, as well as the routing problems of AGVs. Both import and export containers were considered. A bi-level optimization model was developed to minimize makespan. The first level considers the integrated scheduling problem, while the second level handles the AGV routing problem. A GA-based congestion-prevention rule was proposed additionally. The authors highlighted the importance of considering these problems simultaneously. Jonker et al. [40] treated the QCSP, YCSP, and AGV as a hybrid flow shop. The authors proposed a formulation for both import and export containers. In addition, this study considered job pairs, which means that a crane can handle two containers at the same time. SA is the main method used in this research. Yang et al. [41] considered the integrated scheduling problem of the QCs, ALVs, and YCs together with the storage space assignment. The objective was to improve the overall handling efficiency and accelerate shipment. This study assumed that the sequence of containers on the QCs and the yard storage locations were known. The goal was to determine the storage site for outgoing containers and to assign tasks to the various departments involved in the ALVs and YCs, such that the makespan and energy consumption is minimized. However, the assumption of a known sequence for QCs appears to be not practical. The present study is different from the above studies. Given YSP and VSP, the YCSP, YTSP, and QCSP are to be solved simultaneously based on a load-balancing concept and by using simulation-based methods.  [11] v v v v RULE Cao et al. [32] v v GA, RULE Chen et al. [33] v v v CP, RULE, DG Wu et al. [34] v v v MIP, NLMIP, GA Xue et al. [35] v v v ACO He et al. [36] v v v MIP, SIM, GA, PSO Luo et al. [37] v v v MIP, GA Azevedo et al. [10] v v SIM Kizilay et al. [38] v v MIP, CP Yang et al. [39] v v v GA Jonker et al. [40] v v v SA Yang et al. [41] v

Problem Definition
The container assignment problem, container scheduling problem, as well as VSPoriented QCSP, YTSP, and YCSP are formally defined as follows. ||R|| is the total number of resources) be a set of resources, the container assignment is a problem of assigning a resource r (r ∈ R) to handle container j (j ∈ T), denoted as (j,r). . . , ||R||] be a set of resources, the container scheduling problem is a container assignment problem (j,r), with the resource usage duration time being also specified and denoted as [S RT j,r , E RT j,r ], where the S RT j,r and E RT j,r are start and end times of using the resource instance r for container j, respectively. The container scheduling problem contains the container assignment problem implicitly. . . , ||P||] be a set of YCs, the VSP-oriented QCSP, YTSP, and YCSP is an integrated scheduling problem arranging the movements of each container j from the position (x j ,y j ,z j ) to the (b j ,r j ,t j ) by using YC, YT, and QC and considering the VSP constraints. Each usage duration of resource r is denoted as [S RT j,r , E RT j,r ], where the S RT j,r and E RT j,r are start and end usage times, respectively. The r ∈ R RT = P∪ K ∪Q. The RT indicates a resource type. If RT = 1 then r ∈ P; if RT = 2 then r ∈ K; if RT = 3 then r ∈ Q. Minimizing makespan (completion time) is the objective when fulfilling the VSP with sequence constraints taking into consideration.
To find solutions to the simultaneous YCSP, YTSP and QCSP, methodologies, such as mathematical model, heuristic, metaheuristic, and simulation, are available. Figure 4 shows the top view of a storage block that is equipped with two YCs. The operations of the hoist on YC1 are analyzed as follows.

Yard Crane (YC) Operation Analysis
where the , and , are start and end times of using the resource instance r for containe respectively. The container scheduling problem contains the container assignment problem imp itly. To find solutions to the simultaneous YCSP, YTSP and QCSP, methodologies, su as mathematical model, heuristic, metaheuristic, and simulation, are available.   Let containers i, j, and k be in sequence to be picked up by this hoist, then the routi paths 0,1,2,3, and 4 will be used sequentially. After loading container i onto the YT 2 alo path 0, the hoist will move to container j along path 1. After picking up container j, t hoist will move container j to the truck lane along path 2 and finally load it onto a YT. T YSP provides position information of containers in the block. Given ( , , ) and ( , as the storage positions of containers and j, then Equation (3) is used to estimate t total time required for the hoist to load container j onto a YT.

Yard Crane (YC) Operation Analysis
: is the moving velocity of YC along the x-direction. : is the moving velocity of YC along the y-direction. : is the moving velocity of YC along the z-direction. Let containers i, j, and k be in sequence to be picked up by this hoist, then the routing paths 0,1,2,3, and 4 will be used sequentially. After loading container i onto the YT 2 along path 0, the hoist will move to container j along path 1. After picking up container j, the hoist will move container j to the truck lane along path 2 and finally load it onto a YT. The YSP provides position information of containers in the block. Given (x i ,y,z i ) and (x j ,y j ,z j ) as the storage positions of containers i and j, then Equation (3) is used to estimate the total time required for the hoist to load container j onto a YT.
where, v x : is the moving velocity of YC along the x-direction. v y : is the moving velocity of YC along the y-direction. v z : is the moving velocity of YC along the z-direction. d(•): is a function transforming container coordinate into distance.
The right-hand side includes four-time components. The first component indicates the transportation time for the hoist to reach the bay of container j (path 1), where the d( x i − x j ) is the distance for moving from the x i to x j and the d(y j ) is the distance for reaching the row y j of container j. The second component is the roundtrip time to pick up container j, where d(H−z j ) is the moving distance. The third component is the time for moving container j to the truck lane (path 2), where d y j is the moving distance. The fourth component is the roundtrip time to load container j onto a YT, where d(H − 1) is the moving distance.
Given that S 1 j,r is the start time to process container j, then Equation (4) is used to estimate the completion time.

Yard Truck (YT) Operation Analysis
YTs move containers between QCs and YCs. Usually, a fleet of YTs are assigned to a ship, and these YTs run in a loop, one after another. The roundtrip time (PT2 j ) can be assumed as a fixed time. Minimizing the waiting time for YTs is important.
Let min j∈T {a j } be the available time of the earliest available container j (after YC operation) and min k∈K {a k } is the available time of the earliest available YT k, then Equation (5) is the earliest available start time for the YT k to transport the container.
This can minimize the waiting times for both container and YT. Given PT2 j as the roundtrip time, then Equation (6) is used to estimate the end transportation time for container j. Figure 5 shows that two QCs are working for a ship. To solve the QCSP, it needs to consider the following constraints.

•
Each QC handles one first task and one last task.

•
Each task can only be handled by one QC. • Each QC handles a sequence of containers based on their positions given by a VSP.

•
The interference of QCs should be avoided. • Given a set of containers to be loaded into a ship at the same bay and column, the container with the lowest tier should be loaded firstly.
Given as the available time of container j (j ∈ T) and as the available time of the assigned QC q (q ∈ Q), Equations (8) and (9) are used to determine the start and end usage times of the QC q, respectively.

Mathematical Model
This sub-section formulates a mathematical model for the YCSP, YTSP, and QCSP simultaneously. The assumptions, indices, sets, parameters (input data), and decision variables are first introduced as follows. Assumptions  Assume that containers i, j, and k are in a sequence to be loaded by the hoist on QC 1 into the ship. Then, the routing paths 0,1,2,3, and 4 will be used sequentially. After loading container i, the hoist will move along path 1 to access container j. The VSP provides position information of containers in the vessel. Given that (b i ,r i ,t i ) and (b j ,r j ,t j ) are the vessel storage positions of containers i and j, then the loading time for container j is estimated by Equation (7).
where, v b : is the moving velocity of QC along the bay direction. v r : is the moving velocity of QC along the row direction. v t : is the moving velocity of QC along the tier direction. H : is the height for a QC to pick up a container from a YT. d(•): is a function transforming container coordinate into distance.
The right-hand side also includes four time components. The first component is the time for the hoist to reach truck lane of bay b j , where the d b i − b j is the moving distance to reach the bay b j , and the d(r i ) is the moving distance to reach the truck lane. The second component is the roundtrip time for picking up container j from the YT 1. The third component is the moving time to reach the row position r j , where the d(r j ) is the moving distance. The fourth component is the roundtrip time to load container j into the vessel at its tier position t j , where the d(t j ) is the moving distance.
Given a j as the available time of container j (j ∈ T) and a q as the available time of the assigned QC q (q ∈ Q), Equations (8) and (9) are used to determine the start and end usage times of the QC q, respectively.

Mathematical Model
This sub-section formulates a mathematical model for the YCSP, YTSP, and QCSP simultaneously. The assumptions, indices, sets, parameters (input data), and decision variables are first introduced as follows.

•
All export containers have a storage position (x,y,z) in a block and a stowage position (b,r,t) in a vessel.

•
YCs and QCs are homogeneous in moving speed.

Parameters (Input data)
||T|| total number of tasks (containers). (X,Y,Z) dimensions of a storage block. X, Y, Z are the maximum number of bays, rows, tiers, respectively.
x j , y j , z j ) the coordinate of a container j in a block.

Decision variables
X ijp = 1, if container i is loaded before j by YC q; = 0, otherwise. Y ijk = 1, if container i is transported before j by YT k; = 0, otherwise. Z ijq = 1, if container i is loaded before j by QC q; = 0, otherwise. N p total number of containers assigned to YC p. N k total number of containers assigned to YT k. N q total number of containers assigned to QC q. PT1 ij total YC loading time of container j, estimated by Equation (3). PT3 ij total QC loading time of container j, estimated by Equation (7).
S RT j,r start usage time of resource instance r of the type RT for container j (j ∈ T). if RT = 1 then r ∈ P; if RT = 2 then r ∈ K; if RT = 3 then r ∈ Q.
E RT j,r end usage time of resource instance r of the type RT for container j (j ∈ T). if RT = 1 then r ∈ P; if RT = 2 then r ∈ K; if RT = 3 then r ∈ Q.
The mixed-integer linear programming (MILP) model of the simultaneous YCSP, YTSP, and YQSP is formulated as follows.
such that Equation (10) is the objective function minimizing the maximum completion times (makespan) of containers. Equations (1)- (9) are also included in this model. Constraints (11)- (14) relate to the decision variable X ijp . Constraint (11) requires the total number of containers assigned to YCs is the exact number ||T||. Constraint (12) ensures container i has only one successor assigned to the same YC. Constraint (13) ensures container j has only one predecessor assigned to the same YC. Constraint (14) finds the total number of containers assigned to a given YC p (i.e., N p ), and the N p should not be greater than N. Constraints (15)- (18) relate to the decision variable Y ijk . Constraint (15) stipulates that the total number of containers assigned to YTs is the exact number of ||T||. Constraint (16) ensures container i has only one successor j assigned to the same YT. Constraint (17) ensures container j has only one predecessor i assigned to the same YT. Constraint (18) finds the total number of containers assigned to a given YT k (i.e., N k ), and the N k should not be greater than ||T||. Constraints (19)- (22) relate to the decision variable Z ijq . Constraint (19) stipulates that the total number of containers assigned to YQs is the exact number of ||T||. Constraint (20) ensures container i has only one successor j assigned to the same QC. Constraint (21) ensures that container j has only one predecessor i assigned to the same QC. Constraint (22) finds the total number of containers assigned to a given QC q (i.e., N q ), and the N q should not be greater than ||T||. In order to comply with constraint (23), the loading of container j (successor) must be completed before container i (predecessor) assigned to the same YC. It is a requirement of constraint (24) that the transportation of container j (successor) does not take place before the completion of container i (predecessor) assigned to the same YT. According to constraint (25), the start time of YT carrying of a container cannot be earlier than the completion time of YC loading of the same container. This constraint links the YCSP and YTSP. Constraint (26) requires the start transportation time of YT for a container cannot before its YC completion time. This constraint links the QCSP and YTSP. Constraint (27) requires the start time of QC loading of a container cannot be before the completion time of YT transportation. Constraint (28) defines that the values of decision variables S RT j,r and E RT j,r ≥ 0. Constraint (29) defines that X ijp , Y ijk , and Z ijq are binary variables.

Simulation-Based Optimization Methods
Due to NP-hard, the MILP developed in the previous section will become computational intractable when dealing with a big instance. Thus, simulation-based methods are developed in this sub-section as the actual solution means. Figure 6 shows the framework of simulation-based optimization methods for dealing with the integrated problem of YCSP, YTSP, and QCSP. It includes the following steps.  In the simulation-based optimization methods, the load-balancing, sequencing, and simulation models are three core components. Each of them is detailed in the following sub-sections.

The Load-Balancing Heuristic
The load-balancing heuristic aims to balance the workload among equipment (QCs and YCs), considering two constraints. The first constraint requires that containers of the same bay have to be assigned to the same equipment. The second constraint stipulates that equipment cannot cross over each other. Algorithm 1 includes the logic of workload balance applicable to either QCs or YCs.  Step 1: Set parameter values for a container storage block (such as X, Y, Z, Q, H, v x , v y , v z , etc.), methods (such as ||T||, P, T, etc.), and relevant data (x j ,y j ,z j ) of YSP and (b j ,r j ,t j ) of VSP. Set penalty cost PC.
Step 2: Discover the sequence constraints required for VSP based on (b j ,r j ,t j ) data. The cost matrix C ji is generated automatically. Each violation is imposed as a penalty cost PC.
Step 3: Balance workload for available QCs and YCs by using the load-balancing heuristic. This workload-balancing strategy is expected to help achieve a good operational result.
Step 4: Form the loading sequence of containers for each YC by using different methods, such as PSO, MGPSO, GA, SBB, and the ROV technique. Alternative sequences of containers will be created. These sequences of containers will serve as inputs to Step 5.
Step 5: Simulate YC, YT, and QC operations by using the simulation model. The alternative loading sequences of containers determined in Step 4 will be used in this step.
Step 6: Evaluate the simulation results and store the best solution.
Step 7: Check the termination condition at the end of each iteration. If the termination condition is met, then go to step 8; otherwise, go to step 4.
Step 8: Output the best planning result.
In the simulation-based optimization methods, the load-balancing, sequencing, and simulation models are three core components. Each of them is detailed in the following sub-sections.

The Load-Balancing Heuristic
The load-balancing heuristic aims to balance the workload among equipment (QCs and YCs), considering two constraints. The first constraint requires that containers of the same bay have to be assigned to the same equipment. The second constraint stipulates that equipment cannot cross over each other. Algorithm 1 includes the logic of workload balance applicable to either QCs or YCs.

2:
Find the workload limit for QC/YC q by using (30).

3:
Sort containers increasingly by their bay numbers into an ascending list T; set i = 1; q = 1.

END IF 12: END For
In this algorithm, Step (1) first initializes the following parameter values: ||T|| (total number of containers), Q (total number of QCs), and W L q (workload of QC q ∈ Q). Equation (30) estimates the workload limitation for each QC in terms of the number of containers. The INT(•) is a function rounding a real value to its nearest integer. Step (7) checks whether the workload limitation has been reached? if "Yes," it considers the change of QC for assignment. However, Step (8) makes the final decision. Does this step check whether the two consecutive containers, at the i-th and (i-1)-th orders, are in the same bay? if "No" then change the assignment to the next QC; otherwise, no change. This heuristic separates QCs into different working spaces.

The Main Flow Logic of MGPSO
Heuristics and metaheuristics, such as PSO and MGPSO (an improved PSO), are used as sequencing methods, which are detailed as follows.

Particle Swarm Optimization (PSO)
Proposed by Kennedy and Eberhart [42], the PSO searches a solution space by using a group of particles. Let X i (t), V i (t), and P i (t) be the position, velocity, and personal best position vector of the particle i and X g (t) be the position of the global best particle g, then the position of the particle i at the time t+1 is determined by Equation (31).
where w is an inertia weight within the interval (0,1); c1 and c2 are acceleration coefficients usually set to the value 2.0; R1 and R2 are random numbers within the interval (0,1). The velocity V i (t) of a particle is limited within the range (V,V), where V is the maximum velocity while V is the minimum velocity allowed for a particle. Updating velocity in this way enables particles to search in the directions of their personal best and the global best. The next position of particle i is determined by Equation (32). Figure 7 shows the encoding scheme of position for particles, which includes two segments. Each segment contains N dimensions, which corresponds to the total number of containers. The load-balancing heuristic is used to find a solution for the left segment, in which each u i indicates the YC assigned to the container i. Take the vector V to become (1,2,2,1,2,1|0.2,0.4,0.1,0.5,0.3,0.9) as an example, it shows that the Containers 1, 4, and 6 are assigned to YC 1l while Containers 2, 3, and 5 are assigned to YC 2. Figure 7 shows the encoding scheme of position for particles, which include segments. Each segment contains N dimensions, which corresponds to the total nu of containers. The load-balancing heuristic is used to find a solution for the left seg in which each indicates the YC assigned to the container i. Take the vector V to b (1,2,2,1,2,1|0.2,0.4,0.1,0.5,0.3,0.9) as an example, it shows that the Containers 1, 4, are assigned to YC 1l while Containers 2, 3, and 5 are assigned to YC 2. The right segment encodes the sequences of containers to be processed by YCs. ever, the ROV technique is first used to transform the real values (i = 1,…,N), creasing order, into ranking numbers, each corresponding to the loading order o tainer i. For instance, the right segment of the vector V is transformed into (1,3,1 indicating that the sequence 1-4-6 of containers is to be processed by YC 1, while t quence 3-5-2 of containers is to be processed by YC 2. The PSO is used as a sequencing method for the simulation-based optimi method, Method3 (PSO). The right segment encodes the sequences of containers to be processed by YCs. However, the ROV technique is first used to transform the real values u i+N (i = 1, . . . , N), in increasing order, into ranking numbers, each corresponding to the loading order of container i. For instance, the right segment of the vector V is transformed into (1,3,1,2,2,3), indicating that the sequence 1-4-6 of containers is to be processed by YC 1, while the sequence 3-5-2 of containers is to be processed by YC 2.
The PSO is used as a sequencing method for the simulation-based optimization method, Method3 (PSO).

Multiple Groups Particle Swarm Optimization (MGPSO)
The MGPSO, which is an improved PSO, has the basic flow as shown in Figure 8. Each of the steps is detailed as follows.
Step 3: Evaluate and rank particles.
Step 4: Reshuffling particles into groups, which is performed at the beginning of each iteration.
Step 5: Move a particle toward the group best.
Step 6: Check whether the position of this particle is improved? if "Yes" then go to Step 12; if "No" then go to Step 7.
Step 7: Move this particle toward the global best.
Step 8: Check whether the position is better? if "Yes" then go to Step 9; if "No" then go to Step 10.
Step 9: Accept the position for this particle.
Step 10: Give a random position for this particle.
Step 11: Accept this position.
Step 12: Check whether to change to the next particle? if "Yes" then go to Step 5; if "No" then go to Step 13.
Step 13: Check whether to change to the next group? if "Yes" then go to Step 5; if "No" then go Step 14.
Step 14: Check whether to change to the next iteration? if "Yes" then go to Step 3; if "No" then go to Step 15.
Step 15: End. The MGPSO is different from the PSO with the following new features: (1) multiple groups of particles, (2) reshuffling particles at the beginning of each iteration, (3) decreasing number of groups, and (4) adaptive velocity [4].
Feature (1) enables particles to be influenced by more elites, i.e., the best particles of groups. This helps diversify particles. The particles are groped in this way. Given the number of m groups, the best particle is assigned to the 1st group; the second-best to the 2nd group; the (m)th best to the m-th group; the (m+1)th back to the 1st group, the (m+2)th back to the 2nd group, and so on. Feature (2) can intensify the diversity for particles as they can change the belonging group. Feature (3) helps implement the strategy of explorationto-exploitation as more groups enable wide exploration in a solution space while more particles in a group intensify exploitation on the best elite in each group. Equation (33) is used to determine the group numbers at each iteration t.

GN(t) =INT
Step 12: Check whether to change to the next particle? if "Yes" then go to Step 5; if "No" then go to Step 13.
Step 13: Check whether to change to the next group? if "Yes" then go to Step 5; if "No" then go Step 14.
Step 14: Check whether to change to the next iteration? if "Yes" then go to Step 3; if "No" then go to Step 15.

Set parameter values
Initialize particle positions Evaluate and rank particles Move particle toward group best  The function INT(•) rounds a real value to its nearest integer. The P is total number of particles. The T is total number of iterations. At the first iteration, t = 1, a minimum momentum is created to explore the solution space. At the last iteration, t = T, a maximum momentum is created to exploit the best elite of each group. Finally, Feature (4) enables an adaptive fly for a particle based on its current position relative to the position of a target particle (which can be the best particle of the same group or the global best particle in the swarm). Equation (34) is the adaptive velocity used by the MGPSO, which is an improvement of Equation (31) of PSO, in which the random numbers R1 and R2 tend to make random fly.Ṽ In addition, Equation (35) is used in the MGPSO to determine the next position for a particle.
The D o,j (t) is termed as total distance vector and is obtained by Equation (36).
The k-th element of the D o,j (t) is determined by Equation (37).
The ∆V j (t), which is a binary vector, is termed an adaptive velocity vector. Its k-th element is determined by Equation (38).
where the R( ) is a random number and the RB1 j (t) is a rate of binary value 1, a probability used to control the generation of the binary value 1. Both R( ) and RB1 j (t) within the interval (0,1). The RB1 j (t) is determined by Equation (39).
The D o,j (t) is termed as Hamming distance, which is obtained by Equation (40).
The Hamming distance, in fact, counts the number of different elements between two position vectors.
In Equation (41), the operator "⊕" is termed as adaptive distance operator. It works as follows.
In Equation (35), the operator " " uses the following steps: (1) it finds and adopts the first non-zero value out of theṼ j and replaces the value at the same position in the vector X j (t), (2) the replaced value takes the position of the non-zero value in the X j (t), (3) it repeats the Steps (1) and (2) until there is no non-zero value in theṼ j , and (4) the modified X j (t) becomes the next position X j (t + 1).
In addition to the four features mentioned above. The MGPSO also uses Tabu fly and Neighborhood search for moving particles. The Tabu fly aims to stop a particle from flying to the target particle o directly. This will waste one local search as this position has been visited by the target particle. Two steps are used for the Tabu fly. First, it measures the Harming distance between particles j and o. If the Harming distance ≤ 2 then it stops generating the binary value 1 for the vector ∆V j (t) to avoid a direct fly. On the other hand, the Neighborhood search aims to improve the Tabu fly, as this mechanism can keep a particle at its current position, which can also waste one local search. The Neighborhood search uses the neighborhood (p1,p2) function, which switches two elements of a position vector. It initiates a local search for the particle with a Taby fly.
Given P j (t) = (1, 2, 3, 4) and P o (t) = (4, 3, 2, 1) as positions of particles j and o at the time t, respectively, Figure 9 illustrates the movement of the particle j to its next position to approach target particle o. at its current position, which can also waste one local search. The Neighborhood se the neighborhood (p1,p2) function, which switches two elements of a position v initiates a local search for the particle with a Taby fly.
Pj (t+1) 1 Figure 9. Particle i searches around the target particle o within a group.
The following steps are used: 1. Determine the total distance vector by using Equation (34). The total distan is then derived as follows. The following steps are used: 1. Determine the total distance vector by using Equation (34). The total distance vector is then derived as follows. 2.
Determine the RB1 j (t) and D o,j (t) by using Equations (39) and (40).

4.
Determine the adaptive velocity using Equation (34). We can get theṼ j (t) as follow.

The Simulation Method
The simulation is based on a simulation model termed as Timed Predicate/Transition net (TP/T net) [43,44].

Definition 4. Timed Predicate/Transition net (TP/T net) is defined as 6-tuple:
Timed Pr/Tr net = (P,T,A,Σ,L,F) where P: a set of predicates; P = P t ∪ P nt , where P t is a set of timed predicates; P nt is a set non-timed predicates. P i ∈ P t or P nt and P t ∩ P nt = ∅.  Step 4 of the method framework. After a consuming duration at the (timed) predicate Using_R, two scenarios are available. In case that the formula (R_type < 3) in T2 is satisfied, it will trigger the transition T2, which will further transit the container token <C_ID,R_type,E_T> back to the Avail_task predicate, with the value of R_type being added 1 (meaning the change to the next operation); meanwhile the resource token <R_type,R_id,E_T> is returned to the Avail_R predicate. The E_T indicates the end usage time of the resource. In case that the logical formula (R_type = 3) in T3 is satisfied, then it will trigger the transition T3, which will move the container token <C_ID,R_type,E_T> to the Closed_task predicate, meanwhile returning the resource token <R_type,R_id,E_T> to the Avail_R predicate. The Closed_task predicate is used to keep all completed container tokens. The triggers of transitions of T1, T2, and T3 will continue until all containers are staying with the Closed_task predicate. This indicates the completion of all containers, and the simulation model then stops running. During the simulation process, the S_T (including S1, S2, and S3) and E_T of each resource usage is recorded and evaluated. The best solution can be identified and output.  Figure 10. A Predicate/Transition net model for simulation.

Numerical Experiments
Based on the method framework proposed in Section 4.1, different methods, including Method1 (SBB), Method2 (GA), Method3 (PSO), and Method4 (MGPSO), have been developed by Java programming language. Experiments were conducted in a personal notebook with an Intel PENTIUM CPU 2117U (64bits and 1.8GHz) and 4GB DRAM to investigate the effectiveness of these methods.

Parameter Values
To test the effectiveness of alternative techniques employed in stage one, further experiments are carried out in this section. Table 2 shows the parameters for YC, YT, and QC.

Experiment Data Generation
Based on the parameter values defined in Table 3, the computer generates the data of container, YSP, and VSP for experiments. . . ,<3,q,0>} is an initial marking, with the container tokens staying with the Open_task predicate while resource tokens staying with the Avail_R predicate. Each container token contains the three attributes: C_ID (container ID), R_type (resource type), and C_T (available time). Each resource token contains the three attributes: R_type (resource type), R_id (resource id), and R_T (available time). The attribute R_type can have the values 1, 2, and 3, corresponding to a YC, YT, and QC instance, respectively. The summation of labels, <C_ID,R_type,C_T> + <R_type,R_id,R_T> on A3, is a formal sum. To trigger the T1, it requires one container token and one resource token simultaneously. The "R_type < 3" in T2 and "R_type = 3" in T3 are logical formulas (conditions) for triggering the transitions T2 and T3, respectively.
The simulation model works in this way. First, the initial marking M 0 enables the transition T1. The trigger of transition T1 moves one container token <C_ID,R_type,C_T> and one resource token <R_type,R_id,S_T> to the Using_R predicate. The start resource usage time, S_T, is determined by the formula, Max{C_T, R_T}. The assignment of a container to resource instances is determined by the load-balancing heuristic specified in Step 3 of the method framework. The sequence is determined by the sequencing method specified in Step 4 of the method framework. After a consuming duration at the (timed) predicate Using_R, two scenarios are available. In case that the formula (R_type < 3) in T2 is satisfied, it will trigger the transition T2, which will further transit the container token <C_ID,R_type,E_T> back to the Avail_task predicate, with the value of R_type being added 1 (meaning the change to the next operation); meanwhile the resource token <R_type,R_id,E_T> is returned to the Avail_R predicate. The E_T indicates the end usage time of the resource. In case that the logical formula (R_type = 3) in T3 is satisfied, then it will trigger the transition T3, which will move the container token <C_ID,R_type,E_T> to the Closed_task predicate, meanwhile returning the resource token <R_type,R_id,E_T> to the Avail_R predicate. The Closed_task predicate is used to keep all completed container tokens. The triggers of transitions of T1, T2, and T3 will continue until all containers are staying with the Closed_task predicate. This indicates the completion of all containers, and the simulation model then stops running. During the simulation process, the S_T (including S1, S2, and S3) and E_T of each resource usage is recorded and evaluated. The best solution can be identified and output.

Numerical Experiments
Based on the method framework proposed in Section 4.1, different methods, including Method1 (SBB), Method2 (GA), Method3 (PSO), and Method4 (MGPSO), have been developed by Java programming language. Experiments were conducted in a personal notebook with an Intel PENTIUM CPU 2117U (64 bits and 1.8 GHz) and 4GB DRAM to investigate the effectiveness of these methods.

Parameter Values
To test the effectiveness of alternative techniques employed in stage one, further experiments are carried out in this section. Table 2 shows the parameters for YC, YT, and QC.

Experiment Data Generation
Based on the parameter values defined in Table 3, the computer generates the data of container, YSP, and VSP for experiments.
||T||: number of tasks (containers); P: total number of populations; GN(t): number of groups; n: number of particles in each group; n_ls: number of local searches; T: total number of iterations; V: low limit of velocity; Rm: mutation rate; Rc: crossover rate; and -: not used.

An Example of a Small-Sized Problem
A small-sized problem of 10 × 10 (N = 10) is used to illustrate the solution obtained from the Method4 (MGPSO). Table 4 shows the YSP data of export containers in a storage block. Row 1 indicates the container No; Row 2 shows the type of container (type = 2 means an export container); the Rows 3 to 5 show the coordinates of x, y, and z in the block, which correspond to the bay, row, and tier numbers of a container, respectively, of this container.  Table 5 shows the VSP data of export containers in a vessel. Row 1 indicates container No.; Rows 2, 3, and 4 show the bay, row, and tier numbers, respectively, for this container. Table 6 shows the cost matrix C ji automatically generated by the computer based on the VSP data (Table 5). From which, we find the VSP sequence constraints.

3.
Container 4 should be loaded before Container 7.   Table 7 shows the best solution found by Method4 (MGPSO). Row 2 indicates the YC assigned for the container i; Row 3 shows the order of container i on its assigned YC; the Rows 4 and 5, respectively, show the start and end using times on the assigned YCs; Row 6 shows the YT assigned to the container i; the Rows 7 and 8, respectively, show the start and end using time of the YT assigned; Row 9 shows the QC assigned to the container i; Row 10 shows the loading order of the QC assigned; the Rows 11 and 12, respectively, show the start and end using times the assigned QC. Table 7 shows that the best Z found by the MGPSO is 1627.5 s. The best solution is found in Iteration 2. The results are found to be a workload balanced in terms of the number of containers assigned, with constraints being also considered. For YCs, Containers 1,2,4,7, and 8 are assigned to the YC 1 with the loading sequence 2-8-4-1-7, while Containers 3,5,6,9, and 10 are assigned to YC 2 with the loading sequence 6-10-9-5-3. The S1 and E1 are start and end using times for each assigned YC. For the YTs, Containers 2,8,4,1, and 7 are assigned to YTs 1,2,4,1, and 3, respectively; while Containers 6,10,9,5, and 3 are assigned to YTs 3,5,2,4, and 5, respectively. These YTs are assigned and running in a loop; each YT serves two containers. The S2 and E2, respectively, indicate the start and end using times for each assigned YT. For QCs, Containers 1,2,4,5,6,7,8, and 10 are assigned to QC 1 with the loading sequence 2-8-6-4-10-1-7-5 into the ship, while Containers 3 and 9 are assigned to QC 2 with the loading sequence 9-3. The two loading sequences of containers are found conforming to the VSP constraints. The S3 and E3 indicate the start and end using times of each assigned QC, respectively. The assignments for QC 1 and QC 2 appear to be unbalanced. This is because of the need to conform to the constraint, i.e., containers of the same bay are assigned to the same QC. Specifically, it starts assigning Containers 2,8,4,7 (in Bay 2) and Container 10 (in Bay 3) to QC 1. However, it is found that Containers 1, 5, and 6 are all in Bay 3; thus, these containers have to be assigned to QC 1 continuously. Following this, the remaining containers, including Containers 9 and 3 (in Bay 4), are assigned to QC 2. It is noted that the end using times of QC 1 and QC 2 are 1627.5 and 1626.4, respectively, which is a desirable result. As QC 1 works in the bay space consisting of Bays 2 and 3, and QC 2 works in the bay space consisting of Bay 4, the two QCs are not inferencing each other.

Extensive Experiments
Extensive experiments have been conducted, and the results are shown in Table 8

Results Analysis and Discussions
The results are discussed accordingly, as follows.

•
The research aims to improve operational efficiency for a container terminal by coordinating YC, YT, and QC. The VSP constraints are especially respected, and they can affect the QC operation. The penalty costs of violating the VSP constraint is considered in the objective function with the purpose of exclude undesired solutions. The small-sized experiment has demonstrated the found solution, which is found to be free from violating all the VSP constraints. • As QC is bigger and heavier, it requires more setup time and tends to become a bottleneck in the container terminal operation. Therefore, coordinating YC and YT operations to facilitate the QC operation is reasonable.

•
The load-balancing concept enables better utilization of available resources. This concept has been applied to all equipment, including YC, YT, and QC, in this research.
The small-sized experiment shows that the completion of QC 1 and QC 2 in the found solution is close, which implies a good result. • This research has taken two storage positions of containers, one in the storage block and another in the vessel, into consideration. The position in the block is specified by the YSP, while the position in the vessel is specified by the VSP. Changing positions between the block and vessel require the cooperation of YC, YT, and QC; thus, coordination on these MHE is important. In this research, simulation-based optimization methods are used as the planning tool. • Method1 (SBB) is simply due to the use of one single iteration. It generates a simple solution. This method makes a YC work in one direction, from low bay to high bay number. Though it facilitates YC operation, it cannot make sure of the benefit to the QC operation. The experimental results show that Method1 is inferior to other methods in terms of makespan. • Method4 (MGPSO) is found able to find a good solution at an earlier iteration. This advantage is considered owing to its novel features, such as multiple groups, particle reshuffling, adaptive velocity, and decreasing group number. Not only to be more diversified, but particles can also use smart movements, such as Tabu fly and Neighborhood search, to better search a solution space. • Unlike the PSO, the MGPSO employs multiple groups of particles so that particles can be more diversified due to being influenced by more elites (i.e., the best particles of groups). This helps avoid particles from being trapped in local optima. • Figure 11 shows the Z value trends of the problem size 20 × 20 obtained from different methods. Method4 (MGPSO) is found to outperform Method3 (PSO), Method2 (GA), and Method1 (SBB). Method1 (SBB) finds one Z value due to its simplicity. In contrast, Method2 (GA), Method3 (PSO), and Method4 (MGPSO) have evolutionary capabilities. In addition, the evolutionary methods dive to the bottom quickly, and since then, solution improvement is found to not be significant. One likely reason is the use of the load-balancing concept, which gives a good initial foundation for developing solutions so that the contributions from the sequencing methods become non-significant. This finding suggests that a small iteration run can be employed to save the computational times required if necessary.

Conclusions
This research is devoted to solving the YCSP, YTSP, and QCSP simultaneously, with YSP and VSP data taken into account. A MILP has been formulated with the objective aimed at minimizing the makespan. However, due to NP-hard, a framework for developing simulation-based optimization methods has been proposed. Based on this framework, four different methodologies, Method1 (SBB), Method2 (GA), Method3 (PSO), and Method4 (MGPSO), have been developed. A small-sized experiment has been used to illustrate the solution found by Method4 (MGPSO). In addition, more experiments of different problem sizes, including 10 × 10, 20 × 20, 40 × 40, and 80 × 80, have been conducted. The experimental results showed that Method4 (MGPSO) outperforms the others. Statically, t-tests have been performed to test the robustness of the experiments. This research is considered with the following contributions: (1) the formulation of the MILP model; (2) the framework proposed for the development of simulation-based optimization methods; (3) the development of an improved version of standard PSO, i.e., the MGPSO; and (4) the performance of experiments to investigate these developed methods.
The following research directions can be considered: (1) in addition to the export containers, import containers can be considered simultaneously; (2) this study can be extended to include multiple storage blocks or even to multiple ships; and (3)

Conclusions
This research is devoted to solving the YCSP, YTSP, and QCSP simultaneously, with YSP and VSP data taken into account. A MILP has been formulated with the objective aimed at minimizing the makespan. However, due to NP-hard, a framework for developing simulation-based optimization methods has been proposed. Based on this framework, four different methodologies, Method1 (SBB), Method2 (GA), Method3 (PSO), and Method4 (MGPSO), have been developed. A small-sized experiment has been used to illustrate the solution found by Method4 (MGPSO). In addition, more experiments of different problem sizes, including 10 × 10, 20 × 20, 40 × 40, and 80 × 80, have been conducted. The experimental results showed that Method4 (MGPSO) outperforms the others. Statically, t-tests have been performed to test the robustness of the experiments. This research is considered with the following contributions: (1) the formulation of the MILP model; (2) the framework proposed for the development of simulation-based optimization methods; (3) the development of an improved version of standard PSO, i.e., the MGPSO; and (4) the performance of experiments to investigate these developed methods.
The following research directions can be considered: (1) in addition to the export containers, import containers can be considered simultaneously; (2) this study can be extended to include multiple storage blocks or even to multiple ships; and (3)