Robust Scheduling for Berth Allocation and Quay Crane Assignment Problem

Decision makers must face the dynamism and uncertainty of real-world environments when they need to solve the scheduling problems. Different incidences or breakdowns, for example, initial data could change or some resources could become unavailable, may eventually cause the infeasibility of the obtained schedule. To overcome this issue, a robust model and a proactive approach are presented for scheduling problems without any previous knowledge about incidences. This paper is based on proportionally distributing operational buffers among the tasks. In this paper, we consider the berth allocation problem and the quay crane assignment problem as a representative example of scheduling problems. The dynamism and uncertainty are managed by assessing the robustness of the schedules. The robustness is introduced by means of operational buffer times to absorb those unknown incidences or breakdowns. Therefore, this problem becomes a multiobjective combinatorial optimization problem that aims to minimize the total service time, to maximize the buffer times, and to minimize the standard deviation of the buffer times. To this end, a mathematical model and a new hybrid multiobjective metaheuristic is presented and compared with two well-known multiobjective genetic algorithms: NSGAII and SPEA2


Introduction
Within a container terminal, operations related to move containers can be divided into four different subsystems (shipto-shore, transfer, storage, and delivery/receipt) [1]. In each subsystem, terminal operators must deal with with different complex optimization problems that can be overcome by using artificial intelligence techniques. For instance, berthing allocation or stowage planning problems are related to the ship-to-shore area [2][3][4][5], remarshalling problem and transport optimization [6] to the storage and transfer subsystems, respectively, and planning and scheduling hinterland operations related to trains and trucks to the delivery/receipt subsystem [7].
In this paper, we focus on two problems related to the ship-to-shore area, the berth allocation problem (BAP) and the quay crane assignment problem (QCAP). The former is a well-known combinatorial optimization problem [8], which consists in assigning berthing positions and mooring times to incoming vessels. The QCAP deals with assigning a certain number of quay cranes (QCs) to each moored vessel such that all required movements of containers can be fulfilled [9].
A comprehensive survey of BAP and QCAP is given in [9]. These problems have been mostly considered separately, with an interest mainly focused on BAP. An interesting approach for BAP is presented by Kim and Moon [10] where a simulated annealing metaheuristic is compared with a mathematical model. However, there are some studies on the combined BAP + QCAP considering different characteristics of berths and cranes [11][12][13][14][15].
Most of the research in scheduling has been focused on deterministic and complete information, but they are usually not satisfied in real-world environments. Due to the fact that the real world is uncertain, imprecise, and nondeterministic, there might be unknown information, breakdowns, incidences, or changes, which make the initial plans or the obtained schedules become invalid. Thus, there are new trends to cope these aspects in the optimization techniques: proactive and reactive approaches [16]. In this paper, a proactive approach is studied within the berth allocation and 2 Mathematical Problems in Engineering the quay crane assignment problems. The uncertainty within these problems is due to lower movements per time unit than expected or engine failures in quay cranes, among others. Due to the introduction of this new objective in the scheduling optimization problem, a multiobjective optimization approach needs to be taken into consideration.
All the above studies do not take into consideration the uncertainty of the real world to obtain a robust scheduling. Robustness is a measure of the performance characterization of an algorithm in the presence of uncertainties [17]. However, there are some studies that address the robust scheduling. In [18], a robust optimization model for cyclic berthing for a continuous and dynamic BAP is studied by minimizing the maximal crane capacity over different arrival scenarios of a bounded uncertainty given by their arrival agreements. In [19], a proactive approach for a discrete and dynamic model of the BAP is presented taking into account uncertainties in the arrival and handling times given their probability density functions. They propose a mixed integer programming model and a genetic algorithm (GA) for both problems: discrete berth allocation and QC assignment. The objective is to minimize the sum of expected value, the standard deviation of the service time, and the tardiness of the incoming vessels.
Robust scheduling based on operational buffers has already been introduced as a proactive approach in the BAP. An approach to robust BAP is presented in [20]. They presented a feedback procedure for the BAP that iteratively improves the robustness of the initial schedule. This feedback procedure determines the time buffers for each vessel by means of adjustment rules.
In [21], another approach to the robust BAP is solved by a scheduling algorithm that integrates simulated annealing and branch-and-bound algorithms. This study introduces the robustness as an objective to be maximized and an evaluation is carried out by varying the weights of these functions. The robustness is achieved by a constant buffer time assigned to all vessels.
In [22], the robust BAP problem is studied as a proactive strategy as a multiobjective optimization problem. They solved this problem with the squeaky wheel optimization (SWO) metaheuristic. The first objective is to minimize the late departures and the deviation from the desired position; the second objective is to maximize the robustness of the schedule. They tackle the robustness measure as a diminishing return, specifically the exponential function, to capture the decreasing marginal productivity of slacks in a berthing schedule.
However, most of the above approaches consider discrete berths or previous knowledge about the uncertainty in arrival or handling times to produce robust schedules, but usually this knowledge is not available. Furthermore, other approaches propose how to obtain robust schedules by means of operational buffer times, but these buffers are set independently of the handling (or processing) time of the vessels.
Overcoming the above approaches, hybrid metaheuristics for both single and multiobjective combinatorial optimization problems have received a significant interest from the research community [23,24], and also they have been used in a wide range of real-world applications [25].
In this paper, we introduce a robust model to deal with limited incidences with no previous knowledge about them (Section 3) as well as a multiobjective approach to face this problem (Section 4). A formal mixed integer programming (MIP) is presented for the dynamic and continuous robust BAP + QCAP that extends the model presented in [10] (Section 5). Section 6 presents our proposed hybrid multiobjective genetic algorithm based on the scheme NSGAII [26] in order to obtain near-optimal solutions in an efficient way. This hybrid algorithm is used to solve the BAP + QCAP with a continuous quay and dynamic arrivals as well as to provide robust solutions by using operational buffers. As there is no previous knowledge about the incidences, these operational buffers are proportionally distributed among the tasks to be able to absorb as many incidences as possible. Thereby, a new objective function (standard deviation of robustness measures) was introduced to pursue this goal. This algorithm is compared with the mathematical model presented and two well-known multiobjective genetic algorithms: NSGAII and SPEA2+ [27] (Section 7). The development of the technique presented in this paper will provide the terminal operators with different robust berthing plans which are able to absorb limited incidences.
The overall collaboration goal of our group at the Universitat Politècnica de València (UPV) with the Valencia Port Foundation and the maritime container terminal MSC (Mediterranean Shipping Company S.A.) is to offer assistance and help in the planning and scheduling of tasks such as the allocation of spaces to outbound containers, to identify bottlenecks, to determine the consequences of changes, to provide support in the resolution of incidents, and to provide alternative berthing plans. Thus, the development of the technique presented in this paper will provide the terminal operators with different robust berthing plans which are able to absorb limited incidences.

