Heuristic/Metaheuristic-based simulation optimization approaches for integrated scheduling of yard crane, yard truck, and quay crane considering import and export containers

Yard cranes (YCs), yard trucks (YTs), and quay cranes (QCs) are the main material handling equipment (MHE) used to fulfill Yard Storage Plan (YSP) and Vessel Stowage plan (VSP) in a traditional container terminal. However, previous studies have often considered part of these MHE, with either import or export containers being neglected simultaneously. Such studies at most can only achieve local optimality which in a worst-case scenario can cause bottleneck(s) in the downstream operations. In addition, it is noted that the YSP, VSP, and their constraints have been often neglected in previous studies. Such negligence can cause operational problems. To best operate a container terminal, it needs to address the above issues. Our literature review shows that exact methods have difficulty solving a problem of big size to optimality due to NP-hard, thus alternative approaches are still required. Simulation-based optimization approaches are found with the potential to deal with container terminal problems. However, software for the integration of simulation and optimization effectively is found unavailable. In this research, different simulation-based optimization approaches, based on a simulation-based optimization framework (HSBOF), have been developed. These approaches help to schedule YCs, YTs, and QCs in handling both import and export containers for a ship in a traditional container terminal, taking YSP, VSP, and their constraints into consideration. Various heuristics/metaheuristics have been employed in the HSBOF as a sequencing method to arrange alternative operational sequences of containers. The results show that the Simulation(ISFLA) has the best performance.


I. INTRODUCTION
Recently, port congestion and supply chain disruption have introduced problems to the world, which highlights the importance of maritime transport. Maritime container transport is found especially important as it takes a large share of maritime transport. This kind of transport has an essential part of the global transportation system. The use of standard containers leads to the advantage of the ease of loading, unloading, and transportation in the supply chain. The increasing global container shipments highlight the need for continuous improvement in the capacity and efficiency of a container terminal [1].
Improving the terminal capacity, such as by increasing the total number of container terminals, is a hardware approach that can have a significant effect. However, it belongs to a long-term plan that is time-consuming and costly. To achieve an immediate effect, it needs to improve the operational efficiency of the material handling equipment (MHE) used in a container, which is considered one software approach focused in this research. In a conventional container terminal, yard cranes (YCs), yard trucks (YTs), and quay cranes (QCs) 2 are the three essential MHE for handling containers. They affect the overall performance of a container terminal considerably [3] [4]. Coordinating the operations of these MHEs is critical. Fig. 1 shows the layout of a traditional container terminal which is separated into the seaside, yard, and landside areas with different MHE [5][6][7]. For example, the seaside uses QCs to load/unload containers; the yard uses YCs to store/retrieve containers; the YTs transport containers between the QCs and YCs. These MHE are necessary to handle both import and export containers for a ship in a container terminal. The import containers are unloaded firstly, followed by the export containers. Containers can be treated as groups or individuals when planning. This research treats containers as individuals.
A container terminal faces many operational problems, such as berth allocation problem (BAP), quay crane assignment problem (QCAP), quay crane scheduling problem (QCSP), yard crane assignment problem (YCAP), yard crane scheduling problem (YCSP), yard truck assignment problem (YTAP), yard truck scheduling problem (YTSP), YT routing and traffic control problem, storage planning, space planning problem, and pre-reshuffling problem, etc. [8]. These problems are belonging to different areas. Solving these problems in one study is impossible, thus this study focuses on the YCSP, YTSP, and QCSP simultaneously.
The container terminal is a complex system that is hard to optimize its operations. This complexity comes from many factors such as different types of containers, storage positions, storage/stowage plans, and MHE. For instance, a container can have various storage positions such as one in the yard arranged by Yard Storage Plan (YSP) and another in the vessel arranged by Vessel Stowage Plan (VSP). The purpose of the activities in the container terminal can be regarded as to fulfill the two plans through the use of MHE, including YCs, YTs, and QCs.
The [9] reviewed some studies focusing on the seaside and yard side areas. That study helps identify the following research gaps: (1) the only focus on one MHE; (2) the negligence of YSP and VSP and their constraints; and (3) the lack of applying simulation approaches to deal with container terminal problems [9]. Another study [10] reviewed some maritime logistics problems and it pointed out that software for the integration of simulation and optimization effectively is unavailable. This is another research gap that is required to be addressed.
Different kinds of approaches such as exact approaches, evolutionary heuristics, heuristics, and simulation have been applied to deal with container terminal problems [9]. As kinds of exact approaches, Mixed Integer Linear Programming (MILP), Mixed Integer Programming (MIP), Dynamic Programming (DP), Branch and Bound (B&B), and Column Generation (CG) have been often used to find the optimal solution. However, these exact approaches tend to face the computationally intractable problem in dealing with a big instance so Rules, genetic algorithms (GAs), local search-based algorithms, greedy algorithms, particle swarm optimization (PSO), etc. have been used to find an optimal/near-optimal solution [9].
In this research, a heuristic/metaheuristic-based simulation optimization framework is proposed for developing simulation-based optimization approaches in which heuristics and metaheuristics serve as a sequencing method. The developed approaches are used to plan the collaborative operations of YCs, YTs, and QCs in handling both import and export containers, taking YSP and VSP and their constraints into consideration. A simulation model is used in the framework to simulate the planning results. An iterative procedure is used to improve and find the best planning result.
In this research, one heuristic (Random) and some metaheuristics, including GA, PSO, firefly algorithm (GA), shuffled frog-leaping algorithm (SFLA), and an improved SFLA (ISFLA), have been used as alternative sequencing methods respectively used in the framework. The standard SFLA, proposed by Eusuff and Lansey [11], was originally used to design the water distribution network which is not with a discrete solution domain. The standard SFLA is found with the following weaknesses; (1) lack of global strategy for grouping frogs, (2) lack of capability to handle the problem with the discrete domain, and (3) lack of adaptive jump to approach the target frog adequately. The ISFLA improves these weaknesses of the standard SFLA by introducing new features. Then, extensive experiments have been conducted to investigate their effectiveness. The experimental results show that the Simulation(ISFLA) outperforms the others in terms of the total cost consisting of sub-costs of makespan and penalty of constraint violation. The statistical t-test has been used to validate the robustness of the results.
The rest of this paper is organized as follows. Section 2 Yard Berth Rail way Bay 1 Bay 2 Bay n Landside Seaside Yard reviews background knowledge and relevant studies. Section 3 defines the problems of QCSP, YTSP, and YCSP. Section 4 proposes a hierarchical simulation-based optimization framework for developing simulation-based optimization approaches. Section 5 shows and discusses the experimental results. Finally, Section 6 concludes and suggests future research directions.

