Integrated Berth and Crane Scheduling Problem Considering Crane Coverage in Multi-Terminal Tidal Ports under Uncertainty

: In this work, we study the integrated berth and crane scheduling problem in a tidal port with multiple terminals, considering the uncertainties, tides, maximum coverage of cranes and interference between cranes. For coping with the uncertainties, a certain number of randomly generated samples are used to evaluate the solutions, and slack variables are introduced to reduce the impact caused by the variation in vessel arrival and crane operational efﬁciency. A novel nonlinear mixed integer programming model is ﬁrst formulated for the problem to minimize the sum of expectation and variance of costs under all samples. An improved adaptive genetic algorithm, combining a simulated annealing mechanism and greedy construction strategy, is developed and implemented by MATLAB. The feasibility and validity of the algorithm and the beneﬁts of multi-terminal collaborative scheduling strategy under uncertainty are evaluated through numerical experiments. The results show that the algorithm can obtain feasible scheduling solutions with higher quality. Compared to the strategy that considers either the uncertainty or the multi-terminal collaborative mechanism, the resulting solution considering both can effectively reduce the cost and improve the competitiveness of the port.


Introduction
With the development of economic and trade globalization, container transport has been widely adopted worldwide because of its high efficiency, and the volume of container trade is also increasing annually. The enormous trade demands bring great opportunities for the container port while also raising new challenges. Figure 1 shows the flows of vessels, import containers and export containers, and presents the main operations in a port. As seen, the main resources within a container port are anchorages, berths, quay cranes, container trucks, yards and yard cranes. The operations can be classified into landside operations (including yard truck scheduling, storage space allocation and yard crane scheduling) and seaside operations (including berth allocation, quay crane assignment and quay crane scheduling).
Berths and quay cranes as the first link, and core resources play a key role to enhance the operational efficiency of ports. Therefore, this paper studies the berth allocation problem and the specific quay crane assignment problem with the following features: 1.
Container port with multiple terminals.
In practice, many ports have multiple terminals, managed by different operators, such as the Port of Valencia, Keppel Terminal, and the Tianjin Port, et al. Some terminals may also have not only one quay, but a combination of several separate quays. For these ports and terminals, as resources within them cannot be shared, they are prone to the unusual situation that some terminals are congested while some are vacant, which restricts the development of ports.

4.
Maximum coverage of quay cranes and interferences between quay cranes.
In practice, cranes cannot serve everywhere along the quay because of physical limitations. In other words, all quay cranes have their own maximum service range. In addition, when multiple cranes serve one vessel, their operational efficiency is not normally fully utilized.
All of these features have a significant impact on the scheduling solution to berths and quay cranes. However, to the best of the authors' knowledge, no study has considered these features simultaneously. This motivated us to work on this paper. In this paper, a novel mathematical model is proposed for the integrated berth and crane scheduling problem in a tidal port with multiple terminals (MBACAP), considering crane coverage under the uncertain vessel arrival time and quay crane operational efficiency. Berths and quay cranes as the first link, and core resources play a key role to enhance the operational efficiency of ports. Therefore, this paper studies the berth allocation problem and the specific quay crane assignment problem with the following features: 1. Container port with multiple terminals.
In practice, many ports have multiple terminals, managed by different operators, such as the Port of Valencia, Keppel Terminal, and the Tianjin Port, et al. Some terminals may also have not only one quay, but a combination of several separate quays. For these ports and terminals, as resources within them cannot be shared, they are prone to the unusual situation that some terminals are congested while some are vacant, which re- The rest of the paper is organized as follows: Research related to the berth allocation problem and quay crane scheduling problem is reviewed in Section 2. Section 3 presents a novel mathematical formulation and designs an improved adaptive genetic algorithm by combining a local search algorithm based on simulated annealing and a greedy construction strategy. Section 4 presents the numerical experiments for the proposed model and