Berthing Allocation and Quay Crane Assignment: BAP + QCAP
Let be a set of incoming vessels; BAP + QCAP consists in obtaining an optimal (or near-optimal) schedule of the vessels by assigning mooring times, berthing positions, and QCs to each vessel. Our BAP + QCAP model is classified, according to the classification given by Bierwirth and Meisel [9] as follows.
(i) Spatial Attribute: Continuous Layout. We assume that the quay is a continuous line, so there is no partitioning of the quay and the vessel can berth at arbitrary positions within the boundaries of the quay. It must be taken into account that, for a continuous layout, berth planning is more complicated than for a discrete layout, but it better utilizes the quay space [9].
(ii) Temporal Attribute: Dynamic Arrival. Fixed arrival times are given for the vessels, so that vessels cannot berth before their expected arrival times.  (iv) Performance Measure: Wait and Handling Times. The objective is to minimize the sum of the waiting and handling times of all vessels .

Mathematical Problems in Engineering
Following, we introduce the notation used for each vessel ∈ ( Figure 1). The integer data variables are (i) : total number of QCs in the container terminal.
We assume all QCs carry out the same number of movements per time unit (movsQC), given by the container terminal,   (v) and : indexes of the first and last QC assigned to vessel , respectively.
In this study, the following assumptions are considered.
(i) All the information related to the waiting vessels is known in advance (arrival, priority, moves, and length).
(ii) Every vessel has a draft that is lower than or equal to the draft of the quay.
(iii) Movements of QCs along the quay as well as berthing and departure times of vessels are not considered since it supposes a constant penalty time for all vessels.
(iv) Simultaneous berthing is allowed, subject to the length of the berth.
Usually in container terminals, the number of QCs could vary during execution at the quay. This issue has been studied in Rodriguez-Molins et al. [5]. However, in this paper and without loss of generality, we study the robustness of the schedules assuming that the number of QCs assigned to one vessel does not vary along the moored time. Once a QC starts a task in a vessel, it must complete it without any pause or shift (nonpreemptive tasks). Thus, all QCs assigned to the same vessel have the same working time on the vessel ( = ℎ , ∀ = 1, . . . , , = 1). The following constraints must be accomplished.
(i) Moored time of vessel must be at least the same that its arrival time ( ≥ ).
(ii) There is a safe distance between two moored ships. We assume that each vessel has a 2.5% of this length at each side as a safe distance ( ) (Figure 1). This safe distance is added to the length of each vessel : := ℓ + 2 .
(iii) There must be enough contiguous space at berth to moor a vessel of length ( ).
(iv) There must be at least one QC assigned to each vessel. Furthermore, there is a maximum number of QCs that can be assigned to vessel (QC + ). This value, (QC + ), depends on the length of each vessel (ℓ ), since a safe distance is required between two contiguous 4 Mathematical Problems in Engineering QCs (safeQC) and the maximum number of QCs that the container terminal allows per vessel (maxQC) (equtaion (2)). Both safeQC and maxQC parameters are given by the container terminal: Our objective is to allocate all vessels according to several constraints minimizing the total waiting ( ) and handling or processing time ( ℎ ), known as the service time ( ), for all vessels:

Robust BAP + QCAP Model
Uncertainty and nondeterminism of real-world environments may cause difficulties in the initial plans made by the decision makers. In container terminals, the initial obtained schedules for the BAP + QCAP problem might become invalid due to different reasons: breakdowns in QCs, late arrivals of the vessels, extreme weather events, a lower ratio of movements per QC than expected, and so forth. The robustness concept means that, given a schedule, this initial schedule remains feasible when minor incidences occur in its actual scenario.
The usual disruptions to be considered in BAP + QCAP are the following: (i) early or late arrival of a vessel from its expected arrival time ( ); (ii) the handling time of a vessel is larger than its expected handling time (ℎ ).
In this paper, we focus just on the disruptions affecting the handling time which eventually delay the departure time. In case of incidences related to late arrivals, they could also be modeled as delays in the handling time of the vessels which eventually also delay their departure time. Definition 1. Given the possible disruptions, we consider that a schedule is robust if a disruption in one vessel does not affect or alter the mooring times of the other vessels.
The robustness of a schedule of BAP + QCAP might be guaranteed through two periods of time related to each vessel: waiting time of a vessel ( ) and buffer time after the departure of each vessel ( ) [28]. Without loss of generality, early arrivals are not taken into account since they only increase waiting times but they do not alter mooring times.
The schedule could absorb delays or breakdowns that do not exceed the sum of those two periods ( + ). Therefore, both times should be maximized in order to achieve the maximum robustness and ensure that there is no need to reschedule the involved vessels. However, it should be kept in mind that the first objective of the BAP + QCAP is to minimize the total service time of the incoming vessels ( + ℎ ). Therefore, following the proposal given by Davenport et al. [29], we focus on maximizing only the second period of time, buffer times (∑ ), to obtain robust schedules. Let be the vessels that succeed vessel and occupy some berth space of vessel ( ) or use any of QCs assigned to vessel ( ). The buffer time of vessel ( ) is the minimum difference ( ) between the departure time of vessel ( ) and the mooring time of vessel ( , ∈ ). In case there is no vessel in , the maximum buffer time is assigned to (an infinite value). Figure 2 shows an example of the buffer times ( ) assigned for each scheduled vessel as an empty rectangle: In this paper, we assume that the more handling time is, the more likely it is to suffer incidences. Therefore, in general, the larger the buffers are, the more robust the schedules are. Nevertheless, regarding the concept of decreasing productivity (or diminishing returns) presented in [22], there is a certain buffer size beyond which no more robustness is added to the schedule. Thereby, there is no need to assign large buffer times to each vessel. For instance, in Figure 2 suffers some breakdown or delay and so it becomes invalid this schedule. Furthermore, we consider that the magnitude of the incidence is related to the handling time of the vessel. Thus, the robustness measure of each vessel ( ∈ [0, 1]) is related to the buffer time and the average handling time ℎ * (equation (7)). It should be mentioned that other functions, for example, exponential function, could be adopted to define the robustness of each vessel: Given the robustness of each vessel, the robustness of a schedule ∈ [0, | |] is defined by (9), where is a weighting factor ( ≥ 1) which depends on historical data, if available. A = 1 value represents that vessel used to finish its tasks as expected, and > 1 value denotes that vessel used to be delayed: In this paper, we address the BAP + QCAP problem without knowledge of the incidences; thus, the weighting factor is the same for all the vessels ( = 1, ∀ ∈ ).
Example 2. Figure 3 shows two different schedules given the same set of 9 incoming vessels. Each vessel is labeled with its vessel's ID and the assigned QC number in brackets. Furthermore, the buffer time between vessels is also showed.
On the one hand, Figure 3(a) represents a robust schedule since limited incidences over any vessel could be absorbed. On the other hand, Figure 3(b) shows a schedule with the optimal solution according to the objective function . The latter schedule will be highly likely unfeasible if any incidence occurs.