A. VESSEL STOWAGE PLAN
A container ship contains multiple bays. Each bay is composed of many stowage positions for storing containers [12]. The positions of containers in a vessel are planned by the so-called Vessel Stowage Plan (VSP). Fig. 2 shows the example of Bay 6 with stowage positions in a vessel. Each position in the VSP is denoted as ( , , ), with the b, r, and t indicating a bay, row, and tier numbers, respectively.

B. YARD STORAGE PLAN
A container yard consists of multiple blocks [13]. Fig. 3 illustrates a block of Asia type with containers stacked one on one. The assignment of containers to storage positions in a block is specified by the so-called Yard Storage Plan (YSP). Each storage position is denoted as (x,y,z), with the x, y, and z indicating bay, row, and tier number, respectively. Usually, import and export containers are respectively stored in different blocks to avoid problems. These blocks provide temporary storage service for containers before loading/unloading into/from a ship.

C. RELATED STUDIES
(1) Studies focusing on one single problem Some studies focused on QC [14]- [18]. The [14] proposed a coordinated optimization model to deal with the QCSP and QCAP (quay crane assignment problem). The Red Deer Algorithm (RDA) was at its first time applied to solve the two problems simultaneously. Differently, the [15] first decomposed the QCSP into a workload-assignment master problem and operation-sequence slave sub-problem, which were formulated as mathematical models, considering noncrossing constraints. This model was based on logic-based Bender decomposition. The. Logic-based cuts were proposed to ensure the convergence of this approach. The [16] used a two-level dynamic programming (DP) algorithm to solve the QCSP. The solution found by the DP algorithm was compared to a lower bound. The [17] dealt with the QCSP with the uncertainty of loading/unloading times of containers and travel time of quay crane taking into consideration. A simulation-optimization approach was used to find goodquality solutions. An Ant Colony Optimization (ACO) metaheuristic hybridized with a Variable Neighborhood Descent (VND) local search was used to assign tasks to QCs and determine the operational sequences of containers on each QC. The simulation was used to generate scenarios. The [18] was also devoted to solving the QCSP as quay cranes are considered the most important MHE in port container terminals. A MIP mathematical model was formulated. In addition, a B&B method was developed to find the optimal solution. To overcome the computational difficulty of the B&B method, a heuristic search algorithm called the greedy randomized adaptive search procedure (GRASP) was proposed.
Others focused on the YT [19]- [23]. The [19] formulated a mathematical model for reducing emissions generated from YTs idling. The total truck waiting time in the yard was estimated by using discrete event simulation. To solve this model, a GA was proposed and a case study in a Singaporean port was analyzed. The [20] proposed a Look-Ahead Dispatching Method for automated guided vehicles (AGVs) in Automated Port Container Terminals. This study discussed how to dispatch AGVs by utilizing information about locations and times of future tasks. The [21] studied the optimization of dispatching YTs to containers for minimizing the total service time required for a ship. The total time includes the unloading time of import containers as well as the loading time of export containers. Some easy heuristic algorithms were developed for dispatching YTs. The [22] used GA to address the YTSP. A MIP mathematical model was firstly formulated for the YTSP. However, due to NP-hard, a GA was used to solve the problem. The proposed GA has been compared to other GAs with six popular crossover operators. The [23] presented a new MIP model considering two strategies, outsourcing, and collaboration, simultaneously. In literature, the two strategies have been This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ used as alternatives to keep a large number of YTs to alleviate the possible arising of bottlenecks during peak hours. The outsourcing strategy means renting additional YTs for a container terminal while the collaborative strategy means the sharing of YTs among multiple container terminals with adjacent locations. However, the two strategies have been dealt with separately in the literature. The proposed model is solved using a rolling horizon heuristic.
Still, others focused on the YC [24]- [29]. The [24] first formulated a MIP mathematical model for the YCSP with non-interference constraints taken into consideration. However, due to NP-hard, a hybrid GA was further used to obtain near-optimal solutions. The [25] formulated a MIP for the QCSP, considering non-interference constraints. However, due to NP-hard, the authors further proposed an alternative approach. The [26] first formulated an IP mathematical model for the YCSP. The objective was to minimize the total task delay at blocks. However, due to NPhard, a novel dynamic rolling-horizon decision strategy, along with a simulation model, was then proposed for solving the YCSP. The [27] was devoted to solving the YCSP with the trade-off between efficiency and energy consumption considered. The YCSP was regarded as a vehicle routing problem with soft time windows (VRPSTW) and formulated into a MIP. The objectives aimed at minimizing the total completion delay of all task groups and the total energy consumption of all YCs. Subsequently, an integrated simulation optimization method is developed for solving the problem. The simulation is used to evaluate alternative solutions. On the other hand, an optimization algorithm that integrates the genetic algorithm (GA) and the particle swarm optimization (PSO) algorithm was used to explore alternative solutions. The [28] conducted a comparative study of two modeling approaches, centralized and decentralized, for solving the YCSP. The centralized approach was found to outperform the decentralized approach by an average of 16.5 %, due to complete and accurate information about future truck arrivals. Nevertheless, the decentralized approach can better adapt to real-time truck arrivals, making it better suited for real-life operations. Finally, it is suggested to integrate the two approaches as they have complementary features. The [29] studied the YCSP for a given set of loading/unloading jobs with different ready times. The objective was to minimize the total waiting time for jobs. A B&B algorithm was proposed to solve the scheduling problem to optimality. In addition, efficient and effective algorithms are proposed to find lower bounds and upper bounds.
(2) Integrated studies of multiple problems Other studies have been devoted to solving integrated problems. The [30] regarded the integrated scheduling problem of QC, YT, and YC as a hybrid flow shop scheduling problem. The authors highlight that good coordination among the MHE is essential to minimizing the service times of ships. In addition to formulating a mathematical model, a Tabu search (TS) was proposed for dealing with a big instance within an acceptable time. The [31] dealt with the QCSP and YCSP simultaneously with the storage positions in a yard as well as a vessel being considered. A heuristic was proposed as a scheduling means for an automated container handling system with twin 40' cranes. However, that study neglected the YTs. The [32] formulated a MIP model to deal with the YCSP and YTSP simultaneously. However, the mathematical model was found hard to solve a big instance. Two methods, based on Benders' decomposition, were thus further proposed. However, this study neglected the QCs. The [33] formulated a Constraint Programming (CP) model to deal with the QCSP, YCSP, and YTSP simultaneously. The YTs can serve multiple vessels simultaneously. The objective was to minimize empty travel times of the QC, YC, and YT. However, it is found that the CP has difficulty dealing with a big instance and even a small instance. A three-stage approach, based on heuristic and disjunctive graph, was thus proposed. This approach can deal with a big problem of up to 500 containers. However, that study neglected the VSP and YSP constraints. The [34] proposed an integrated model including YC and AGV, considering yard storage position. The objective was to optimize the yard operation. This study considered the loading operation as short-term planning and the YSP as long-term planning. In addition, it assumed that the processing sequences of containers for QCs were known. A MIP model was formulated to minimize the total berthing time for ships. However, a GA was another proposed to solve a big instance. The [35] studied an integrated problem including yard storage allocation, QCSP, and YTSP. A MIP was formulated, which minimizes two weighted objectives. The first objective is makespan and the second objective is the total transportation distance of YTs. However, this study neglected the non-crossing constraint required for QCs. The yard storage assignment considered block, instead of container slots assignments. A two-stage heuristic algorithm was proposed. The first stage uses ant colony optimization (ACO) for yard storage allocation. The second stage uses a greedy algorithm and a local search algorithm to solve the YCSP and QCSP simultaneously. The [36] studied the integrated problem of QCSP, YTSP, and YCSP. A MIP was firstly formulated to minimize the overall departure delay of ships and energy consumption of tasks. Then, a simulation-based optimization approach, which integrates GA with PSO, was proposed to solve the integrated problem. The GA performs a global search while the PSO performs a local search. The results showed the optimal tradeoff between time-saving and energy-saving. The [37] studied the YTSP together with the yard storage problem for import containers in an automated container terminal. A MIP model was formulated for dealing with small instances while a GA was further used to deal with big instances. However, this study neglected QC interference and export containers. The [12] proposed a framework to deal with the integrated problem of QCSP and 3D stowage planning problem. Simulation and GA were This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.  [38] formulated a MIP model and a CP model to optimize the assignment and scheduling of QC and YC. The containers were treated as groups (or shipments). Each group has a port destination and belongs to a customer. Then, QC and YC are assigned to handle these shipments. Containers of the same group are stored in the same vessel bay and storage block. The CP model was found more efficient than the MIP model. The [39] studied integrated scheduling of QCs, YCs, and AGVs as well as the routing problems. Both import and export containers were considered. A bi-level optimization model was proposed to minimize makespan. The top level considers the scheduling problem and the second level handles the AGV routing problem. A GA-based rule was proposed to prevent congestion. The authors highlighted the impotence to solve these problems at the same time. The [40] considered the QCSP, YCSP, and AGV simultaneously as a hybrid flow shop problem and proposed a formulation to process both import and export containers. In addition, this study allows a QC to handle two containers at the same time. SA was used as the main approach. The [41] presented a planning model for the integration of BAP, QCSP, and YTSP. The developed model contributes to the existing literature by including energy-saving goals and some realistic factors such as shortages of internal trucks and handling time estimations. A Lagrangian relaxationbased method was proposed in that study. The [42] integrated the assignment of QCs in container terminals and internal truck sharing assignments among them. A biobjective optimization model was proposed. This model includes multiple assignment phases, including the assignments of the vessel to container terminals, cranes to terminals, cranes to vessels, and trucks to cranes. The model also seeks to improve the efficiency and effectiveness of YTs by sharing them among multiple terminals. Two metaheuristic multi-objective algorithms, including modified non-dominated sorting genetic algorithm-II (MNSGA-II) and modified multi-objective particle swarm optimization (MMOPSO), are presented. The [43] presented a mathematical model for solving the QCAP, the specific QCAP, and the assignment of YTs to each QC simultaneously. The proposed model considers important practical aspects such as the limited availability and operational cost of the YTs. A Lagrangian relaxation and subgradient optimization procedure-based heuristic were proposed. The [13] developed integrated scheduling for YC and YT in the handle of export containers collaboratively in the yard side area of a container terminal. However, this study neglected QCs and the VSP constraints in the seaside area. The [44] extended the [13] to further consider the QC operations as a whole, in addition to taking the VSP constrain into consideration. However, this study and many of the past studies, including the [13], have neglected either import or export containers whose operations are found can also affect the whole operations in a container terminal. Both operations of inbound and outbound containers should be equally esteemed.
(3) Research gap The above literature review shows that many past studies have been devoted to one single MHE problem, which at most can only achieve local optimality. The local optimality in the worst case can introduce bottleneck(s) to the downstream operations if the whole operations in the container terminal are not well synchronized. Coordinating the whole operations in a container terminal is essential. A   Table 1 compares the present research with some previous studies in terms of problems domain, export/import (I/E) container, and methodology. The present research has a wider scope which includes YC, YT, QC, VSP, YSP, VSP and YSP constraints, and both import and export containers. The consideration of both import and export containers has complicated the integrated scheduling problem considerably. In addition, the present research investigates more methodologies (SIM, ISFLA, SFLA, FA, GA, PSO, and RA). Various heuristic/metaheuristic-based optimization approaches have been developed to solve this integrated scheduling problem.
The [10] reviewed studies with simulation and optimization. A total number of 107 papers published within the last two decades were reviewed and classified into three areas: terminal, shipping line, and hinterland transport. Of these papers, 20 papers were classified into the "yard" area while 8 papers were classified into the "terminal" area. This shows more simulation and optimization-based studies are required for the container terminals. Meanwhile, the [10] pointed out that software for the integration of simulation and optimization effectively was relatively unavailablewhich is regarded as one research gap in this research.