Literature Review
The seaside operations in a container port involve the berth allocation problem (BAP), quay crane assignment problem (QCAP) and quay crane scheduling problem (QCSP). These problems have received significant attention in the literature. For an overview of them, Bierwirth and Meisel [1,2] and Kizilay and Eliiyi [3], who provided a comprehensive description of the relevant research, are recommended.
Research was initially focused on the separate scheduling of berths or cranes. As the research progressed, more realistic constraints were addressed and more algorithms were developed. De et al. [4] considered the cost of fuel consumption in their study for the dynamic BAP and designed a chemical reaction optimization algorithm. Liu et al. [5] developed a robust optimization model for the BAP under uncertainty by introducing an uncertainty set and designed a rolling cycle heuristic algorithm. Wang and Guo [6] modeled the vessel arrival times as a probability distribution function, so as to build an optimization model, and designed an improved genetic algorithm. Bierwirth and Meisel [7] noticed that the operational efficiency decreases when multiple cranes serve one vessel and studied the QCSP under this consideration. Al-Dhaheri and Diabat [8] considered that the distribution of containers on vessels changes during handling, which affects the stability of the vessel, and studied the QCSP under the situation. Abou Kasm and Diabat [9] studied the QCSP under the application of a new generation quay crane, capable of handling four containers in two bays simultaneously. Al-Dhaheri et al. [10] developed a genetic algorithm based on the Monte Carlo simulation for the QCSP under uncertainty.
Meanwhile, the strong correlation between BAP and QCSP makes them inseparable. The result of BAP is a prerequisite for QCSP, while the result of QCSP affects the duration of vessels in port and, thus, the result of BAP. Therefore, more research focuses on the joint scheduling of berths and cranes. Li et al. [11] integrated the BAP and QCAP, in consideration of the maximum service range of cranes. Malekahmadi et al. [12] studied the BAP and QCSP simultaneously, after considering the constraints of water depth influenced by tides. Abou Kasm et al. [13] studied the joint scheduling of berths and cranes under different assignment strategies for cranes, such as whether to allow task preemption between cranes, static or dynamic crane assignment, and analyzed the applicability of these strategies. Wang et al. [14] studied the joint scheduling problem of berths and cranes with the consideration of carbon taxation, thus, providing theoretical guidance for container port operators in a green context. Han et al. [15] jointly studied the discrete BAP and QCSP with random arrival time and random handling time, and designed a simulation-based genetic algorithm. Iris and Lam [16] studied the BAP and QCSP under uncertainty from the perspective of recoverable robustness, and designed an adaptive large neighborhood algorithm for solving it. Iris et al. [17] studied the BACAP considering both time-variant and time-invariant QC assignment policies and proposed novel Generalized Set Partitioning formulations for them. Besides, they also presented a set of column and variable reduction techniques for solving the models. They, furthermore, proposed a number of valid inequalities and variable fixing methods for enhancing the model proposed by Meisel and Bierwirth [18] and developed an adaptive large neighborhood search heuristic in [19]. Shang et al. [20] proposed a deterministic model and two robust optimization models for the integrated BACAP. In these models, the QC setup time and the handling time influenced by berth deviation and interferences of QCs were considered. They also presented a GA and an insertion heuristic algorithm to obtain near optimal solutions.
However, there are only a few studies on the collaborative scheduling of resources in ports with multiple terminals. Imai et al. [21] investigated the BAP when berths at another terminal can be shared, with two container terminals at Colombo port as an example. Despite considering two terminals, their study was still conducted mainly with one terminal and did not consider how to schedule resources, including berths and cranes, within another terminal. Similarly, Cho et al. [22] considered the BACAP by allowing the reassignment of vessels to other terminals and proposed a linear methodology model. For solving it efficiently, they adopted a sequential approach, consisting of an FBS-based heuristic, a GRASP-based heuristic, and an iterative approach.
Frojan et al. [23] studied the continuous BAP with multiple terminals and designed a heuristic algorithm based on a genetic algorithm with local search. As far as we know, they extended the single-terminal berth allocation problem to a multi-terminal situation for the first time. Nevertheless, many vital factors for the scheduling solution of berths and cranes are not considered in their study, such as uncertainty and water depth in the terminal. Therefore, in a further study, Budipriyanto et al. [24] investigated the BAP under uncertainty using the simulation approach, applying the multi-terminal collaborative mechanism. Gutierrez et al. [25] studied the static BAP and dynamic BAP for ports with two quays. They expressed the uncertain vessel arrival time as triangular fuzzy numbers and developed a fuzzy mixed-integer linear programming model and a fully fuzzy linear programming model, respectively. Venturini et al. [26] integrated the discrete BAP and the speed optimization problem for multiple ports under the emission considerations. They developed novel integer linear programming to minimize the cost of waiting, handling, delay and fuel consumption.
Krimi et al. [27] studied the multi-terminal BAP jointly with the QCAP and designed a heuristic algorithm based on the rolling strategy. However, the subject of their study was bulk ports, which makes their results unsuitable for container terminals. Grubisic et al. [28] studied the BACAP for a medium-sized terminal with a multi-quay layout and developed a mixed integer programming model to minimize the duration of vessels in port. Bouzekri et al. [29] integrated laycan and BAP and time-invariant QCAP in tidal ports with multiple quays and developed an integer programming model, after considering tides and the maximum outreach of quay cranes. Lujan et al. [30] expressed the imprecise vessel arrival time, berth and departure time as triangular fuzzy numbers and developed a fuzzy optimization model for the BACAP in port with multiple quays. The scheduling problem of vessels, trains and trucks on a multi-terminal, multi-modal maritime container port was studied by Schepler et al. [31]. They formulated an optimization model to minimize the weighted turnaround time for vessels and trains with the limitation of the inter-terminal transport of containers and proposed a restrict-and-fix heuristic based on the decomposition for solving it. Table 1 shows the comparison of studies related to seaside operations in ports with multiple terminals with our paper. " √ " indicates that the paper considered the factor indicated by the column. Based on these studies mentioned above, there is still a gap in the field of joint scheduling for berths and cranes in ports with multiple terminals under uncertainty. With this consideration, after considering the constraints of tide-influenced water depth, interference between cranes, transshipment cost within terminals and maximum service range of cranes, an optimization model is presented for the MBACAP under uncertainty.  [22] √ √ [23] √ [24] √ √ [25] √ √ [26] √ [27] √ √ [28] √ √ [29] √ √ √ √ [30] √ √ √ Tides.
The traditional way to cope with the tides is to set different time windows so as to bound the berthing and departure time of large vessels, e.g., Bouzekri, Alpan and Giard [29], Jiao et al. [32]. However, the water depth in the terminals is dynamically affected by tides. Figure 2 reveals the time windows and the actual water depth at the terminal. The red dashed line indicated the high and low tide periods. 0~T 1 and T 2~T3 are the low tide periods; T 1~T2 are T 3~T4 are the high tide periods. The blue solid line shows the actual water depth influenced by the tide. It can be found that the scheduling plan obtained by using the time window is correct only when the draft of the vessel is df2. When the draft of the vessel is df1 (less than df2), the actual reliable berthing time is t4 − t1 (greater than the high tide time T2 − T1). Therefore, the plan obtained under the time window is likely not to be optimal because it does not utilize the two periods t1~T1 and T2~t4, while when the draft is df3 (greater than df2), the actual reliable berthing time becomes t3 − t2 (less than the high tide time T2 − T1). In this case, the scheduling plan utilizes two infeasible time periods, T1~t2 and t3~T4, therefore it is likely to be infeasible. Based on the analysis above, the model refines tidal influence on the water depth to maximize the utilization of tides, i.e., each terminal has different water depths at different times. Specifically, ℎ is introduced according to the tide table to indicate the water depth of terminal m at the moment t. By comparing ℎ with the vessels' draft, we can get the reliable time for all vessels to berth.
2. Maximum coverage of quay cranes.
Since quay cranes have to occupy a certain quay segment, each crane has a maximum service range, and the cranes on the same terminal are same, yet the cranes on different terminals are different. Figure 3 depicts a terminal with four cranes. Vessel 1 can only be served by quay cranes 1 to 3. However, if the constraint of maximum service range is ignored, it is likely that quay crane 4 will be assigned to the vessel. Clearly, this is incorrect. Therefore, the model considers the coverage of cranes to ensure that all quay cranes assigned to vessels can provide the service. For coping with this constraint, we represent the starting and ending points of quay crane q located on terminal m by introducing It can be found that the scheduling plan obtained by using the time window is correct only when the draft of the vessel is df 2 . When the draft of the vessel is df 1 (less than df 2 ), the actual reliable berthing time is t 4 − t 1 (greater than the high tide time T 2 − T 1 ). Therefore, the plan obtained under the time window is likely not to be optimal because it does not utilize the two periods t 1~T1 and T 2~t4 , while when the draft is df 3 (greater than df 2 ), the actual reliable berthing time becomes t 3 − t 2 (less than the high tide time T 2 − T 1 ). In this case, the scheduling plan utilizes two infeasible time periods, T 1~t2 and t 3~T4 , therefore it is likely to be infeasible. Based on the analysis above, the model refines tidal influence on the water depth to maximize the utilization of tides, i.e., each terminal has different water depths at different times. Specifically, h t m is introduced according to the tide table to indicate the water depth of terminal m at the moment t. By comparing h t m with the vessels' draft, we can get the reliable time for all vessels to berth.