Figures 3(a) and 3(b)
are an example of the well-known trade-off between optimality and robustness. However, a robust schedule is not only achieved by extending an optimal schedule over the time. A robust schedule must also consider an optimized allocation of vessels to achieve the maximum sum of buffer sizes with a proper distribution among all vessels. Note that the optimality is not directly the makespan of the schedule but the total service time (waiting and handling times).
An important issue in this paper is that there is no available information about how likely the incidences or breakdowns occur. Therefore, it is interesting that these buffers are proportionally distributed among all the vessels. Thereby, a third objective is introduced into the model in order to improve the robustness of a schedule: minimizing the standard deviation ( ) of the robustness measures of all vessels ( ∀ ∈ ): where is the average of the buffers of the schedule and | | is the number of incoming vessels. Both measures presented above, robustness of a schedule ( ) and standard deviation of these values ( ), represent 6 Mathematical Problems in Engineering the actual robustness of a schedule R to be maximized (see (11)). This measure guarantees the absorption of incidences that imply at most a delay of a R% of the weighted average handling time ( ℎ * ): Example 3. Figure 4 shows two different schedules with a similar robustness value ( = 0.7) but different standard deviations, 1 = 0.17 and 2 = 0.45. With these values, the first schedule has an actual robustness value of R 1 = 0.7 − 0.17 = 0.53. Thus, in average, this schedule guarantees that it could absorb incidences that imply at most a delay of the 53% of the average handling time of the vessels. In contrast, the second schedule has an actual robustness value of R 2 = 0.7 − 0.45 = 0.25; thus, it is able to absorb only incidences that imply at most a 25% of the average handling time of the vessels.
Surico et al. [30] presented a close function to measure the robustness of a schedule (avg( ) − ( ), ∀ ∈ ). avg( ) and ( ) denote the average and the standard deviation of the waiting times, respectively; is a constant weighting factor that must be set. However, this measure does not reflect the relationship between the handling or processing time of the task and the buffer times. Thus, to our best knowledge, there is no other study which, considering the BAP + QCAP with a continuous quay and dynamic arrivals, tackles the robustness without any previous knowledge about the incidences. Figure 4 shows two different schedules of 10 vessels with the same high value of the robustness measure. However, schedule of Figure 4(a) has a greater value for the standard deviation ( ) than schedule of Figure 4(b). Thereby, it is important to note that buffer times from schedule in Figure 4(a) are not equally distributed and this schedule will fail if an incidence which delays the departure time just 1 time unit over vessels 4 or 6 occurs or more than 3 time units over vessel 1. However, in the schedule in Figure 4(b), it is highly unlikely to be invalid since all vessels have enough buffer time after its schedule departure time.

Multiobjective Approach for the Robust BAP + QCAP
Three different objectives must be optimized to solve the robust BAP + QCAP: the service time ( ) (equation (5)), the robustness ( ) (equation (9)), and the standard deviation of the robustness measures ( ) (equation (10)). These objective functions must be normalized in order to apply the search process correctly. Equation (14) shows how to normalize the service time objective into the interval [0, 1] (̂) and it implies to normalize both the waiting timê(equation (12)) and the handling timê(equation (13)). On the one hand, the handling time is just a linear normalization since the maximum (ℎ + ) and minimum (ℎ − ) times are known by assigning the minimum (1) and the maximum number of QCs to vessel (QC + ). On the other hand, normalizing the waiting time requires to determine a maximum total waiting time ( ). In this case, value is the total waiting time of the incoming vessels when a first-come, first-served (FCFS) policy is applied, assigning 2 QCs to each vessel, and just one vessel is allowed in the berth at the same time (see, for example, Figure 5). The maximum total waiting time ( ) could also be obtained by assigning just one QC to each incoming vessel, but in that case, value would be too large and all the normalized waiting times would be close to zero: (2) Robustness objective function must also be normalized into the interval [0, 1] (̂) as defined by (15). The third objective, standard deviation of robustness measures, is already normalized due to the fact that values are already in the interval [0, 1] (equation (10)): Thereby, the objective function for the robust BAP + QCAP is to minimize the function (equation (16)). Each coefficient (0 ≤ ≤ 1) assigns different weights to each component or objective function in order to establish an aggregate function: These coefficients are subject to ∑ = 1. In a multiobjective optimization problem, usually there is no single solution wherein all its objectives are simultaneously optimized. However, there may exist a set of Pareto optimal solutions with different trade-offs between their objective functions. Pareto efficiency, or Pareto optimality, is a solution in which it is impossible to make any one criterion better off without making at least one criterion worse off [31]. Pareto optimal solutions are defined by means of the dominance concept. Considering the robust BAP + QCAP, let and be two different solutions; dominates if at least one of the following conditions is satisfied: Given a set of feasible solutions , a solution ∈ is Pareto optimal solution if it is nondominated by any other solution ∈ . The Pareto optimal set is the set of all the solutions that are Pareto optimal solutions [31].
In general, generating the Pareto optimal set is expensive computationally and it is often impracticable. Therefore, algorithms try to find a good approximation of the Pareto optimal set. In this work, we refer that each approximation as Pareto front, which contains solutions that, although are nondominated among them, could be dominated by other solutions not found by our algorithms.