A. PROBLEM DEFINITION
The problem to be dealt with in this research is defined as an integrated scheduling problem of QC, YT, and YC in a traditional container terminal. The traditional container terminal is a non-automated container terminal that uses single-spreader-single-trolley QCs and tired-type YCs and YTs to complete the changes of storage positions of import and export containers of a ship between the yard and vessel. Each import/export container is assigned with a yard storage position and a vessel stowage position defined in the YSP and VSP, respectively. This problem also considers the YSP and VSP constraints. All of the QCs, YTs and YCs work with the single-cycle operation mode. Fig. 4 (a) shows the moving paths for the YC1 to retrieve the container j at the storage position ( , , ), after loading the container i with the storage position ( , , ) along the path 1. The hoist first moves to the container j along the path 2, then picks up this container. Subsequently, the hoist moves the container j to the truck lane along the path 3. Finally, it loads the container j onto the YT2. Eq. (1) is used to estimate the total retrieval time for the container j ( ).

B. OPERATIONS IN THE CONTAINER YARD
The right-hand side includes four terms. The first term is the moving time for the hoist to reach the position above the target container j from the position of the previous container i, where the ( − ) is a moving distance between the two bays of containers i and j; the is a moving speed along the x (bay) direction; the d( ) is the distance between the row of container j and the truck lane; the is a moving speed along the y (row) direction. The second term estimates the roundtrip time for the hoist to lower down and pick up the container j, where d(H-) is the one-trip distance; the is a moving speed along the z (tier) direction. The third term is the time required for the hoist to move the container j to the truck lane. The fourth term side is the roundtrip time for the hoist to load the container j onto the YT2, where the d(H-1) is a one-trip distance to load the container. The fourth time components form the cycle time of retrieval.
Let and be the available times of container j and its assigned YC p (p=1 in this case), respectively. Then, Eq. (2) is used to determine the start retrieval ( ) time of the YC p.
= Max , (2) The end retrieval ( ) time of the export container j is estimated by Eq. Fig. 4 (b) shows the moving paths for the YC1 to store the container j to its yard storage position ( , , ). After storing the container along the path 1, the hoist starts storing the container j. After arriving at the YT2 along the path 2, the hoist starts picking up the container j. Following this, the YC1 moves and stores the container j into its yard storage position. Eq. (4) is used to estimate the total storage time for the container j ( ).
Let and be available times of the container i and its assigned YC p, respectively. Eq. (5) is used to estimate the start storage ( ) time of the import container j for the YC p.
= Max , (5) The end storage ( ) time is estimated by Eq. (6). = + (6) When storing or retrieving containers, the sequence of these containers is important. Given a set of containers to be retrieved or stored, the following constraints should be taken into account.  YCs should not cross over each other.  Containers of the same bay are assigned to the same YC.  For some containers with the same bay and the same row number, a container with a lower-tier storage position is stored earlier and retrieved later. The last one is regarded as YSP constraints whose violation is given a penalty cost (PC) which is an additional time required for recovering this unsuitable order. Eq. (7) shows the cost matrix (M1) of storing operations. The 1 is set with the value 0 if there is no problem storing container j after container i. Otherwise, the 1 is set with a penalty cost value PC (PC > 0).
Eq. (8) shows the cost matrix ( M2 ) of container retrievals. The 2 is set with the value 0 if there is no problem retrieving container j after the container i. Otherwise, the 2 is set with the penalty cost PC (PC > 0).
The penalty costs are taken into account in the objective function to rule out those solutions violating the VSP constraints seriously.