2.
Maximum coverage of quay cranes.
Since quay cranes have to occupy a certain quay segment, each crane has a maximum service range, and the cranes on the same terminal are same, yet the cranes on different terminals are different. Figure 3 depicts a terminal with four cranes. Vessel 1 can only be served by quay cranes 1 to 3. However, if the constraint of maximum service range is ignored, it is likely that quay crane 4 will be assigned to the vessel. Clearly, this is incorrect. Therefore, the model considers the coverage of cranes to ensure that all quay cranes assigned to vessels can provide the service. For coping with this constraint, we represent the starting and ending points of quay crane q located on terminal m by introducing s m q and e m q . For any quay crane q, the values of s m q and e m q are related to the number of quay cranes to its left and right. To prevent the solution of quay crane assignment from being infeasible, it must be ensured that the service intervals of all quay cranes assigned to a vessel must be able to cover the berthing interval of the vessel.
dicate the water depth of terminal m at the moment t. By comparing ℎ with the vessels' draft, we can get the reliable time for all vessels to berth.
2. Maximum coverage of quay cranes.
Since quay cranes have to occupy a certain quay segment, each crane has a maximum service range, and the cranes on the same terminal are same, yet the cranes on different terminals are different. Figure 3 depicts a terminal with four cranes. Vessel 1 can only be served by quay cranes 1 to 3. However, if the constraint of maximum service range is ignored, it is likely that quay crane 4 will be assigned to the vessel. Clearly, this is incorrect. Therefore, the model considers the coverage of cranes to ensure that all quay cranes assigned to vessels can provide the service. For coping with this constraint, we represent the starting and ending points of quay crane q located on terminal m by introducing and . For any quay crane q, the values of and are related to the number of quay cranes to its left and right. To prevent the solution of quay crane assignment from being infeasible, it must be ensured that the service intervals of all quay cranes assigned to a vessel must be able to cover the berthing interval of the vessel.

3.
Vessels already served at the beginning of the planning period.
For these terminals, vessels arrive during the planning period with different lengths, drafts, expected arrival times and expected departure times. The arrival of vessels is a continuous event, so it is rare for a terminal to be empty and without any vessels. However, few studies have considered this issue. On this basis, this paper assumes that several vessels are being served by cranes on the quay of each terminal at the beginning of the planning period.
According to Figure 4, it can be seen why we consider vessels already served at the beginning of the planning period. In the figure, vessel 1 indicates the vessel already served at the beginning of the planning period, vessel 2 indicates the vessel to be scheduled, and L is the length of the quay. Vessel 1 occupies the section of quay from x 1 to x 2 and cranes from No. 7 to No. 9. Restricted by vessel 1, vessel 2 can only berth in the section of quay from 0 to x 1 , and can only be served by cranes from No. 1 to No. 6. Therefore, although these vessels have been served at the beginning of the planning period and are not involved in the subsequent scheduling process, they cannot be ignored due to their need to occupy the quay and cranes during the planning period. 3. Vessels already served at the beginning of the planning period.
For these terminals, vessels arrive during the planning period with different lengths, drafts, expected arrival times and expected departure times. The arrival of vessels is a continuous event, so it is rare for a terminal to be empty and without any vessels. However, few studies have considered this issue. On this basis, this paper assumes that several vessels are being served by cranes on the quay of each terminal at the beginning of the planning period.
According to Figure 4, it can be seen why we consider vessels already served at the beginning of the planning period. In the figure, vessel 1 indicates the vessel already served at the beginning of the planning period, vessel 2 indicates the vessel to be scheduled, and L is the length of the quay. Vessel 1 occupies the section of quay from x1 to x2 and cranes from No. 7 to No. 9. Restricted by vessel 1, vessel 2 can only berth in the section of quay from 0 to x1, and can only be served by cranes from No. 1 to No. 6. Therefore, although these vessels have been served at the beginning of the planning period and are not involved in the subsequent scheduling process, they cannot be ignored due to their need to occupy the quay and cranes during the planning period. 4. Uncertain arrival times and operational efficiency of quay cranes for vessels.
According to Han, Lu and Xi [15], it is known that, in general, the actual vessel arrival time and the actual operation efficiency of cranes follow normal distributions. To solve these uncertainties, the berthing time of vessels is delayed by introducing buffer time ∆ and buffer efficiency ∆ , thus reducing the probability that the scheduling solution is impractical. Besides, inspired by Han, Lu and Xi [15], to evaluate the quality of the scheduling solution, a certain number of samples are randomly generated based on the distribu-

4.
Uncertain arrival times and operational efficiency of quay cranes for vessels.
According to Han, Lu and Xi [15], it is known that, in general, the actual vessel arrival time and the actual operation efficiency of cranes follow normal distributions. To solve these uncertainties, the berthing time of vessels is delayed by introducing buffer time ∆ A i and buffer efficiency ∆ v i , thus reducing the probability that the scheduling solution is impractical. Besides, inspired by Han, Lu and Xi [15], to evaluate the quality of the scheduling solution, a certain number of samples are randomly generated based on the distribution functions to simulate the actual vessel arrival time and crane operational efficiency. Finally, the scheduling solution is evaluated by the sum of the expectation and standard deviation of the total cost for all ships under all samples.