Mathematical Formulation
A mixed integer programming (MIP) model is presented to solve the robust BAP + QCAP. The objective function of this model is to minimize (16). This mathematical model is based on the model presented in [10,28].
In the proposed model, denotes a sufficiently large number (as it is used in MIP). Furthermore, there are four auxiliary binary variables.
is a decision variable that indicates if vessel is located to the left of vessel on the berth ( = 1); = 1 indicates that vessel is moored before vessel in time. The auxiliary variable indicates whether the QC works (1) or not (0) on vessel ; = 1 denotes that the number of QCs assigned to vessel is : In the proposed model, there are four auxiliary binary variables.
is a decision variable that indicates if vessel is located to the left of vessel on the berth ( = 1); = 1 indicates that vessel is moored before vessel in time. The auxiliary variable indicates whether the QC works (1) or not (0) on vessel ; = 1 denotes that the number of QCs assigned to vessel is .
The constraints of this mathematical model are detailed below. Constraint (19) ensures that vessels must moor afer they arrive at the terminal: Constraints (20) and (21) establish the waiting and departure times according to and ℎ : Constraint (22) guarantees that a moored vessel does not exceed the length quay: The number of QCs to the vessel are assigned by means of constraints (23)-(28) as follows: Mathematical Problems in Engineering Constraints (29)-(31) establish the minimum handling time needed to load and unload their containers according to the number of assigned QCs: Constraint (32) ensures that QCs that are not assigned to vessel have = 0: Constraint (33) forces all assigned QCs to vessel working the same number of hours: Constraint (34) avoids that one QC is assigned to two different vessels at the same time: Constraints (35) and (36) force the QCs to be contiguously assigned (from up to ): The safety distance between vessels is taken into account by constraint (37) as follows: Constraint (38) avoids that one vessel uses a QC which should cross through the others QCs: Constraint (39) avoids that vessel moors while the previous vessel is still at the quay: Constraint (40) establishes the relationship between each pair of vessels avoiding overlaps: Constraint (41) ensures that the total waiting time of the schedule does not exceed the maximum total waiting time ( ): Constraints (42)-(44) assign the time between the departure time of vessel and the mooring time of vessel . For those vessels so that ̸ = 1, they are assigned as a value representing an unbounded time for the robustness: Constraints (45) and (46) set the value of the available buffer time after vessel and its robustness value, respectively: The decision variable (see constraint (47)) indicates if a vessel moors later than and, at the same time, the vessel intersects with the berth length occupied by vessel ( ):

Multiobjective Genetic Algorithms: MOGA + SA
Commonly approximations of the Pareto optimal sets of a multiobjective optimization problem are obtained by means of multiobjective evolutionary algorithms [31]. Furthermore, nowadays, metaheuristics are usually hybridized with other techniques or algorithms in order to enhance their effectiveness and performance [23,24]. One of the most common forms of hybrid genetic algorithm involves incorporating local search to a canonical genetic algorithm. Genetic algorithm is used to perform global exploration among the population, and local search is used to perform local exploitation around the chromosomes. Because of the complementary properties of genetic algorithms and local search methods, the hybrid approach often outperforms either methods operating alone [32]. Thereby, a hybrid multiobjective genetic algorithm (MOGA) has been implemented in this paper. The NSGAII schema has been extended with a multiobjective local search based on the multiobjective simulated annealing proposed by Bandyopadhyay et al. [33] (AMOSA), hereinafter named as MOGA + SA (see Algorithm 1). Moreover, two different schemes from the literature have been assessed NSGAII [26] and SPEA2+ [27].
The same chromosome structure is used in these three MOGAs. This chromosome has as many genes as incoming vessels (| |). Each gene consists of three values (see Figure 6): (1) the ID of the next vessel to dispatch ( ); (2) the number of QCs assigned ( ); (3) the buffer size after this vessel ( ).
It should be noted that each gene must be composed of feasible values with respect to vessel . That is, according to the problem constraints, each vessel can be assigned at most QC + cranes. Therefore, 1 ≤ ≤ QC + . Likewise, if the berth length is , then ≤ ≤ − − .
In the following subsections, genetic operators that are used by the implementations of NSGAII and SPEA2+ are described.