C. OPERATIONS IN THE CONTAINER SHIP
Let ( , , ) and ( , , ) be vessel stowage positions of the containers and j, respectively. Fig. 5(a) shows the paths for the QC1 to load the container j into a ship. After loading the container into the ship, the hoist starts to pick up the container j on the YT2 in the truck lane along the path 1. After picking up the container j, the hoist starts moving the container j to its stowage position along the path 2. Finally, the container j is stored into its stowage position ( , , ). Eq. (9) estimates the total loading ( ) time for the container j.
On the right-hand side, the first term is the time for the hoist to reach above the container j, from the position of previous container i, where the ( − indicates a distance between the two bays of containers i and j; the is a moving speed along the bay direction; the d( ) indicates a distance between the bay and the truck lane; the is a moving speed along the row direction. The second term is the time for the QC1 to pick up the container j, where the CH is the one-trip distance to reach the container j on the YT2; the is the moving speed of the hoist along the tier direction.
The third term is the time for the QC1 to reach the row position of the container j. The d( ) is the moving distance to reach the row position of the container j. The fourth term is the roundtrip time to put the container j at its stowage position ( , , ), where d(H-) is the one-trip distance to reach the tier position . Let be the available time of the container j and be the available time of the assigned QC q, Eq. (10) is the start loading ( ) time of the QC q.  Fig. 5(b) shows the paths for the QC1 to unload container j from the ship. After unloading the container from the ship to the YT1 along the path 0, the hoist starts to handle the container j. After reaching above the row position of the container j along the path 1, the QC1 starts picking up the container. After picking up this container j, the QC1 starts moving this container to the YT2 in the truck lane along the path 2. Finally, the container j is loaded onto the YT2. Eq. (12) On the right-hand side, the first term is the time for the hoist moving from the container i to the container j. The second term is the pick-up time of the container. The third term is the time for the hoist moving the container j to the truck land. The fourth term is the round-trip loading time for the hoist to load the container onto the YT2.
Let be the available time of the container j and be the available time of the quay crane q assigned to the container j, Eq. (13) is the start unloading ( ) time for the QC q.
For the QC operation, the following constraints are considered [10].  QCs should not cross over each other.  Containers of the same bay are assigned to the same QC.  For some containers with the same bay and row numbers, the container with a lower-tier number is loaded earlier and unloaded later.
The last constraint is identified as a VSP constraint. Each violation of the VSP constraint is given a penalty cost. Eq. (15) shows the matrix of penalty cost (M3) for loading the container j into a ship. If the container j is loaded after the container i and this causes no problem then the 3 is set with the value 0; Otherwise, 3 is set with the penalty cost PC (PC > 0).
Eq. (15) shows the matrix of penalty cost ( M4 ) for unloading the container j from a ship. If the container j is unloaded after the container i and this causes no problem then the 4 is set with the value 0; Otherwise, 4 is set with a penalty cost PC (PC > 0).
The penalty costs will be taken into account in the objective function to rule out those solutions violating the VSP constraints seriously.

D. OPERATION OF THE YT
YTs are used to move containers between QCs and YCs. After loading an export container j onto an assigned YT, the YT will transport this container to an assigned QC at the seaside. Let be the available time of container j and Min ∈ { } be the available time of the first available YT k in the block, then the start transportation time of the export container j ( ) by the YT k is determined by Eq. (17).
Given an average one-trip transportation time for an export container (ETK), then the end transportation time ( ) of the container j is determined by Eq. (18) For an import container j, after being loaded onto a YT k, the YT then starts transporting the container j to its assigned YC in the container yard. Let be the available time of container j and Min ∈ { } be the available time of the first available YT k at the seaside, then the start transportation time of the container j ( ) by the YT k is determined by Eq. (19).
Given an average one-trip transportation time for an import container (ITK), the end transportation time of the import container j ( ) by the YT k is determined by Eq.  No reshuffling of containers.

The objective function
Eq. (21) shows the objective function used in this research, which aims to minimize the makespan as well as the cost of additional time caused by violating the YSP and VSP constraints.
, if the import container i is stored before the container j on the YC p; =0, otherwise 2 =1, if the export container i is retrieved before the container j on the YC p; =0, otherwise 1 =1, if the import container i is transported before the container j on the YT k; =0, otherwise 2 =1, if the export container i is transported before the container j on the YT k; =0, otherwise 1 =1, if the import container i is unloaded before the container j on QC q; =0, otherwise 2 =1, if the export container i is loaded before the container j on the QC q; =0, otherwise

A. A HIERARCHICAL SIMULATION-BASED OPTIMIZATION FRAMEWORK (HSBOF)
A hierarchical simulation-based optimization framework (HSBOF), shown in Fig. 6, is used to develop simulationbased optimization approaches. It includes the three main layers, (1) load-balancing, (2) sequencing, and (3) simulation, which are introduced as follows: (1) Load-balancing layer: this layer aims to balance the workload for both QCs and YCs by using a loadbalancing heuristic. Workload balance among available equipment is expected to provide a good foundation for the generation of good solutions. (2) Sequencing layer: this layer aims to form operation sequences for containers assigned to the same equipment. Heuristic/metaheuristics such as PSO, Random, SFLA, ISFLA, GA, and FA are used as the sequencing tools. The operational sequences of containers obtained in this layer will be used as inputs to the simulation layer. (3) Simulation and evaluation layer: based on the obtained operational sequences, this layer aims to simulate the YC, YT, and QC operations, meanwhile evaluating and finding the best one. Also, the HSBOF includes functions such as constraints discovery, input, and output. The constraints discovery aims to discover operational constraints from YSP and VSP. The input is a basic function to get the necessary data for experiments. The output function is a basic function to output the best planning result. A termination check is used to terminate the simulation or routing the process back to the sequencing layer again. A total number of iterations is used as the termination condition.
The three core layers are further detailed as follows.

B. THE LOAD-BALANCING HEURISTIC
This layer uses a load-balancing heuristic (Algorithm 1) to balance the workload among equipment. In Algorithm 1, Step (1) initializes parameter values for N (the total number of import and export containers), E (the total number of equipment), and (current workload of equipment e), where e ∈ P ∪ Q.
Step (2) determines the upper workload limit of each equipment. The function RD() rounds a real value up to an integer. Step (5) assigns the first container in L to the equipment e. Step (7) checks whether the upper limit (WL_limit) has been achieved?
Step (8) checks whether the next container to be assigned is located at a different bay; If "yes" then change the assignment to the next equipment by using Step (9). This algorithm ensures containers of the same bay are assigned to the same equipment. It also separates equipment into different working spaces to avoid collision(s).

C. SEQUENCING METHODS
Heuristics/Meta-heuristics are used as different sequencing methods to generate alternative simulation-based optimization approaches. The standard SFLA and an improved SFLA (ISFLA) are introduced below. 1) Standard SFLA Proposed by Eusuff and Lansey [11], the SFLA is a population-based metaheuristic mimicking the foraging behavior of frogs, which aims to find the most amount of food (the best position in a solution space). The goodness of a position can be measured by a predefined objective function. The jumps of frogs lead to search (exploration or exploitation) in the solution space. These frogs are separated into groups (termed memeplexes) and each frog can be influenced by the elites (i.e., the best frog of the same group or the global best frog). Given m groups and F frogs, the best frog is assigned to the memeplex 1; the 2 nd best to the