5.
Generation of feasible solutions. Figure 5 represents the scheduling process of vessels in a port with multiple terminals. When a vessel arrives, it needs to berth at a terminal with sufficient length of quay and sufficient number of cranes. If such a terminal is not available, the vessel must wait, anchored. Based on this process, we propose a heuristic method to sequentially determine the berthing time, berthing terminal, berthing position and the assigned cranes and thus to generate a feasible solution. The detailed steps can be found in Section 3.4.

Model Assumptions and Notations
The programming model presented in the paper is based on Assumptions (1)-(6). Assumption 1. The berths are continuous, which means that vessels can berth anywhere on the quay for handling containers as long as the requirements are met. Assumption 2. Different vessels have different drafts due to their container capacity, ship length and so on, while terminals also have different water depth which varies periodically due to tides. Assumption 3. The minimum number of cranes assigned to each vessel is limited by the contract between the vessel and the terminal operator, and the maximum number of cranes assigned to each vessel at the same time is limited by the length of the vessel. Assumption 4. When multiple cranes serve one vessel, the handling efficiency of each crane will decrease due to the safety distance and the cranes also have the maximum coverage range that can be reached.
Assumption 5. Due to some restrictions, the actual arrival time of vessels and the actual handling efficiency of cranes are both random variables, and they are assumed to follow the normal distribution. Assumption 6. The change in the draft caused by the container handling operation on the ship is not under consideration, i.e., the drafts of vessels are constant during the whole operation.
The indices, sets, parameters and variables in the model are defined as follows.

Model Assumptions and Notations
The programming model presented in the paper is based on Assumptions (1)-(6). Assumption 1. The berths are continuous, which means that vessels can berth anywhere on the quay for handling containers as long as the requirements are met. Assumption 2. Different vessels have different drafts due to their container capacity, ship length and so on, while terminals also have different water depth which varies periodically due to tides. Assumption 3. The minimum number of cranes assigned to each vessel is limited by the contract between the vessel and the terminal operator, and the maximum number of cranes assigned to each vessel at the same time is limited by the length of the vessel. Assumption 4. When multiple cranes serve one vessel, the handling efficiency of each crane will decrease due to the safety distance and the cranes also have the maximum coverage range that can be reached.
Assumption 5. Due to some restrictions, the actual arrival time of vessels and the actual handling efficiency of cranes are both random variables, and they are assumed to follow the normal distribution. Assumption 6. The change in the draft caused by the container handling operation on the ship is not under consideration, i.e., the drafts of vessels are constant during the whole operation.
The indices, sets, parameters and variables in the model are defined as follows.
k, m: Container terminal. i, j: Vessel to be served. t: Time period.
n: Vessel that has already been served at the beginning of planning period. q, a, b, c: Quay crane. s: Sample representing actual arrival time and operational efficiency of quay cranes for all vessels.

3.
Parameters related to terminals.

4.
Parameters related to vessels that have already been served at the beginning of planning period.
L n 0 : Length of the dummy vessel n, including safety distance, n ∈ V 0 . x n 0 : Berthing position of the dummy vessel n, n ∈ V 0 . z nm 0 : If the dummy vessel n berthed at terminal m at the beginning of planning period, z nm 0 = 1; otherwise, z nm 0 = 0, n ∈ V 0 , m ∈ T. d n 0 : Departure time of the dummy vessel n, n ∈ V 0 . θ 0 nqm : If the crane q at terminal m serves the dummy vessel n, θ 0 nqm = 1; otherwise, Parameters related to vessels to be scheduled.  Other parameters. c x : Unit deviation cost between vessel's berthing position and its desired position. c y : Operating cost per unit time for each crane. c km : Transshipment cost per unit from terminal k to m, k, m ∈ T. p i : Penalty cost per unit when the actual departure time of vessel i exceeds its expected departure time, i ∈ V. p before : Waiting cost per unit when actual arrival time earlier than scheduled berthing time.
p after : Penalty cost per unit when actual arrival time later than scheduled berthing time.
M: A large enough number.
Relaxation with the operational efficiency of cranes serving vessel i, i ∈ V. 8.
Auxiliary decision variables.
ε r ni : If the vessel i berths at the right of the vessel n, ε r ni

Objective Function
Variations in vessel arrival time and crane operation efficiency are likely to make the existing scheduling solution non-optimal, or even infeasible. Therefore, to minimize the impact of such variation on the scheduling solution, slack variables ∆ A i and ∆ v i are added to the expected vessel arrival time µ A i and the expected crane operational efficiency µ v m , respectively. Such a proactive strategy allows the scheduling solution to remain feasible and cost-neutral regardless of these variations. For evaluating the scheduling solution, a certain number of samples are randomly generated to simulate the possible vessel arrival time and crane operation efficiency. The sum of the expectation and variance in the total cost of the scheduling solution under all samples is used as the objective function of the model, as shown in Equations (1) and (2).
The lower the expected value, the lower the cost required when the scheduling solution is applied; the lower the variance, the more robust the scheduling solution is against variations of vessel arrival time and crane operation efficiency. In Equation (2), represents the operational cost of cranes. p a f ter W i U + W i D (A s i − y i ) + means the delay costs that are required when the actual arrival time is later than the scheduling berthing time. c km z ik u im W i U indicates the transshipment cost arisen by the difference between the actual terminal and the pre-assigned terminal. p be f ore W i D (y i − A s i ) + represents the waiting cost when the vessel's actual arrival time is earlier than the time in scheduling solution. p i W i U (d s i − ED i ) + indicates the delay cost when the actual departure time is later than the expected departure time. c x W i U + W i D |x i − bp i |z im u im means the deviation cost when the berthing position is different with the desired position.
To simplify our objective function and make it more efficient, the following process is to be performed for Equation (2). For For To linearize it, a new variable eu im is introduced and eu im = e i u im . Constraints (7) to (11) are also needed. For After the above transformation, Equation (2) can be replaced by Equation (16).