Decoding and Evaluation of One Chromosome/Solution.
The structure of the chromosome, specifically the order of the vessels, is used as a dispatching rule. Hence, we use the following decoding algorithm: the genes are visited from left to right in the chromosome sequence. For each gene ( , , and ), the vessel is scheduled at the earliest mooring time with consecutive QCs available, so that none of the constraints are violated. In case there are several positions available at the earliest mooring time, the one closest to the berth extremes is selected. After the departure of the vessel ( ), it is ensured that there are time units where no other vessel (∀ ∈ , ̸ = ) uses the QCs assigned to vessel nor moors where vessel does [ , + ).
Once a valid mooring time ( ) and initial position ( ) have been assigned to each vessel , the fitness of the chromosome (equation (16)) is obtained by computing each one of the objective functions: total service time (̂), robustness (̂), and standard deviation of the robustness ( ).

Generation of Initial Population. Construction of initial population (
procedure) is performed so that the service time of a percentage of the initial population (GA parameter) is at least as good as the solution provided by the FCFS policy. The other chromosomes (or solutions) are constructed by instantiating each gene in the following way.
(i) Vessel identifier ( ): an integer, between 1 and , is randomly chosen. Two genes of the same chromosome cannot have the same vessel identifier.
(iii) Buffer size ( ): the initial buffer size is 0 for all genes of the initial population.
Once all chromosomes in the initial population have been instantiated, their fitness values are obtained as described in Section 6.1. Furthermore, the Pareto front X is updated considering all these chromosomes. Let be a chromosome (or solution); is added to the Pareto front X if there is no other solution ∈ X such that dominates . If is added to X, then all solutions dominated by are removed from X.

Evolution of One Population.
In each iteration of the MOGA, a new population is built from the previous one (or the initial) by applying the genetic operators of selection, reproduction, and replacement. The proposed approach follows the scheme: (i) selection: all chromosomes in the actual population are randomly grouped into pairs; (ii) reproduction: (1) each one of these pairs is mated or not according to the crossover probability generating two offspring; (2) each offspring, or parent if the parents were not mated, undergoes mutation in accordance with the mutation probability ; (iii) replacement: after evaluating the chromosomes previously generated, a tournament selection (4 : 2) is carried out among each pair of parents and their offspring as a replacement.

Crossover.
The crossover operator receives one pair of chromosomes ( 1 and 2 ), which are in the current population and have been randomly selected. The objective of this operator is to construct two offspring chromosomes ( 1 and 2 ). For that, each time the crossover operation is performed, the following steps are made.
(2) Each gene in chromosomes 1 and 2 which is in position , 1 ≤ < 2 , is copied to the same position in chromosomes 1 and 2 , respectively.
(3) Each gene in chromosomes 1 and 2 which is in position , 1 ≤ < 1 , is copied to the same position in chromosomes 1 and 2 , respectively.
(4) Each gene in chromosomes 1 and 2 which is in position , 2 ≤ ≤ , is copied to the same position in chromosomes 1 and 2 , respectively. Figure 7 is a graphical representation of the procedure that is used to perform the crossover operation, which is based on the technique generalized position crossover [34] that is commonly used in permutation based encodings.
In one chromosome there cannot be two genes with the same vessel identifier. Therefore, if the vessel identifier in the gene that will be copied already exists in the offspring ( 1 / 2 )  Once the vessel identifier of the selected gene does not exist in the offspring, then the gene is copied to the offspring in the corresponding position.

Mutation.
Mutation operation is performed on one chromosome, following these steps.
(2) Genes that are in positions between 1 and 2 (both included) are shuffled.
(3) The number of QCs in each gene located between 1 and 2 (both included) is modified by a feasible random value with respect to the vessel in the same gene.
(4) The buffer size in each gene located between 1 and 2 , both included, is modified by a random value that  Figure 7: Crossover operation. is between 0 and the average handling time ℎ * , of the vessel in the same gene. Figure 8 shows how the offspring , which has been obtained after the crossover operation, is mutated. First, two values 1 = 2 and 2 = 4 are selected randomly. Then, all genes between both positions are shuffled. Gene 2 is moved to position 4, gene 3 to position 1, and gene 4 to position 3. Finally, the number of QCs and the buffer size of each gene in position , 2 ≤ ≤ 4, are modified by selecting feasible random values for each one.
6.6. Local Search. The multiobjective simulated annealing presented by Bandyopadhyay et al. [33] has been modified and included into the MOGA as a local search to solve the robust BAP + QCAP. The neighborhood structure of a solution takes advantage of the chromosome structure and their neighbors are obtained by changing the values of the variables presented in its genes ( , , and ). Algorithm 2 describes how to modify a given solution in order to create a neighbor . In this process, two operations are applied to solution : (i) interchanging the position of two vessels (V 1 and V 2 ) in the chromosome, randomly chosen, (ii) changing the values of number of QCs assigned ( ) and buffer size ( ) of a vessel randomly chosen. This multiobjective simulated annealing algorithm (mosa function) is computed every iterations. It receives, as parameter, the Pareto front X of the actual iteration of the MOGA + SA. As a result, it returns two different sets of solutions or schedules: (i) a Pareto front X where the solutions from X have been improved to obtain a local optimal following the AMOSA scheme, (ii) a new set consisting of the new nondominated solutions found in the search which are part of X .
Unlike AMOSA [33], the simulated annealing algorithm implemented in this paper makes use of a different clustering method (see Algorithm 3) based on the crowding distance used in the NSGAII algorithm. This clustering method selects the representative solutions of the population according to the density of solutions surrounding a particular solution. After this local search process is performed, the solutions in set replace the solutions in population whose crowding distances are the lowest ones. The solution to be replaced is chosen by means of the same clustering method used in the simulated annealing. The purpose is to improve the quality of the population by keeping the solutions that are most spread around the search space.   Minimum temperature ( min ) 0.001 Annealing factor ( ) 0.9