WL_limit=RD(N/E)
(3) SORT containers into a list L according to their bay numbers increasingly; set j=1; e=1. (4) REPEAT (5) Assign the i-th containers in L to the e; (6) = +1; // add 1 container (7) IF ( > = WL_limit) // workload limit is reached or over (8) IF ( ≠ ) //bay numbers are different (9) e=e+1; // change to the next equipment (10) END IF (11) END IF (12) UNTIL j=N This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3180752 memeplex 2; the m th best to the memeplex m; the m+1 th best back to the memeplex 1, and so on. Each group finally contains F/m frogs. However, only a specified number of q frogs (q<F/m) in a group are allowed to evolve. These frogs are sampled from the original memeplex to a sub-memeplex by using the triangular distribution shown in Eq. (22).
The R is a random number within the range [0,1]; the is the maximum leaping step allowed. With a leaping distance, the next position of frog f is determined by Eq. (24).
( + 1) = ( ) + ∆ ( ) (24) If the next position is better then accept this position; Otherwise, replace the ( ) with the ( ) in Equation (23) and jump this frog again. If this jump finds a better position then accept this position; Otherwise, give a random jump for this f. Repeat the above procedure until the last frog in the sub-memeplex is processed. Then, continue the next round of local searches for frogs. When all rounds of local searches are completed, go to the next iteration. Before starting the next iteration, a reshuffle of frogs into groups and then the sub-groups is necessary. The above procedure repeats until the satisfaction of the termination condition [44].

