# Robust Scheduling for Berth Allocation and Quay Crane Assignment Problem

Authors: M. Rodriguez-Molins; M. A. Salido; F. Barber
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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+.


---

## Body

## 1. Introduction

Within a container terminal, operations related to move containers can be divided into four different subsystems (ship-to-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–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–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 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 (Section3) 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.


## 2. Berthing Allocation and Quay Crane Assignment: BAP + QCAP

LetV be a set of incoming vessels; BAP + QCAP consists in obtaining an optimal (or near-optimal) schedule of the vessels V 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.(iii)
Handling Time Attribute: Unknown in Advance. The handling time of a vessel depends on the number of assigned QCs (QCAP) and the moves required.(iv)
Performance Measure: Wait and Handling Times. The objective is to minimize the sum of the waiting and handling times of all vessels V.Following, we introduce the notation used for each vesseli∈V (Figure 1). The integer data variables are(i)
K: 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,(ii)
L: total length of the berth in the container terminal,(iii)
ai: arrival time of the vessel i at port,(iv)
ci: number of required movements to load and unload containers of vessel i,(v)
li: vessel length.Figure 1
Data related to one vessel.The decision variables are(i)
mi: mooring time of i. Thus, waiting time (wi) of vessel i is calculated as wi=mi-ai,(ii)
pi: berthing position where vessel i moors,(iii)
qi: number of assigned QCs to vessel i,(iv)
uik: indicates whether the QC k (1≤k≤K) works (1) or not (0) on the vessel i,(v)
nik: denotes that the number of QCs assigned to vessel i is k QCs (nik=1), For instance, if vessel 3 has been assigned 4 QCs, then n34=1 and the others QCs n3k=0, ∀k=1,…,K, k≠4.The variables derived from the previous ones are(i)
Hik: loading and unloading time at quay (handling time) of vessel i using k QCs (1≤k≤K). This handling time depends on ci and it is defined by
(1)Hik=cik*movsQC∀i∈V,∀k=1,…,K,(ii)
hi: required handling time of vessel i when qi QCs are assigned to it. This value is set by means of Hiqi,(iii)
tik: working time of the kth QC (1≤k≤K) that is assigned to vessel i,(iv)
di: departure time of vessel idi=mi+hi,(v)
si and ei: indexes of the first and last QC assigned to vessel i, 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 i have the same working time on the vessel tik=hi,∀k=1,…,K,uik=1.The following constraints must be accomplished.(i)
Moored time of vesseli must be at least the same that its arrival time mi≥ai.(ii)
There is a safe distance between two moored ships. We assume that each vesseli has a 2.5% of this length at each side as a safe distance ηi (Figure 1). This safe distance is added to the length of each vessel i: li∶=li+2ηi.(iii)
There must be enough contiguous space at berth to moor a vesseli of length li.(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 vesseli (QCi+). This value, QCi+, depends on the length of each vessel (li), since a safe distance is required between two contiguous 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:
(2)QCi+=min⁡maxQC,max⁡1,lisafeQC∀i∈V.Our objective is to allocate all vessels according to several constraints minimizing the total waiting (Tw) and handling or processing time (Th), known as the service time (Ts), for all vessels:
(3)Tw=∑i∈V‍wi,(4)Th=∑i∈V‍hi,(5)Ts=Tw+Ts.


## 3. 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 vesseli from its expected arrival time (ai);(ii)
the handling time of a vesseli is larger than its expected handling time (hi).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 (wi) and buffer time after the departure of each vessel (bi) [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(wi+bi). 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 (wi+hi). Therefore, following the proposal given by Davenport et al. [29], we focus on maximizing only the second period of time, buffer times (∑‍bi), to obtain robust schedules.Letφi be the vessels that succeed vessel i and occupy some berth space of vessel i (φip) or use any of QCs assigned to vessel i (φiq). The buffer time of vessel i (bi) is the minimum difference (τij) between the departure time of vessel i (di) and the mooring time of vessel j (mj,j∈φi). In case there is no vessel in φi, the maximum buffer time is assigned to bi (an infinite value). Figure 2 shows an example of the buffer times bi assigned for each scheduled vessel as an empty rectangle:
(6)φip=∀j∈V,mj≥di∧pi,pi+li∩pj,pj+lj≠⌀∀i∈Vφiq=∀j∈V,mj≥di∧∃k,1≤k≤K∧uik=1∧ujk=1∀i∈Vφi=φip∪φiq∀i∈Vτij=mj-di∀i∈V,∀j∈φibi=+∞,|φi|=0min⁡j∈φiτij,otherwise∀i∈V.Figure 2
Buffer timesbi given an example schedule.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, vessel 1 would not need 8 time units of buffer time (b1) since its handling time is only 3 time units. It is not likely that this vessel would suffer a delay of that magnitude. However, vessel 2, with a handling time of 8 time units, has only 2 time units of buffer time (b2). In this case, it is highly likely that this vessel 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 vesseli (ri∈[0,1]) is related to the buffer time bi and the average handling time hi* (equation (7)). It should be mentioned that other functions, for example, exponential function, could be adopted to define the robustness of each vessel:
(7)hi*=ci(1+QCi+)/2movsQC.Given the robustness of each vessel, the robustness of a scheduleR∈[0,|V|] is defined by (9), where ωi is a weighting factor (ωi≥1) which depends on historical data, if available. A ωi=1 value represents that vessel i used to finish its tasks as expected, and ωi>1 value denotes that vessel i used to be delayed:
(8)ri=min⁡1,biωihi*,∀i∈V,(9)R=∑i∈V‍ri.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 (ωi=1,∀i∈V).Example 2.
Figure3 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 Ts. The latter schedule will be highly likely unfeasible if any incidence occurs.Two possible schedules given the same incoming vessels.
(a)
Robust schedule
(b)
Optimal schedule according toTsFigures3(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 (ri∀i∈V):
(10)σ=1|V|∑i∈V‍(ri-r¯),
where r¯ is the average of the buffers of the schedule and |V| is the number of incoming vessels.Both measures presented above, robustness of a schedule (R) and standard deviation of these values (σ), represent 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 (ωihi*):
(11)R=R-σ.Example 3.
Figure4 shows two different schedules with a similar robustness value (R=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 R1=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 R2=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.Two different schedules with similar robustness and different standard deviation.
(a)
High standard deviation
(b)
Low standard deviationSurico et al. [30] presented a close function to measure the robustness of a schedule (avg(wi)-ασ(wi), ∀i∈V). avg(wi) and σ(wi) 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.Figure4 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.


## 4. Multiobjective Approach for the Robust BAP + QCAP

Three different objectives must be optimized to solve the robust BAP + QCAP: the service time(Ts) (equation (5)), the robustness (R) (equation (9)), and the standard deviation of the robustness measures σ(R) (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] (T^s) and it implies to normalize both the waiting time T^w (equation (12)) and the handling time T^s (equation (13)). On the one hand, the handling time is just a linear normalization since the maximum (hi+) and minimum (hi-) times are known by assigning the minimum (1) and the maximum number of QCs to vessel i(QCi+). On the other hand, normalizing the waiting time requires to determine a maximum total waiting time (WF). In this case, WF 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 (WF) could also be obtained by assigning just one QC to each incoming vessel, but in that case, WF value would be too large and all the normalized waiting times would be close to zero:
(12)T^w=1WF∑i∈V‍mi-aiT^w∈0,1,(13)T^h=1|V|∑i∈V‍hi-hi-hi+-hi-T^h∈0,1,(14)T^s=T^w+T^h2T^s∈0,1.Figure 5
Schedule generated to obtain the maximum value of waiting timeWF.Robustness objective function must also be normalized into the interval[0,1] (R^) as defined by (15). The third objective, standard deviation of robustness measures, is already normalized due to the fact that ri values are already in the interval [0,1] (equation (10)):
(15)R^=R|V|R^∈[0,1].Thereby, the objective function for the robust BAP + QCAP is to minimize the functionF (equation (16)). Each coefficient λi0≤λi≤1 assigns different weights to each component or objective function in order to establish an aggregate function:
(16)F=λ1T^s-λ2R^+λ3σ.These coefficientsλi are subject to ∑i‍λi=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 x and y be two different solutions; x dominates y if at least one of the following conditions is satisfied:
(17)T^sx<T^sy∧R^x≥R^y∧σx≤σy,T^sx≤T^sy∧R^x>R^y∧σx≤σy,T^sx≤T^sy∧R^x≥R^y∧σx<σy.Given a set of feasible solutionsD, a solution x∈D is Pareto optimal solution if it is nondominated by any other solution x′∈D. 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.


## 5. 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,M denotes a sufficiently large number (as it is used in MIP). Furthermore, there are four auxiliary binary variables. zijx is a decision variable that indicates if vessel i is located to the left of vessel j on the berth (zijx=1); zijy=1 indicates that vessel i is moored before vessel j in time. The auxiliary variable uik indicates whether the QC k works (1) or not (0) on vessel i; nik=1 denotes that the number of QCs assigned to vessel i is k:
(18)∀i,j∈Vi≠j∀k=1,…,Kzijx,zijy,uik,nik0/1integer.In the proposed model, there are four auxiliary binary variables.zijx is a decision variable that indicates if vessel i is located to the left of vessel j on the berth (zijx=1); zijy=1 indicates that vessel i is moored before vessel j in time. The auxiliary variable uik indicates whether the QC k works (1) or not (0) on vessel i; nik=1 denotes that the number of QCs assigned to vessel i is k.The constraints of this mathematical model are detailed below. Constraint (19) ensures that vessels must moor afer they arrive at the terminal:
(19)∀i∈Vmi≥ai.
Constraints (20) and (21) establish the waiting and departure times according to mi and hi:
(20)∀i∈Vwi=mi-ai,(21)∀i∈Vdi=mi+hi.
Constraint (22) guarantees that a moored vessel does not exceed the length quay:
(22)∀i∈Vpi+li≤L.
The number of QCs to the vessel i are assigned by means of constraints (23)–(28) as follows:
(23)∀i∈Vqi=∑k=1K‍uik,(24)∀i∈V∑k=1K‍nik=1,(25)∀i∈V∑k=1K‍nikk=qi,(26)∀i∈V1≤qi≤QCi+,(27)∀i∈V1≤si≤ei≤K,(28)∀i∈Vqi=ei-si+1.
Constraints (29)–(31) establish the minimum handling time needed to load and unload their containers according to the number of assigned QCs:
(29)∀i∈V∑k=1K‍tikmovsQC≥ci(30)∀i∈V∑k=1K‍nikHik=hi(31)∀i∈Vhi=max⁡∀k=1,…,Ktik.
Constraint (32) ensures that QCs that are not assigned to vessel i have tik=0:
(32)∀i∈V∀k=1,…,Ktik-Muik≤0.
Constraint (33) forces all assigned QCs to vessel i working the same number of hours:
(33)∀i∈V∀k=1,…,Khi-M(1-uik)-tik≤0.
Constraint (34) avoids that one QC is assigned to two different vessels at the same time:
(34)∀i,j∈V∀k=1,…,Kuik+ujk+zijx≤2.
Constraints (35) and (36) force the QCs to be contiguously assigned (from si up to ei):
(35)∀i∈V∀k=1,…,KM(1-uik)+(ei-k)≥0,(36)∀i∈V∀k=1,…,KM(1-uik)+(k-si)≥0.
The safety distance between vessels is taken into account by constraint (37) as follows:
(37)∀i,j∈Vi≠jpi+li≤pj+M1-zijx.
Constraint (38) avoids that one vessel uses a QC which should cross through the others QCs:
(38)∀i,j∈Vi≠jei+1≤sj+M1-zijx.
Constraint (39) avoids that vessel j moors while the previous vessel i is still at the quay:
(39)∀i,j∈Vi≠jdi≤mj+M1-zijy.
Constraint (40) establishes the relationship between each pair of vessels avoiding overlaps:
(40)∀i,j∈Vi≠jzijx+zjix+zijy+zjiy≥1.
Constraint (41) ensures that the total waiting time of the schedule does not exceed the maximum total waiting time (WF):
(41)∑i∈V‍wi≤WF.
Constraints (42)–(44) assign the time between the departure time of vessel i and the mooring time of vessel j. For those vessels j so that zijt≠1, they are assigned M as a value representing an unbounded time for the robustness:
(42)∀i,j∈Vi≠jzijt=zijx+zjix+zijy,(43)∀i,j∈Vi≠j∧(zijt=0∨zijt=2)τij=M,(44)∀i,j∈Vi≠j∧zijt=1di+τij=mj+M1-zijy.
Constraints (45) and (46) set the value of the available buffer time after vessel i and its robustness value, respectively:
(45)∀i∈Vbi=min⁡min⁡j∈Vi≠j⁡τij,hi*,(46)∀i∈Vrihi*=bi.
The decision variable zijt (see constraint (47)) indicates if a vessel j moors later than i and, at the same time, the vessel j intersects with the berth length occupied by vessel i (zijt):
(47)∀i,j∈Vi≠j0≤zijt≤2


## 6.  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].Algorithm 1:Implementation of MOGA + SA.
Input:  I: instance of robust BAP + QCAP;
popsize: number of chromosomes;
k: number generations for local search
Output:  X: set of nondominated schedules
{Generate the parent populationP0}
P0← generateInitialPopulation(popsize, I)
{Generate the offspring populationQ0}
Q0← evolvePopulation(P, popsize)
t←0
while No termination criterion is satisfied do
unionSet← makeUnionSet(Pt, Qt)
F← fastNondominatedSort(unionSet)
X← updateParetoFront(F0)
{Create the next parent populationPt+1}
i←0
Pt+1←∅
while  i<|F|∧|Pt+1|+|Fi|<popsize  do
{Add theith nondominated front (Fi) into the parent population Pt+1}
Pt+1←Pt+1∪Fi
i←i+1
end while
if  |Pt+1|<popsize  then
{SortFi according to the crowding distance measure}
crowdingDistance(Fi)
sort(Fi)
{Add the firstpopsize-|Pt+1| elements of Fi}
j←0
while  j<|Fi|∧|Pt+1|<popsize  do
Pt+1←Pt+1∪{Fi[j]}
j←j+1
end while
end if
{Use selection, crossover, and mutation to create a new populationQt+1}
Qt+1← evolvePopulation(Pt+1, popsize)
{Perform the local search eachk generations.}
if  t%k=0  then
X′,Sn← mosa(X)
{Assign the new Pareto front toX}
X←X′
Qt+1← clustering(Qt+1, popsize-|Sn|)
{The new solutions found by the local search are kept in the population}
Qt+1←Qt+1∪Sn
end if
{Increase the number of iterations}
t←t+1
end while
return  Schedule of each element of the Pareto front XThe same chromosome structure is used in these three MOGAs. This chromosome has as many genes as incoming vessels (|V|). Each gene consists of three values (see Figure 6): (1) the ID of the next vessel to dispatch (i); (2) the number of QCs assigned (qi); (3) the buffer size after this vessel (bi).Figure 6
Structure of one gene of a chromosome.It should be noted that each gene must be composed of feasible values with respect to vesseli. That is, according to the problem constraints, each vessel i can be assigned at most QCi+ cranes. Therefore, 1≤qi≤QCi+. Likewise, if the berth length is L, then ηi≤pi≤L-li-ηi.In the following subsections, genetic operators that are used by the implementations of NSGAII and SPEA2+ are described.


### 6.1. 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 (i, qi, and bi), the vessel i is scheduled at the earliest mooring time with qi 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 i (di), it is ensured that there are bi time units where no other vessel j (∀j∈V,j≠i) uses the QCs assigned to vessel i nor moors where vessel i does [pi,pi+li).Once a valid mooring time (mi) and initial position (pi) have been assigned to each vessel i, the fitness of the chromosome (equation (16)) is obtained by computing each one of the objective functions: total service time (T^s), robustness (R^), and standard deviation of the robustness (σ).


### 6.2.  Generation of Initial Population

Construction of initial population (generateInitialPopulation 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 (i): an integer, between 1 and N, is randomly chosen. Two genes of the same chromosome cannot have the same vessel identifier.(ii)
Number of QCs (qi): an integer, between 1 and QCi+, is randomly chosen.(iii)
Buffer size (bi): 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 Section6.1. Furthermore, the Pareto front X is updated considering all these chromosomes. Let x be a chromosome (or solution); x is added to the Pareto front X if there is no other solution y∈X such that y dominates x. If x is added to X, then all solutions dominated by x are removed from X.


### 6.3.  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 Pc generating two offspring; (2) each offspring, or parent if the parents were not mated, undergoes mutation in accordance with the mutation probability Pm;(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.


### 6.4.  Crossover

The crossover operator receives one pair of chromosomes (P1 and P2), which are in the current population pop and have been randomly selected. The objective of this operator is to construct two offspring chromosomes (O1 and O2). For that, each time the crossover operation is performed, the following steps are made.(1)
Two cross points are randomly chosen,k1 and k2(1≤k1<k2≤N).(2)
Each gene in chromosomesP1 and P2 which is in position p, k1≤p<k2, is copied to the same position in chromosomes O1 and O2, respectively.(3)
Each gene in chromosomesP1 and P2 which is in position p, 1≤p<k1, is copied to the same position in chromosomes O1 and O2, respectively.(4)
Each gene in chromosomesP1 and P2 which is in position p, k2≤p≤N, is copied to the same position in chromosomes O1 and O2, respectively.Figure7 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.Figure 7
Crossover operation.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 (O1/O2) during steps 2 or 3, a new gene must be selected from the chromosome parent (P1/P2).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.


### 6.5.  Mutation

Mutation operation is performed on one chromosome, following these steps.(1)
Two positions (k1 and k2) of the chromosome are randomly chosen (1≤k1<k2≤N).(2)
Genes that are in positions betweenk1 and k2 (both included) are shuffled.(3)
The number of QCs in each gene located betweenk1 and k2 (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 betweenk1 and k2, both included, is modified by a random value that is between 0 and the average handling time hi*, of the vessel i in the same gene.Figure8 shows how the offspring oi, which has been obtained after the crossover operation, is mutated. First, two values k1=2 and k2=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 p, 2≤p≤4, are modified by selecting feasible random values for each one.Figure 8
An example of mutation operation.


### 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 (i, qi, and bi). Algorithm 2 describes how to modify a given solution cur in order to create a neighbor new. In this process, two operations are applied to solution cur:(i)
interchanging the position of two vessels (v1 and v2) in the chromosome, randomly chosen,(ii)
changing the values of number of QCs assigned (qi) and buffer size (bi) of a vessel i randomly chosen.Algorithm 2:createNeighbor.
Input:  cur: Actual solution/chromosome;
Output:  new: Neighbor solution/chromosome;
Copy chromosomecur into new
{Interchange the order of two vessels}
Vesselsv1 and v2 randomly chosen from new
Interchange position ofv1 and v2 in new
{Change the buffer assigned after its departure}
Randomly choose other vesselv from new
r←randInteger(0,2)
if  r=0  then
bv=0
else if  r=1  then
bv=hi*
else
if  rand(0,1)≤0.5  then
{Decrease a 10% the buffer size}
bv←max⁡(0,bv-hi*10)
else
{Increase a 10% the buffer size}
bv←min⁡(hi*,bv+hi*10)
end if
end if
if  rand(0,1)≤0.5  then
{Change number of assigned QCs}
qv←rand(2,QCv+)
end if
return Neighbour chromosome/solution newThis multiobjective simulated annealing algorithm (mosa function) is computed every k 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 frontX′ where the solutions from X have been improved to obtain a local optimal following the AMOSA scheme,(ii)
a new setSn 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 Sn set replace the solutions in population tmpPop 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.Algorithm 3:Clustering.
Input:  P: population;
Nm; maximum number of chromosomes;
Output:  P′: population with exactly Nm chromosomes;
{calculate the crowding distance for each element/chromosome inP}
crowdingDistance(P)
{sort the elements in ascending order}
sort(P)
{choose the biggestNm elements according to the crowding distance}
P′←Nm elements of P with the biggest crowding distances
return  P′


## 6.1. 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 (i, qi, and bi), the vessel i is scheduled at the earliest mooring time with qi 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 i (di), it is ensured that there are bi time units where no other vessel j (∀j∈V,j≠i) uses the QCs assigned to vessel i nor moors where vessel i does [pi,pi+li).Once a valid mooring time (mi) and initial position (pi) have been assigned to each vessel i, the fitness of the chromosome (equation (16)) is obtained by computing each one of the objective functions: total service time (T^s), robustness (R^), and standard deviation of the robustness (σ).


## 6.2.  Generation of Initial Population

Construction of initial population (generateInitialPopulation 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 (i): an integer, between 1 and N, is randomly chosen. Two genes of the same chromosome cannot have the same vessel identifier.(ii)
Number of QCs (qi): an integer, between 1 and QCi+, is randomly chosen.(iii)
Buffer size (bi): 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 Section6.1. Furthermore, the Pareto front X is updated considering all these chromosomes. Let x be a chromosome (or solution); x is added to the Pareto front X if there is no other solution y∈X such that y dominates x. If x is added to X, then all solutions dominated by x are removed from X.


## 6.3.  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 Pc generating two offspring; (2) each offspring, or parent if the parents were not mated, undergoes mutation in accordance with the mutation probability Pm;(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.


## 6.4.  Crossover

The crossover operator receives one pair of chromosomes (P1 and P2), which are in the current population pop and have been randomly selected. The objective of this operator is to construct two offspring chromosomes (O1 and O2). For that, each time the crossover operation is performed, the following steps are made.(1)
Two cross points are randomly chosen,k1 and k2(1≤k1<k2≤N).(2)
Each gene in chromosomesP1 and P2 which is in position p, k1≤p<k2, is copied to the same position in chromosomes O1 and O2, respectively.(3)
Each gene in chromosomesP1 and P2 which is in position p, 1≤p<k1, is copied to the same position in chromosomes O1 and O2, respectively.(4)
Each gene in chromosomesP1 and P2 which is in position p, k2≤p≤N, is copied to the same position in chromosomes O1 and O2, respectively.Figure7 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.Figure 7
Crossover operation.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 (O1/O2) during steps 2 or 3, a new gene must be selected from the chromosome parent (P1/P2).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.


## 6.5.  Mutation

Mutation operation is performed on one chromosome, following these steps.(1)
Two positions (k1 and k2) of the chromosome are randomly chosen (1≤k1<k2≤N).(2)
Genes that are in positions betweenk1 and k2 (both included) are shuffled.(3)
The number of QCs in each gene located betweenk1 and k2 (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 betweenk1 and k2, both included, is modified by a random value that is between 0 and the average handling time hi*, of the vessel i in the same gene.Figure8 shows how the offspring oi, which has been obtained after the crossover operation, is mutated. First, two values k1=2 and k2=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 p, 2≤p≤4, are modified by selecting feasible random values for each one.Figure 8
An example of mutation operation.


## 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 (i, qi, and bi). Algorithm 2 describes how to modify a given solution cur in order to create a neighbor new. In this process, two operations are applied to solution cur:(i)
interchanging the position of two vessels (v1 and v2) in the chromosome, randomly chosen,(ii)
changing the values of number of QCs assigned (qi) and buffer size (bi) of a vessel i randomly chosen.Algorithm 2:createNeighbor.
Input:  cur: Actual solution/chromosome;
Output:  new: Neighbor solution/chromosome;
Copy chromosomecur into new
{Interchange the order of two vessels}
Vesselsv1 and v2 randomly chosen from new
Interchange position ofv1 and v2 in new
{Change the buffer assigned after its departure}
Randomly choose other vesselv from new
r←randInteger(0,2)
if  r=0  then
bv=0
else if  r=1  then
bv=hi*
else
if  rand(0,1)≤0.5  then
{Decrease a 10% the buffer size}
bv←max⁡(0,bv-hi*10)
else
{Increase a 10% the buffer size}
bv←min⁡(hi*,bv+hi*10)
end if
end if
if  rand(0,1)≤0.5  then
{Change number of assigned QCs}
qv←rand(2,QCv+)
end if
return Neighbour chromosome/solution newThis multiobjective simulated annealing algorithm (mosa function) is computed every k 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 frontX′ where the solutions from X have been improved to obtain a local optimal following the AMOSA scheme,(ii)
a new setSn 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 Sn set replace the solutions in population tmpPop 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.Algorithm 3:Clustering.
Input:  P: population;
Nm; maximum number of chromosomes;
Output:  P′: population with exactly Nm chromosomes;
{calculate the crowding distance for each element/chromosome inP}
crowdingDistance(P)
{sort the elements in ascending order}
sort(P)
{choose the biggestNm elements according to the crowding distance}
P′←Nm elements of P with the biggest crowding distances
return  P′


## 7. 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 (L) 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 Tables1(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.Table 1
Setting of the algorithms.
(a)
NSGAII and SPEA2+ settings
Parameter
MOGA scheme
Number of generations
500
Number of chromosomes
100
Mutation probability (Pm)
0.1
Crossover probability (Pc)
0.9
(b)
Multiobjective local search settings from MOGA + SA
Parameter
MOGA + SA
Initial temperature (Tmax⁡)
20
Minimum temperature (Tmin⁡)
0.001
Annealing factor (α)
0.9
Termination criterion
T
<
T
min
⁡
Cycle length
2CPLEX 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.Algorithm 4:Pareto front from the mathematical model.
Input:  I: Instance;
To: timeout;
Output:  X: Set of nondominated solutions;
Initialize set of solutionsS={}
for  λ∈{0,1} (steps of 0.1)  do
Solve the mathematical model for the instanceI and λ value given the timeout To
Add the schedule with the tuple of the objective functions (R^,T^s) to the set S
end for
X← nondominated solutions from S.Figure9 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.Pareto fronts of GA and CPLEX varying in the number of incoming vessels (|V|).
(a)
|
V
|
=
10 vessels
(b)
|
V
|
=
20 vessels
(c)
|
V
|
=
30 vessels
(d)
|
V
|
=
40 vesselsMultiobjective 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 p 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 (Pvalue<0.01).Table 2
Kruskal-Wallis test over the hypervolumes obtained in 5 different instances.
Instance
P value
Average hypervolume
NSGAII
SPEA2+
MOGA + SA
1
<2.2e − 16
0.175315
0.207921
0.236660
2
<2.2e − 16
0.174893
0.204749
0.252922
3
<2.2e − 16
0.178945
0.216683
0.261251
4
<2.2e − 16
0.180473
0.221138
0.258453
5
3.699e − 15
0.160139
0.173000
0.223231As 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, Table3 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.Table 3
Wilcoxon test over the hypervolumes obtained for a given instance.
NSGAII
SPEA2+
SPEA2+
0.00011
—
MOGA + SA
<2e − 16
<2e − 16Figure10 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 was able to produce a Pareto front with higher quality. Figures 10(a) and 10(d) show the relationship between R^ and T^s. MOGA + SA algorithm turned out to achieve schedules with greater robustness and lower T^s values (R^=1; T^s=0.3831) than NSGAII algorithm (R^=0.9227; T^s=0.4196). Furthermore, taking into account the relationship between T^s and R^ (see Figures 10(c) and 10(f)), MOGA + SA algorithm achieved schedules with lower standard deviation of robustness (σ(R^)=0.0) than the ones obtained by the NSGAII algorithm (σ(R^)=0.1557).Pareto fronts obtained by using or not local search.
(a)
NSGAII
(b)
NSGAII
(c)
NSGAII
(d)
MOGA + SA
(e)
MOGA + SA
(f)
MOGA + SAThe 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 vesseli is modeled as a delay d in the handling time of vessel i. This incidence is absorbed if there is enough buffer time behind vessel i 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.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 Table4(a)): the one with the minimum robustness (Rmin⁡), the one with the maximum robustness (RMax⁡), and one intermediate robust schedule (Ri where Rmin⁡<Ri<RMax⁡). 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)).Table 4
Percentages of incidences absorbed in schedules of 100 vessels.
(a)
Delays applied to schedules with different levels ofR^
Range
R
min
⁡
R
i
R
Max
⁡
d
∈
[
1
,
0.2
h
i
]
21.78
95.51
99.95
d
∈
[
1
,
0.5
h
i
]
18.73
93.58
99.85
d
∈
[
1
,
0.8
h
i
]
16.64
90.49
98.96
d
∈
[
1
,
1.0
h
i
]
13.92
87.10
98.01
d
∈
[
1
,
1.2
h
i
]
12.93
85.31
97.15
(b)
Delays applied to schedules with different levels ofT^s
Range
T
s
min
⁡
T
s
i
T
s
Max
⁡
d
∈
[
1
,
0.2
h
i
]
20.17
79.91
99.94
d
∈
[
1
,
0.5
h
i
]
16.16
73.88
99.75
d
∈
[
1
,
0.8
h
i
]
13.89
65.17
98.94
d
∈
[
1
,
1.0
h
i
]
11.85
60.25
98.08
d
∈
[
1
,
1.2
h
i
]
10.32
56.76
96.74
(c)
Delays applied to schedules with different levels ofσ(R^)
Range
σ
(
R
)
min
⁡
σ
(
R
)
i
σ
(
R
)
Max
⁡
d
∈
[
1
,
0.2
h
i
]
99.96
75.48
61.85
d
∈
[
1
,
0.5
h
i
]
99.95
72.23
54.42
d
∈
[
1
,
0.8
h
i
]
99.34
67.88
48.72
d
∈
[
1
,
1.0
h
i
]
98.87
65.34
47.77
d
∈
[
1
,
1.2
h
i
]
97.39
62.38
45.75The incidences (or delays,d) 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 (hi) 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.Table4(a) shows the percentage of incidences absorbed by each type of schedule. It can be observed that the more robust the schedule is, the more incidences absorbed. For instance, with delays d∈[1,0.5hi], the Rmin⁡ schedule only absorbed 18.73% of incidences in average, but the RMax⁡ schedule absorbed up to 99.85.%. Note that as the delay became larger, fewer schedules can absorb the incidences. With delays in the range of [1,0.2hi], the RMax⁡ schedule can absorb 99.95% of incidences in average. However, with larger ranges, the incidences absorbed decreased down to 97.15% in average. This pattern was also repeated in Table 4(b). The lower Ts, the lower incidences absorbed due to the fact that either there would not be buffers among vessels or there would be small buffers.Table4(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.8hi]. 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 Table5. 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.0hi], NSGAII was able to absorb 95.33% of incidences, whereas the MOGA + SA was able to absorb 97.34% of incidences.Table 5
Percentages of incidences absorbed in schedules of 100 vessels obtained using or not LS (timeout 30 secs).
Range
R
Max
⁡ no LS
R
Max
⁡ with LS
d
∈
[
1
,
0.2
h
i
]
99.53
99.88
d
∈
[
1
,
0.5
h
i
]
99.40
99.70
d
∈
[
1
,
0.8
h
i
]
97.62
98.51
d
∈
[
1
,
1.0
h
i
]
95.33
97.34
d
∈
[
1
,
1.2
h
i
]
94.04
96.00


## 8. 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.



---

*Source: 834927-2014-12-31.xml*