Termination criterion < min
Cycle length 2

Evaluation
As no benchmark is available in the literature, the experiments were performed in a corpus of 100 instances randomly generated, where parameters (maxQC, safeQC, etc.) follow the suggestions of container terminal operators. All these benchmarks are freely available at http://gps.webs.upv.es/ bap-qcap/. Each one is composed of a queue from 100 vessels. These instances follow an exponential distribution for the interarrival times of the vessels (scale parameter = 20). The number of required movements and length of vessels are uniformly generated in [100, 1000] and [70, 400], respectively. In all cases, the berth length ( ) was fixed to 700 meters; the number of QCs was 7 (corresponding to a determined MSC berth line) and the maximum number of QCs per vessel was 5 (maxQC); the safe distance between QCs (safeQC) was 35 meters and the number of movements that QCs carry out was 2.5 (movsQC) per time unit. The approaches developed in this paper, NSGAII, SPEA2+, and MOGA + SA, were coded using C++; their settings are showed in Tables 1(a) and 1(b). Due to the stochastic nature of the GA process, each instance was solved 30 times and the results show the average obtained values. The mathematical model was coded and solved by using IBM ILOG CPLEX Optimization Studio 12.5. Due to the fact that the square root function defines concave region, standard deviation function could not be introduced into the objective function in the mathematical solver. They were solved on an Intel i7-2600 3.4 Ghz with 8 Gb RAM.
CPLEX is able to obtain a schedule of an instance for a given value. Algorithm 4 describes how to obtain a Pareto front using CPLEX solver for a given instance in order to be compared with the MOGA. Figure 9 shows the Pareto fronts obtained of a representative instance by both the MOGA + SA and CPLEX. In this experiment, the timeout for the CPLEX solver was set to 1000 seconds for each value. It is important to note that the greater the incoming vessels are, the fewer the solutions obtained by CPLEX solver are. Given this timeout, CPLEX was only able to get optimal solutions when = 0.0 and the incoming vessels were set to 10 and 20. Considering the Pareto fronts obtained by MOGA + SA and CPLEX, they were very similar with a queue of 10 vessels (see Figure 9(a)). However, for instance, with a queue 20 vessels, the solutions obtained by CPLEX were not able to reach the quality of the Pareto front of MOGA + SA (see Figure 9(b)). Furthermore, it turned out that for 40 incoming vessels just one nonoptimal solution was obtained (see Figure 9(d)) and even more there was no solution with 50 vessels.
Multiobjective optimization algorithms are not comparable directly since there is no a unique optimal solution. Zitzler et al. [35] propose different measures to compare Pareto front approximations. Among these measures, the size of the dominated space or the hypervolume is one of the most used measures to differentiate two algorithms [36]. This measure is related to a reference point and it is set according to the suggestion of While et al. [36]. To this end, for each objective, the worst value from any of the sets being compared is chosen and increased by an value.
Comparison among these different schemes has been performed by using the Kruskal-Wallis' nonparametric statistical test, according to Zitzler et al. [35]. This test assesses whether there are significant differences among different sets of values: in this case, sets of hypervolume measures. Table 2 shows the values obtained for this test given the results after solving five instances of 50 vessels with the three different algorithms. Kruskal-Wallis test revealed a significance effect of the algorithms on the hypervolumes ( value < 0.01).  As there was a significance difference among them, a post-hoc test using a pairwise comparison test (Wilcoxon) with Bonferroni correction was carried out and showed the significant differences between the different algorithms. As an example, Table 3 shows the results of the Wilcoxon test for the fifth instance. Note that, MOGA + SA algorithm produces Pareto fronts which are statistically different with respect to the other algorithms. Examining the average values in Table 2, it can be determined that MOGA + SA obtained better Pareto front approximations. Figure 10 shows the Pareto fronts obtained by NSGAII and MOGA + SA algorithms. The schedules with the minimum and maximum values for each objective are highlighted by circles. It is important to note that MOGA + SA algorithm   The performance of the schedules obtained by our approach (MOGA + SA) was evaluated by generating actual scenarios with some incidences in the actual handling time of the vessels. An incidence over a vessel is modeled as a delay in the handling time of vessel . This incidence is absorbed if there is enough buffer time behind vessel as to not alter the mooring time of the subsequent vessels. For each instance, the vessels that vary their handling times were uniformly chosen among all the scheduled vessels.