2) ISFLA
The improved SFLA (ISFLA) is employed by Simulation(ISFLA) to generate alternative sequences of containers.  Position Representation of a frog Fig. 7 shows the position scheme used for frogs in this research. The position scheme has a variable dimension depending on the number of containers (n) assigned to specific equipment (YC or QC). In this position scheme, each indicates the container No. assigned to the container i. Given the position vector [4,3,2,1,5,6] with n=6 as an example, it indicates that container 4 is assigned with order 1; container 3 with order 2, and so on.  The novel features of ISFLA The ISFLA includes the following features: (1) Decreasing the group number with increasing group size An exploration-to-exploitation strategy is employed by the ISFLA. Eq. (25) implements this strategy by changing the number of groups at each iteration.
The RD() is a function rounding a decimal value up to its nearest integer. The indicates the total number of frogs in the swarm, which decreases iteratively until to a bottom line. The most number of groups appears at the first iteration (t=1), giving the maximum momentum for exploration. The last iteration has the least momentum for exploration, but the maximum momentum for exploitation, as each group has the most number of forges in exploiting the target frog which can be the best in the same group or the global best frog.
(2) Discrete operators with a self-adaptive jump To make a better jump, discrete operators "~", "⊕" and "⊗" are used by the ISFLA to operate on discrete vectors.
Based on their current positions, frogs jump toward the target frog adaptively. Three steps are initiated for this. First, measure the total distance between the two frogs; Second, determine the adaptive leaping distance for the frog; Third, jump this frog to the next position for local search. The (t) is a binary vector, with its k-th element being determined by Eq. (29).
The ( ) is a random number and the 1 ( ) is a threshold controlling the generation probability of the binary value 1 for the frog f and it is determined by Eq. (30).
The HD , ( ) is termed Hamming Distance which counts the total number of different elements between two vectors and can be determined by Eq. (31). After determining ( ), the next position of the frog f is determined by Eq. (33).

) Advanced movement mechanisms for frogs
Also, the ISFLA employs two mechanisms, Tabu jump and neighborhood search, to better move a frog. The Tabu jump is used to prevent a frog f from jumping to the target frog o with one step. This will waste one local search as this position has been already visited by the target frog. This mechanism works in this way. First, it measures the , ( ) between the frogs o and f. Then, the following rule is employed.
Rule: If the , ( ) ≤ 2 then stop the frog from jumping toward the target frog o. This rule is implemented by stopping generating binary value 1 for the vector ( ). However, as Tabu jump stops jumping for a frog, the frog can also lose one local search. For improvement, the ISFLA allows such a frog to make a neighborhood search, which exchanges two randomlyselected elements in the position vector of the frog f.

D. THE SIMULATION MODEL
Referring to [45] [46], a simulation model termed Timed Predicate/Transition net (TP/T net) is developed for discreteevent simulation.
Each transition in the TP/T net model corresponds to a discrete event that can be triggered by a rule [43]. The rules used to trigger transitions T1 to T6 are detailed as follows:  predicate and the resource token will be returned to the Avail_R predicate. The following simulation procedure is used. 1) Make initial marking . Let export and import container tokens stay with the Open_task predicate and all resource tokens stay with the Avail_R predicate. 2) Trigger the enabled transitions (rules) in this model. The T1 is used to trigger operations (YC, YT, or QC) for an import container token. The R_TP specifies the operation to perform. An import container uses QC, YT, and YC sequentially. The T2 is used for export container tokens which will use YC, YT, and QC sequentially. The trigger of T1 or T2 will send one container token <C_TP, CID,R_TP,ST> and one resource token <R_TP,RID,SM, ST> to the Using_R predicate and the start event time is Max{CT,RT}, where the CT is container available time and RT is resource available time. The load-balancing heuristic and sequence method are used to determine the assignment and sequence of containers on available resources. 3) After a resource usage duration at the Using_R predicate, the transitions T3, T4, T5, or T6 will be triggered, depending on the conditions (Rules). If T3 or T4 is triggered then the container token <C_TP,CID,R_TP,ET> will go back to the Avail_task predicate and the resource token <R_TP,RID,ET> will return to the Avail_R predicate. Meanwhile, the R_TP will be decreased by 1 for an import container token to indicate the next resource type (operation) required (for an export container the R_TP is increased by 1). If Rule 5 is satisfied then it will trigger the T5 which will transition the import container token to the Closed_task predicate and return the resource token <R_TP,RID,SM,ET> to the Avail_R predicate. The T6 is used for export containers. 4) Check whether all container tokens are staying with the Closed_task predicate? When all import and export container tokens stay with the Closed_task predicate it means all operations have been completed. The STs and ETs of resource instances will be recorded during the simulation process. If "yes" go to Step 5); Otherwise, go to step 2). 5) Stop running and output the best panning result.