Constraints
In addition to the objective function, there are constraints that need to be considered in our model.
Constraint (17) indicates that all vessels must and can only berth at one terminal to load and unload containers. Besides, because of the tides, it must be ensured that the water depth of the selected terminal for any vessel is not lower than its draft during its berthing period, as shown in Constraint (18).
In Constraint (18), the t is within a variable range which makes it hard to be formulated by a solver. Hence, we introduce the σ itm and u itm to transform this constraint. σ itm is a binary constant and can be obtained after preprocessing the two constants df i and h t m . When the water depth of terminal m meets the draft of vessel i at moment t, σ itm = 1; otherwise, σ itm = 0. u itm is a binary variable. If the vessel i berths at the terminal m at the moment t, u itm = 1; otherwise, u itm = 0. To replace Constraint (18), Constraints (19) to (23) are required.
Constraint (19) ensures that vessel i can berth on the terminal m at moment t only when its draft is met. Constraint (20) implies that u itm must be 0 when the berthing terminal of vessel i is not m.
Constraints (21) and (22) ensure that u itm = 0 outside the berthing time of vessel i, i.e., t < y i or t > d i . Therefore, Constraint (23) allows u itm to be 1 when Since the quay has a length limitation, all vessels are guaranteed to berth within the quay of their berthing terminal by Constraint (24). In addition, to prevent temporal conflicts when vessels berth, that is, they occupy the same section of quay at the same time, Constraints (25) to (27) must also be satisfied. Constraint (25) is the temporal constraint for vessels to be scheduled, while Constraints (26) and (27) are the temporal constraints between the vessels to be scheduled and the vessels berthed along quay at the beginning of planning period.
In addition to the berthing terminal and berthing position, the berthing time of the vessel also has to be determined. Firstly, the vessel must berth after its arrival time. And to deal with the stochastic vessel arrival time, as shown in Constraint (28), it is imposed that the vessel must wait for an additional buffer time after its expected arrival time before berthing. The vessel's departure time is determined by the actual berthing time, the volume of containers to be handled, the operational efficiency of cranes and the number of cranes assigned. Constraint (29) defines the relationship between them. Constraints (30) and (31) define the relationship between the berthing time of any vessel and the departure time of other vessels.
Together with Constraints (25) and (30), Constraints (32) and (33) define the relationship between any two vessels to be scheduled in space and time. Combining Constraints (26), (27) and (31), Constraints (34) and (35) define the relationship between vessels to be scheduled and vessels already berthed in space and time. These constraints prevent the situation that more than one vessel occupies the same section of quay at the same time.
If the vessel operational time in Constraint (29) is considered as constant, the multiterminal BAP can be defined by Constraints (3) to (35). However, the scheduling result of cranes is related to the vessel operational time in berth allocation. Hence, the assignment of the cranes is implemented by Constraints (36) to (45). Constraint (36) ensures that the number of cranes assigned to each vessel is between its minimum number and maximum number of assignable cranes. Constraint (37) defines the relationship between the number of cranes assigned and the specific crane number and for avoiding a scenario where one crane serves two vessels at the same time, Constraints (38) and (39) are introduced. Constraint (40) keeps the numbering of cranes serving the same vessel consecutive. In other words, when any two cranes a and c serve a vessel at the same time, all cranes located between them must also serve the vessel. Constraints (41) to (43) make it impossible for a crane to cross over to provide service when more than one vessel exists in the quay at the same time, that is, when a crane a serves a vessel i, the other cranes on right of a cannot serve other vessels on the left of i. Constraints (44) and (45) ensure that the vessel berths where the crane assigned to it can serve.

The Improved Adaptive Genetic Algorithm
The MBACAP with uncertainties is more complex than the BAP which is already an NP-hard problem. Furthermore, many constraints such as the maximum service range of cranes and the water depth affected by tides need to be considered. Hence, it is hard to find the optimal solution in a polynomial time for large-scale problems. On the basis of this, an adaptive improved genetic algorithm (SA-AGA) by combining the simulated annealing mechanism and greedy construction strategy is developed. Figure 6 is the flow chart of the algorithm SA-AGA. First, we generate the initial population by combining random generation and greedy construction strategies and apply the gene repair algorithm to make the solution feasible. The children population is then produced by the combination of the elite retention strategy and the roulette wheel strategy. The search of the solution space is performed using crossover and mutation. Finally, a local search algorithm based on the simulated annealing mechanism is employed to seek a better solution.

Encoding and Decoding Rules
The encoding and decoding rules of the solution are critical for genetic algorithm. As far as the problem in our paper, representing solutions using the real number is the most effective way because of the existence of continuous variables.
A chromosome has a length of 7 × and contains a total of 7 types of information, which is represented as in Figure 7. The first three segments indicate the terminal , position , and time to berths for each vessel, respectively. The fourth segment represents the number of cranes assigned and the fifth segment is the serial number of the first quay crane . According to Figure

Population Initialization and Fitness Function
In the solution space, if the initial solution is close to the optimal solution, then few

Encoding and Decoding Rules
The encoding and decoding rules of the solution are critical for genetic algorithm. As far as the problem in our paper, representing solutions using the real number is the most effective way because of the existence of continuous variables.
A chromosome has a length of 7 × N V and contains a total of 7 types of information, which is represented as in Figure 7. The first three segments indicate the terminal T V i , position x i , and time y i to berths for each vessel, respectively. The fourth segment represents the number of cranes assigned C i and the fifth segment is the serial number of the first quay crane q start

Population Initialization and Fitness Function
In the solution space, if the initial solution is close to the optimal solution, then few iterations may be needed for the algorithm to converge on the optimal solution. Otherwise, the algorithm may require more stochasticity and longer time for iteration to close to the optimal solution. By combining the greedy construction strategy and the random generation strategy, we obtain an initial population that is both uniformly distributed in the solution space and has high quality.
Specifically, half of the initial population is generated by the greedy construction strategy. Let the berthing terminal be the pre-assigned terminal and regard the desired