Mathematical Problems in Engineering
In this experiment, 100 instances of 100 vessels were evaluated. For each instance, three different schedules were chosen from the Pareto front according to their robustness (see Table 4(a)): the one with the minimum robustness ( min ), the one with the maximum robustness ( Max ), and one intermediate robust schedule ( where min < < Max ). Likewise, three schedules were chosen according to their service time (see Table 4(b)) and other three schedules according to their standard deviation of the robustness (see Table 4(c)).
The incidences (or delays, ) introduced were randomly chosen from different ranges. These ranges vary from a minimum value (1) to a maximum value, which is related to the handling time (ℎ ) of the vessel affected by the incidence (see first column in Table 4). For each range, 100 incidences were uniformly created and applied to the three schedules of each instance.  Table 4(b). The lower , the lower incidences absorbed due to the fact that either there would not be buffers among vessels or there would be small buffers.
Table 4(c) shows the percentage of incidences absorbed choosing three schedules by their standard deviation values. As expected, the highest percentages of incidences absorbed were obtained with the lowest values of standard deviation, for example, 99.34% with delays in the range [1, 0.8ℎ ]. In general, schedules with the lowest standard deviation are related to those schedules with the greatest buffers proportionally distributed among all vessels (the most robust schedules).
The percentage of incidences absorbed by the most robust schedules using or not the local search algorithm are showed in Table 5. In this experiment, a timeout of 30 seconds was set for both algorithms. It is important to note that adding the local search to the multiobjective genetic algorithm allowed to increase the incidences absorbed in all the ranges. For instance, in range [1, 1.0ℎ ], NSGAII was able to absorb 95.33% of incidences, whereas the MOGA + SA was able to absorb 97.34% of incidences.

Conclusions
The competitiveness among container terminals causes the need to improve the efficiency of each one of the subprocesses or scheduling problems that are performed within them. However, this efficiency is affected by the uncertainty of the environment. This uncertainty might provoke delays in the arrivals of the vessels or handling times greater than expected due to extreme weather events, breakdowns in engines, delays, and so forth. Furthermore, these scheduling problems are even harder since they are interrelated and sometimes there is no previous knowledge about these incidences. To this end, we introduce the robustness into these scheduling problems. In this paper, we introduce the robustness into one of the main scheduling problems in container terminals: berth allocation and quay crane assignment problems. Its objective function is to minimize the total service time of the incoming vessels. The robustness, as second objective function, has been modeled as a measure related to the likelihood of a schedule to absorb incidences. This robustness has been related to the operational buffers found after each assigned vessel. The greater the operational buffers are, the higher the robustness of the schedule is. However, due to the lack of the knowledge about incidences, operational buffers should be distributed among vessels proportionally,  and thus the third objective managed in this way is to minimize the standard deviation of the robustness measurements. In this paper, a mixed integer programming (MIP) model and a new hybrid multiobjective genetic algorithm (MOGA + SA) were developed for the dynamic and continuous robust BAP + QCAP. They were compared with two well-known multiobjective genetic algorithms (MOGAs): NSGAII and SPEA2+. In multiobjective optimization problems there is no a unique optimal solution, and it is necessary to assess the trade-off between all the objectives by using the Pareto front. Visualizing Pareto fronts provides container terminals operators with a helpful system to decide which schedule is better depending on the actual state of the container terminal.
The results showed that the MIP model was able to obtain robust and efficient schedules up to 10 incoming vessels. However, MOGA + SA achieved better Pareto fronts than the MIP model for queues of incoming vessels greater than or equal to 20 vessels. Thereby, the schedules obtained by MOGA + SA were more efficient and robust than the schedules obtained by the MIP model. Furthermore, the MIP model was unable to find any solution with a given timeout for a queue of 50 incoming vessels. Additionally, differences between the MOGAs have been assessed by means of nonparametric statistical tests. It turned out to be that MOGA + SA obtained better Pareto fronts according to the hypervolume measures. Furthermore, different sets of incidences were simulated into the schedules obtained by the NSGAII and the MOGA + SA. The results returned that the schedules obtained by MOGA+SA were more robust due to the fact that they could absorb more incidences.