V. NUMERICAL EXPERIMENTS
Based on the proposed framework, we have developed six heuristic/metaheuristic-based simulation optimization approaches, namely Simulation(Random), Simulation(SFLA), Simulation(FA), Simulation(GA), Simulation(PSO), and Simulation(ISFLA), by using Java programming language. To investigate their effectiveness, experiments were conducted on a personal computer with an Intel PENTIUM CPU 2117U (64 bits and 1.8 GHz) and 4GB DRAM.  The N=20,100,200,400 as well as P=121 are common parameter values used for all approaches. The N=20 is used for a small-sized example. In addition, different approaches may have some different parameter values. For example, Simulation(Random) is set with IT=1000; Simulation(GA) is set with IT=1000, Rm=0.3, and Rc=0.4; Simulation(PSO) is set with IT=1000, w=0.5, =2, =−2; Simulation(ISFLA) is set with n=P/G(t), n_ls=4, and IT= 250.

B. EXPERIMENT DATA GENERATION
Experimental instances were generated with the parameter values defined in Table 2 considered. For example, in Table  4 the data in rows 3, 4, and 5 constitute the ( , , ) which represents the stowage position of container j in the VSP (vessel). These values are within these ranges 1≤ ≤ ‖B‖, 1≤ ≤ ‖R‖, 1≤ ≤ ‖T‖. Similarly, the data in rows 17,18, and 19 constitute the ( , , ) which represents the storage position of container j in the YSP (block). The , , are 1 M/s 1 M/s CH 30 M X: the total number of bays in a block; Y: the total number of rows in a block; Z: the total number of tiers in a block; H: the tier used to transfer a container H=Z+1; ‖P‖: the total number of YCs in a block: ‖K‖: the total number of YTs; ‖Q‖: the total number of QCs; CH: the height of QC (meters) to reach YT; ‖B‖: the total number of bays in a vessel: ‖R‖: the total number of rows in a bay; ‖T‖: the number of tiers in a vessel.   Table 2.
The position data of each container is generated by the computer automatically. For an import container, the container type is set to 1 (C_TP=1). In Tables 4 and 5, these marked rows contain the data generated by the computer automatically, while other unmarked rows 6-16 are those relating to decision variables whose values are determined by the approach used. Table 5 shows the data of export containers. The rows 3,4 and 5 constitute the current storage position ( , , ) in the YSP while the rows 17, 18, and 19 represent the future stowage position ( , , ) in the VSP for the container j to be transported by the ship.

C. A SMALL-SIZED EXAMPLE
This subsection illustrates the inputs and outputs of a smallsized problem (N=20) obtained from the Simulation(ISFLA).

2) VSP and YSP constraints
Based on the input VSP data, the VSP constraints for both import and export containers and the YSP constraints for both import and output containers are generated automatically. These constraint data are shown in Appendix A and B with PC=600 seconds used as the penalty cost. Appendix A is the cost matrix M4 of VSP constraints for import containers. It shows no constraint is imposed.
Appendix B shows the cost matrix M1 of YSP constraints for import containers, which includes the following constraints:  Container 2 should be unloaded before container 3;  Container 16 should be loaded before container 10.
Appendix C shows the cost matrix M2 of YSP constraints for export containers, which suggests the following constraint:  Container 20 should be unloaded before container 15.
Appendix D shows the cost matrix M3 of YSP constraints for export containers, which includes the following constraint:  Container 9 should be loaded before container 13.  For each violation of the YSP or VSP constraints, a penalty cost of 600 is imposed. Table 4 shows the best solution found by Simulation(ISFLA) in the process of import containers. Row 6 shows the QC assigned to the container j; row 7 shows the unloading order of container j on its assigned QC; rows 8 and 9 respectively show the start and end completion time of the QC; row 10 shows the YT assigned to the container j; rows 11 and 12 respectively show the start and completion times of the YT; row 13 shows the YC assigned to the container j; row 14 shows the storage order of the container j on the assigned YC; rows 15 and 16 show the start and completion time of the YC, respectively. The best Z is 2456.9 (s) found by the Simulation(ISFLA).

3) Output result
The import containers 11,7,9,14,17,19,5,8,4, and 6 are assigned to the QC1 with this unloading sequence while the import containers 2,18,15,3,13,1,20,16,10, and 12 are assigned to QC2 with this unloading sequence. The two unloading sequences are feasible due to no violation of any of the VSP constraints. In addition, the two QCs are workload balanced. Following the QC operations, the YTs are subsequently used to transport these containers to their assigned YCs. The workload of these YTs is found to be balanced, each transporting 4 containers. At the container yard, the import containers 7,18,9,13,19,1,5,20,4,6 and 12 are assigned to the YC1 and loaded with this sequence while the import containers 11,2,15,14,3,17,8,16 and 10 are assigned to YC2 and loaded with this sequence. The two YCs are workload balanced and subject to the YSP constraints. Table 5 shows the best solution found by Simulation(ISFLA) in the handle of export containers. The best Z=5397.4 (s) is found at the 245th iteration, at the computational cost of 6.1 seconds. This table shows that the export containers 8,17,11,6,12,20,3,15,13,19 and 16 are assigned to the YC1 with this retrieval sequence while the export containers 18,5,14,7,9,1,2,4 and 10 are assigned to YC2 with this retrieval sequence. The two sequences are conforming to YSP constraints and the two YCs are almost balanced and conform to the requirement that containers of the same bay are assigned to the same YC. Also, the YTs are workload balanced, each transporting 4 containers. The YTs run in a loop and take turns transporting the import containers. For the QCs, the export containers 17,11,6,9,20,1,15,13,10, and 16 are assigned to QC1 with this loading sequence while the export containers 8,18,5,14,7,12,3,2,4, and 19 are assigned to QC2 with this loading sequence. Both QC1 and QC2 are workload balanced, each loading 10 containers. Table 6 shows the data of the best solutions found by different approaches. It is found that the Simulation(ISFLA) approach finds the best solution, but with the most computational cost.