Population Initialization and Fitness Function
In the solution space, if the initial solution is close to the optimal solution, then few iterations may be needed for the algorithm to converge on the optimal solution. Otherwise, the algorithm may require more stochasticity and longer time for iteration to close to the optimal solution. By combining the greedy construction strategy and the random generation strategy, we obtain an initial population that is both uniformly distributed in the solution space and has high quality.
Specifically, half of the initial population is generated by the greedy construction strategy. Let the berthing terminal be the pre-assigned terminal and regard the desired berthing position bp i as the berthing position for each vessel. Further, set the serial number of the first crane assigned to 1. The other half of the initial population are generated by the random generation strategy. The vessel's berthing terminal is randomly generated within [1, N T ] and the vessel's berthing position is randomly generated within 0, L m T − L i V . The serial number of the first crane assigned is randomly generated within 1, N m q − C i + 1 . For both strategies, randomly generate the number of cranes assigned within C min i , C max i and set the berthing time, the relaxation of arrival time, the relaxation of operational To evaluate the quality of chromosomes, the inverse of the objective function is taken as the fitness function, as shown in Equation (46). The larger the fitness value, the smaller the objective function value; in other words, the scheduling solution is more robust and less costly to apply.

Selection, Crossover and Mutation
Genetic algorithms generally search the accumulated information in depth and accomplish the iteration of the population through the selection, searching new zones of the solution space through crossover and mutation. In selection, a combination of elite retention and roulette strategy is used to determine the selection probability based on the fitness value and generate a new population. The best 10% chromosomes from the parent population are directly retained into the new population.
Algorithm 1 shows the pseudo code for the crossover and mutation. In crossover, the fusion operator based on the fitness value is applied to obtain new chromosomes. The operator considers the structure and fitness value of parent chromosomes simultaneously. In mutation, the uniform mutation strategy is used to change the berthing information in the chromosome. Furthermore, the solution converges gradually with the number of iterations. Thus, we gradually increase the mutation probability as the number of generations for which optimal solution remains unchanged increases to explore more areas in the solution space.  for j = 1 to popsize do Generate two random integers r 1 , r 2 within [1, popsize], and r 1 = r2

Gene Repair Algorithm
The population initialization can only satisfy some simple constraints such as the number of cranes and berthing position, but is powerless for complex constraints such as the maximum service range of cranes. Besides, the crossover and mutation will destroy the structure of chromosomes and make them infeasible. Therefore, the gene repair algorithm is designed to modify the infeasible solution.
Algorithm 2 shows the steps in the gene repair algorithm. Considering a chromosome, the set of scheduled vessels and the set of vessels to be scheduled are first obtained (Line 1). All vessels are then deployed in the order of their scheduled berthing times. In the gene repair algorithm, the berthing terminal, the berthing time (Lines 2~5), the berthing position (Lines 6~10), the number of cranes assigned and the serial number of cranes (Lines 11~13) are determined in turn to meet the requirements. The infeasible solutions become feasible by adjusting the berthing time, berthing position, and the crane assigned (Lines 14~16).

The Local Search Based on Simulated Annealing Mechanism
The genetic algorithm which approaches the optimal solution by starting from multiple positions in the solution space is good in global search but poor in local search and has low efficiency at the later stage of evolution. In contrast, the simulated annealing algorithm is less likely to fall into a local optimum by accepting poorer solutions with a probability. Therefore, we apply a local search algorithm based on the simulated annealing mechanism to improve the quality of solutions obtained by genetic algorithm.

4:
If ∃ h ∈ H use f ul , so that y i is not lower than the start time t start h of h, and the departure time d i is not greater than the end time t end h of h, turn to 6; otherwise, go to 5; 5: If ∃ h ∈ H use f ul , so that t end h − t start h is greater than the minimum handling time of i, go to 6; otherwise, change the berthing terminal of i and turn to 1; 6: Let the set of vessels that berthing when scheduling V now = ∅, ∀ j ∈ V m already , if d j > y i , add j to V now ; 7: Obtain the set of idle quay segments when scheduling Q idle according to the berth position x k , length L k V , the number of cranes assigned C k and the serial number of the first quay crane q start k of all vessels k in V now ; 8: Let the set of useful quay segments when scheduling Q use f ul = ∅, ∀ r ∈ Q idle , if the length of idle segment l r is not lower than L i V , and the available number of cranes NC r is not lower than C min i , add r to Q use f ul ; 9: If Q use f ul = ∅, let y i be the minimum departure time for all vessels in V now and turn to 1; 10: If ∃ s ∈ Q use f ul , so that x i ≥ x start s , and Algorithm 3 shows the steps of the local search algorithm. It is introduced to jump out of the local optimum and improve the quality of the solutions when the maximum fitness value does not change for a while. For a solution P, randomly choose a vessel and modify its berthing information including berthing terminal, berthing time, berthing position and so on to any value in their defined domain. We define the set of all solutions generated in the above way as the neighborhood Ω P of the solution.

Algorithm 3.
Steps of the local search algorithm based on simulated annealing mechanism (Source: compiled by authors).
The Local Search Algorithm Based Simulated Annealing Mechanism: Step 1: Initialize the current generation h = 1 and the system temperature T = T 0 ; Step 2: If the number of generations for which the optimal solution remains unchanged g uc ≥ 10 and h is less than the maximum iterations H max , go to Step 3; otherwise, stop the algorithm; Step 3: Obtain a new solution P new from the neighborhood Ω P of the optimal solution P max ; Step 4: Repair P new , and calculate its fitness value f new ; Step 5: If f new is greater than the maximum fitness value f max , go to Step 6; otherwise, go to Step 7; Step 6: Let P max = P new , f max = f new , g uc = 1, and stop the algorithm; Step 7: Calculate the acceptance probability p accept = e ( fnew − fmax )/T , and generate a random number r within [0, 1]. If r < p accept , go to Step 8; otherwise, go to Step 9; Step 8: Modify the worst solution of current population P min = P new and its fitness value f min = f new ; Step 9: Let T = ηT, h = h + 1, and turn to Step 2.