E. EXTENSIVE EXPERIMENTS
Extensive experiments of different problem sizes have been further conducted to investigate the effectiveness of different approaches. Table 7 shows the results obtained from these different approaches. They are summarized as follows:  Fig. 10

F. STATISTICALLY TEST
The statistical t-test is used to test the Hypotheses H0 and H1 at the significance level of =0.05. The Ho is a null hypothesis assuming that the average Z values obtained from comparing approaches have no significant difference; H1 assumes that these average Z values are different. Table 8 shows test results of comparing each approach to Simulation(ISFLA). The symbol "+" that Simulation(ISFLA) is better; the symbol "-" indicates that Simulation(ISFLA) is inferior; the symbol "N" indicates no difference. It shows the p-values of each comparison are less than 0.005. For example, at the problem size N=100, the comparisons of Simulation(ISFLA) with Simulation(Random), Simulation(FA), Simulation(GA), and Simulation(SFLA), Simulation(PSO) have the following p-values, 1.34839E-6, 3.91071E-8, 4.84051E-5, 0.0002207, 5.40429E-7, respectively, which lead to rejection of H0 and the acceptance of H1. Simulation(ISFLA) is concluded to outperform the others.

G. ANALYSIS AND DISCUSSIONS
 While balancing workload, the load-balancing heuristic meanwhile can separate YCs and QCs into different working spaces. This avoids crossover of each other. Subsequently, a heuristic/metaheuristic is used to initiate alternative operational sequences for containers assigned on each equipment. Finally, the simulation model is used to simulate and evaluate these sequences and identify the best one.  For the small-size experimental instance, it shows that the solution found by Simulation(ISFLA) has a good quality and is free from the violation of the VSP and YSP constraints. Specifically, the completion times of import containers for the YC1 and YC2 are 2456.9 and 2455.8, respectively. It shows that import containers are balanced and parallel handled by the two YCs. The completion times of export containers for the QC1 and QC2 are 5390.2 and 5397.4, respectively, which shows that the two QCs are also with balanced workload. In addition, the YTs are also with a balanced workload and well utilized with no spare time.  Fig. 13 shows the experimental results of sensitivity analysis using different numbers of YTs, together with 2 QCs and 2 YCs in handling 100 import and 100 export containers for a ship. The use of one YT leads to a huge cost (Z=119662.5) that can be further reduced by increasing the number of YTs from 2 to 8. However, the numbers of YTs exceeding 8 are found unable to reduce the Z value (19442.8) further, meaning that the best number of YTs to cooperate with the 2 QCs and 2 YCs for this ship is 8. If the number of YTs is below 8 then improvement is available by increasing this number.  In our experiments, Simulation(ISFLA) is found better than the other approaches in terms of makespan. The advantage is believed to result from the novel features of the ISFLA, such as decreasing the number of frog groups, self-adaptive leaps, and the advanced movement mechanism. In addition, the ISFLA is with more population advantage due to allowing all frogs to evolve [47].  Containers can be treated as groups or individuals when planning [35] [44]. Though the treatment of containers as groups makes the planning easier, it cannot facilitate real operations due to rough resolution. By treating containers as individuals, our approach can offer a high-resolution solution that makes implementation easier.  Effective software for the integration of simulation and optimization effectively is still unavailable [10], which is considered one research gap in the existing literature. To address this gap, we have proposed a hierarchical simulation-based optimization framework (HSBOF) for developing alternative heuristic/metaheuristic and simulation-based optimization approaches to deal with the integrated scheduling problem of QCs, YTs, and YCs. This helps container terminal planners to pursue global optimality, instead of local optimality.  Managers should note that local optimality is not certainly beneficial as it may form bottleneck(s) in the downstream operations. It is important to consider all MHE as well as both import and export containers. To achieve global optimality, the simulation-based optimization approach is useful. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. This research proposed a hierarchical simulation-based optimization framework (HSBOF) for developing heuristic/metaheuristic-based simulation optimization approaches to deal with the simultaneous YCSP, YTSP, and QCSP, taking YSP, VSP, and their constraints into consideration.

VI. CONCLUSION
Specifically, Simulation(Random), Simulation(ISFLA), Simulation(PSO), Simulation(SFLA), Simulation(GA), and Simulation(FA) have been developed in this research. The objective is to minimize costs of makespan as well as penalty costs of constraint violation. A small-size instance has been used to illustrate the planning results, and extensive experiments of different sizes (N=100,200, and 400) have been conducted to investigate the effectiveness of different approaches. Finally, the statistical t-tests have been used to validate the experimental results. The Simulation(ISFLA) is found to outperform the others.
The contributions of this research are listed as follows: (1) we have proposed a hierarchical simulation-based optimization framework (HSBOF) for the integration of simulation and optimization. Based on this framework, various heuristic/metaheuristic-based simulation optimization approaches have been developed; (2) these approaches have been successfully applied to solve the YCSP, YTSP, and QCSP simultaneously, taking the YSP and VSP and their constraints into consideration; (3)  He has published 80 journal papers and 40 patents. He achieves around 25 projects of research and academy and has some successful cases of technology transfer. He also achieves awards from many international invention competitions (2 Gold Medal Awards, 5 Silver Medal Awards, and several Bronze Medal Awards). In addition, he has plenty of practical experience in many industries and provides consultant services for several stock companies in Taiwan. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3180752