Numerical Experiments
To verify the algorithm, numerical experiments were carried out on a computer configured with a CPU with Intel Core i5-8300 H, 2.30 GHz and 16 GB running memory. Besides, the benefits of three different scheduling strategies were analyzed. They are the multi-terminal collaborative scheduling strategy under uncertainty (denoted as M + U), the multi-terminal independent scheduling strategy under uncertainty (denoted as S + U), and the multi-terminal collaborative scheduling strategy under certainty (denoted as M + C).

Data Design
The experimental data are necessary for running the experiments. The data involved in our experiments include the terminals data, the vessels data and the algorithm parameters.

Terminals Data Design
In our experiments, the number of terminals is three and the length of quay, the number of cranes, the operational efficiency of cranes, the interference coefficient and the cost coefficient are shown in Table 2.  Table 3 illustrates the water depth influenced by tides for each terminal during the planning period. When determining the berthing position for vessels, it must be ensured that the water depth is not less than the vessel's draft while in port; if not, the vessel may be unable to enter or leave the port.  Table 4 gives the service range of cranes for each terminal, which is related to the number of cranes on its left and right. For example, the No.4 crane of terminal 1 has three cranes on its left, so the start of its service range is 150 (50 × 3) m, while there are seven cranes on its right, so the end of its service range is 750 (1100 -50 × 7) m, where 1100 is the length of the quay, and the same goes for other cranes.

Vessels Data Design
It is essential that all terminals are free at the beginning of scheduling. Commonly, the scheduling will start along with several vessels that have berthed on the quay. Hence, the data of vessels include the data of scheduled vessels and the data of vessels to be scheduled. Table 5 indicates the range of the data for vessels, where U represents the uniform distribution and Z is the integer set. Table 5. Vessels data design.

Scheduling Solution Analysis under Three Strategies: S + U, M + C and M + U
To verify the feasibility of the SA-AGA algorithm, an instance with a scale of 40 is randomly generated. By simply adjusting the SA-AGA algorithm, we can obtain the solution for the instance under S + U and M + C scheduling strategy. In particular, keep the sixth and the seventh segments of the chromosome to 0 under the M+C scheduling strategy and keep the first segments of the chromosome that is the berthing terminal as the preassigned terminal under the S + U scheduling strategy. These algorithms are implemented via MATLAB. In these algorithms, the population size popsize is 100, the maximum number of generations G max is 500, and the initial mutational probability p 0 m is 0.6. In the local search algorithm, the initial temperature T 0 is 100, the cooling rate η is 0.8 and the maximum number of iterations H max is 100.
The scheduling solutions of the instance under S + U, M + C and M + U are shown in Figures 8-10       Similarly, we can find that the drafts of all vessels are satisfied during their berthing period, and the assigned cranes can also serve the vessels. According to Figures 9 and 10, it is observed that vessels are able to berth in any terminal, as long as the terminal has enough quays and cranes. In summary, the SA-AGA algorithm can yield feasible scheduling solutions under all three scheduling strategies.

Comparison Analysis of Algorithm before and after Improvement
Compared with the original genetic algorithm, the SA-AGA algorithm combines the greedy construction strategy and the simulated annealing mechanism into the genetic algorithm. To evaluate the performance of the SA-AGA algorithm, experiments using the algorithm without the greedy construction strategy (NG-AGA), the algorithm without the simulated annealing mechanism (NSA-AGA) and the algorithm without consideration of either (AGA) are designed under the scale of 20, 30 and 40. reveal that the objective function value can be reduced by about 50~60% by combining the simulated annealing mechanism. Comparing the curves of SA-AGA and NG-AGA, the objective function value can be reduced by about 10~20% by combining the greedy construction strategy. greedy construction strategy and the simulated annealing mechanism into the genetic algorithm. To evaluate the performance of the SA-AGA algorithm, experiments using the algorithm without the greedy construction strategy (NG-AGA), the algorithm without the simulated annealing mechanism (NSA-AGA) and the algorithm without consideration of either (AGA) are designed under the scale of 20, 30 and 40. Figure 11a-c shows the curves of the objective function values obtained by SA-AGA, NG-AGA, NSA-AGA and AGA with the iterations for the randomly generated instances with the scale of 20, 30 and 40, respectively. The curves of SA-AGA, NSA-AGA and AGA reveal that the objective function value can be reduced by about 50~60% by combining the simulated annealing mechanism. Comparing the curves of SA-AGA and NG-AGA, the objective function value can be reduced by about 10~20% by combining the greedy construction strategy.

Analysis of Experimental Results under Different Strategies and Scales
To analyze the benefits of the three scheduling strategies, S + U, M + C and M + U, 10 instances were generated randomly with the scale of 20, 30 and 40, respectively, for a total of 30 instances. For each instance, the algorithm was run five times and the average value of the solutions was used as the result, for a total of 450 runs. Table 9 illustrates the objective function values for each instance under S + U and M + U, where the Gap is the difference between them. The gap between the total costs under the two scheduling strategies is up to 35.46% (instance 3) and about 17.71% on average when the number of vessels is 20. When the number of vessels is 30, the gap is up to 27.19% (instance 14) and about 18.36% on average. When the scale is 40, the gap is up to 24.88% (instance 28) and about 13.98% on average. It is not difficult to find that the cost required for the independent scheduling of multiple terminals is generally higher than that for the collaborative scheduling of multiple terminals. The reason for the high difference in costs between the two strategies for some instances (instance 2, 3, etc.) is the large variation in the number of vessels pre-assigned to each terminal, which leads to bottlenecks at certain terminals, while by the collaborative scheduling of multiple terminals, they can equally share the workload of the whole port, thus, reducing the cost. For example, the number of vessels pre-assigned to each terminal in instance 3 is 1, 6 and 13 respectively, while it is 6, 6 and 8 for the three terminals in instance 1. Thus, the gap of instance 3 is much higher than that of instance 1.

Analysis of Experimental Results under Different Strategies and Scales
To analyze the benefits of the three scheduling strategies, S + U, M + C and M + U, 10 instances were generated randomly with the scale of 20, 30 and 40, respectively, for a total of 30 instances. For each instance, the algorithm was run five times and the average value of the solutions was used as the result, for a total of 450 runs. Table 6 illustrates the objective function values for each instance under S + U and M + U, where the Gap is the difference between them. The gap between the total costs under the two scheduling strategies is up to 35.46% (instance 3) and about 17.71% on average when the number of vessels is 20. When the number of vessels is 30, the gap is up to 27.19% (instance 14) and about 18.36% on average. When the scale is 40, the gap is up to 24.88% (instance 28) and about 13.98% on average. It is not difficult to find that the cost required for the independent scheduling of multiple terminals is generally higher than that for the collaborative scheduling of multiple terminals. The reason for the high difference in costs between the two strategies for some instances (instance 2, 3, etc.) is the large variation in the number of vessels pre-assigned to each terminal, which leads to bottlenecks at certain terminals, while by the collaborative scheduling of multiple terminals, they can equally share the workload of the whole port, thus, reducing the cost. For example, the number of vessels pre-assigned to each terminal in instance 3 is 1, 6 and 13 respectively, while it is 6, 6 and 8 for the three terminals in instance 1. Thus, the gap of instance 3 is much higher than that of instance 1.   19) and about 50.39% on average. When the scale is 40, the gap is up to 51.38% (instance 28) and about 39.18% on average. As can be seen, the total cost of the solution grows as the scale increases, regardless of which scheduling strategy is applied. When the collaborative scheduling strategy is applied, the cost of the solution obtained by considering the uncertainty is lower than when it is not considered, while the gap between the two scheduling strategies decreases gradually as the scale increases. The gap between M + C and M + U is much larger than the gap between S + U and M + U. Therefore, for the container port, the highest priority should be given to the impact of the uncertainty when scheduling, which can reduce operational costs, leading to 40% to 60% savings. Besides, the integration of all terminals and the collaborative scheduling of their internal resources are also important, which can lead to 10% to 20% savings.
uling strategy under the scale of 20, 30 and 40, respectively. For all instances, the cost of the solution is largest under M + C, the second-largest under S + U, and the smallest under M + U. The gap between M + C and M + U is much larger than the gap between S + U and M + U. Therefore, for the container port, the highest priority should be given to the impact of the uncertainty when scheduling, which can reduce operational costs, leading to 40% to 60% savings. Besides, the integration of all terminals and the collaborative scheduling of their internal resources are also important, which can lead to 10% to 20% savings.

Conclusions and Future Work
In this paper, we have proposed a novel integrated berth and crane scheduling model for tidal ports with multiple terminals. The proposed model considers the uncertain vessel arrival time, the imprecise operational efficiency of cranes, the interference between cranes and the maximum coverage of cranes. Due to the NP-hardness of the model, the adaptive improved genetic algorithm, combining a simulated annealing mechanism and greedy construction strategy, was designed for solving it.
In addition, numerical experiments, with 30 instances under the scale of 20, 30 and 40, were conducted to verify the feasibility and effectiveness of the algorithm and to analyze the benefits of three scheduling strategies: S + U, M + C and M + U. By running the algorithm in the case where the number of vessels to be scheduled is 40, the final scheduling solutions were obtained under the S + U, M + C and M + U, respectively. The solutions show that the algorithm can produce a feasible scheduling plan.
We also performed experimental comparisons of the results for the improved algorithm with the unimproved algorithm or partially improved algorithm under the same instance. The results show that the objective function value can be improved 50% to 60% after the application of the simulated annealing mechanism, while the objective function value can be reduced about 10~20% by combining the greedy construction strategy. Besides, the results also show that the cost of the scheduling solution under M + U is 10% to 20% less than that under S + U and 40% to 60% less than that under M + C.
Further research could be carried out in future works on several points. Firstly, the optimization model proposed in this paper is nonlinear, which makes it difficult to use the solver to find the optimal solution or the upper and lower bounds for the problem. With this in mind, the model can be linearized by referring to the literature [17,19] in fur-

Conclusions and Future Work
In this paper, we have proposed a novel integrated berth and crane scheduling model for tidal ports with multiple terminals. The proposed model considers the uncertain vessel arrival time, the imprecise operational efficiency of cranes, the interference between cranes and the maximum coverage of cranes. Due to the NP-hardness of the model, the adaptive improved genetic algorithm, combining a simulated annealing mechanism and greedy construction strategy, was designed for solving it.
In addition, numerical experiments, with 30 instances under the scale of 20, 30 and 40, were conducted to verify the feasibility and effectiveness of the algorithm and to analyze the benefits of three scheduling strategies: S + U, M + C and M + U. By running the algorithm in the case where the number of vessels to be scheduled is 40, the final scheduling solutions were obtained under the S + U, M + C and M + U, respectively. The solutions show that the algorithm can produce a feasible scheduling plan.
We also performed experimental comparisons of the results for the improved algorithm with the unimproved algorithm or partially improved algorithm under the same instance. The results show that the objective function value can be improved 50% to 60% after the application of the simulated annealing mechanism, while the objective function value can be reduced about 10~20% by combining the greedy construction strategy. Besides, the results also show that the cost of the scheduling solution under M + U is 10% to 20% less than that under S + U and 40% to 60% less than that under M + C.
Further research could be carried out in future works on several points. Firstly, the optimization model proposed in this paper is nonlinear, which makes it difficult to use the solver to find the optimal solution or the upper and lower bounds for the problem. With this in mind, the model can be linearized by referring to the literature [17,19] in further research, so as to make the model more widely useable and easily solvable.
Secondly, this paper assumes that a crane cannot serve other vessels across the whole operational time, once it has been assigned to a vessel, even if its tasks were all completed. Although it reduces the complexity of the problem, it also decreases the utilization of cranes. Therefore, the integrated problem of BAP and time-variant QCAP on a tidal port with multiple terminals can be considered in the next study.
In addition, other resources within the port, such as container trucks, yards and yard cranes, also affect the scheduling results of the berth and cranes. Only by integrating all the operations, the efficiency of the port can be optimized.