The quay crane scheduling problem with nonzero crane repositioning time and vessel stability constraints

Accepted Manuscript The Quay Crane Scheduling Problem with Vessel Stability and Nonzero Crane Moving time Noura Al-Dhaheri, Aida Jebali, Ali Diabat PII: DOI: Reference:

S0360-8352(16)30004-3 http://dx.doi.org/10.1016/j.cie.2016.01.011 CAIE 4234

To appear in:

Computers & Industrial Engineering

Received Date: Revised Date: Accepted Date:

1 May 2015 16 January 2016 18 January 2016

Please cite this article as: Al-Dhaheri, N., Jebali, A., Diabat, A., The Quay Crane Scheduling Problem with Vessel Stability and Nonzero Crane Moving time, Computers & Industrial Engineering (2016), doi: http://dx.doi.org/ 10.1016/j.cie.2016.01.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

The Quay Crane Scheduling Problem with Vessel Stability Consideration: Formulation and Heuristic Solution Approach

Noura Al-Dhaheri, Aida Jebali, Ali Diabat* Department of Engineering Systems & Management, Masdar Institute of Science & Technology, Abu Dhabi, UAE *Corresponding author e-mail [email protected]; Tel +971-2-810-9101; Fax +971-2-810-9901

The Quay Crane Scheduling Problem with Vessel Stability and Nonzero Crane Moving time

Abstract The continuous growth in worldwide container terminals’ traffic resulted in an increasing interest for the Quay Crane Scheduling Problem (QCSP) in research and practice. This problem consists of scheduling the discharge and load operations of the containers of a vessel by a set of quay cranes; the objective is to minimize the completion time in an attempt to increase container terminal throughput. In the literature, most of the proposed studies focus on improving model formulation and solution methods with a trend, in most recent papers, to incorporate more realistic features of the problem. Despite the importance given by practitioners to vessel stability in scheduling discharge and load operations, there is little research that considers this constraint. This paper presents a novel MIP formulation of the QCSP that takes into account vessel stability constraints. Furthermore, the proposed model is very flexible in handling various settings of the QCSP, such as those related to crane traveling time, task preemption and unidirectional quay crane operating mode. In order to tackle problem complexity, a Genetic Algorithm (GA) is proposed. Computational results validate the MIP formulation on small sized problems and highlight the performance of the proposed GA. Keywords: Quay Crane Scheduling, Container terminal, Mixed Integer Programming, Genetic algorithm

1.

Introduction

Globalization entails a steep increase in the containerized freight flow over the last decades (UNCTAD 2011) and this increase is expected to continue in the future. As a result, container terminals are continuously challenged to accommodate the growth of container volume by maximizing their operational efficiency, and investing in new handling technology and extended terminal infrastructure. This effort is additionally motivated by the need of container terminals to face the fierce competition with other terminals. The availability and the speed of service operations are henceforth the keys to attracting and retaining vessel operators. Providing container terminals with models and methods that lead to shortened vessel handling time and increase terminals’ throughput capacity is undeniably essential to help seaports respond to the incrementing container streams through the universal supply chain system and achieve a competitive advantage. Therefore, the last two decades witnessed an increasing number of research papers that aim to advance seaport operations with the use of quantitative methods. Recent classification of the existing models can be found in (Bierwirth & Meisel, 2010) that has been extended in (Carlo et al., 2013) and (Bierwirth & Meisel, 2015). Container terminals can be divided into four main areas: the berth, the quay, the yard and the gate. The berth and the quay areas are considered seaside, while the yard and the gate are considered landside (Carlo et al., 2013). Once a vessel is berthed, three main operations are processed: loading/unloading containers between vessels and landside trucks, transporting containers between berths and the storage yard and, loading/unloading containers between landside trucks and storage yard. The vessel is divided into several bays and each container is unloaded from or loaded onto a given bay according to the stowage plan; at any time, at most one QC can perform operations on a bay. The cost of constructing berths and handling equipment such as the quay cranes (QCs), used for loading and unloading containers from and onto vessels, is extremely high. Therefore, major emphasis in the literature has been placed on

optimizing the utilization of these two critical resources: the berths and the QCs. Three seaside operation problems are distinguished in the literature (Bierwirth and Meisel, 2010): (1) the Berth Allocation Problem (BAP), (2) the Quay Crane Assignment Problem (QCAP), and (3) the Quay Crane Scheduling Problem (QCSP). This paper tackles the QCSP. The latter aims at finding the schedule for the QCs serving a vessel that minimizes the total handling time. The QCSP includes both the assignment of container handling operations to specific QCs and the time schedule for these operations. The existing literature provides different ways of modeling the QCSP, as summarized by the classification scheme presented by (Bierwirth & Meisel, 2010). The authors classify the models according to four attributes: task, crane, interference and, performance measure. The task attribute specifies task definition and related constraints. Tasks of QCs are defined on the basis of container groups belonging to the same bay, bays, and bay areas (Meisel & Bierwirth, 2011). The crane attribute describes the assumptions made regarding QCs, such as their initial position and movement speed. The interference attribute is used to describe the spatial constraints of QCs’ movements, namely non-crossing and safety margin constraints. QCs move along the quay on a single rail track and are not allowed to cross one another. Moreover, a safety distance should be kept between adjacent QCs. The fourth attribute specifies the objective function of the QCSP. Carlo et al. (2013) extend this QCSP classification scheme to include further details on the QCSP, such as the safety margin for indented berths and performance measures related to the utilization of landside equipment. In the vast majority of QCSP models, tasks are defined on the basis of bays or container groups. Non-crossing constraints are commonly considered, while safety margin constraints and QC traveling time have been less involved. Task preemption is generally not allowed. Since QCSP with non-crossing constraints is an NP-hard problem (Guan et al., 2013), most of the proposed methods are still not fit to solve large-sized problem instances. Moreover, some realistic features of the QCSP have been seldom incorporated in the proposed models. For example, despite the importance given by practitioners to vessel stability in crane scheduling, there is little research that considers this constraint. In the current paper, we present a novel MIP formulation of the QCSP that takes into account vessel stability constraints. Furthermore, the proposed model captures various practical constraints of the QCSP, such as non-crossing, safety margin and crane traveling time constraints. In order to tackle problem complexity, a Genetic Algorithm (GA) is designed to solve medium and large sized problems. The remainder of this paper is organized as follows: section 2 reviews relevant literature on the QCSP with a focus on works which defining tasks on the basis of bays and container groups. Section 3 describes the addressed QCSP problem and the proposed mathematical model. Section 4 introduces the developed GA and the lower bound algorithm proposed to assess the performance of the GA. Section 5 reports computational experiments and results. Finally, section 6 concludes this work with important findings as well as directions and recommendations for future research on this topic. 2.

Literature Review

The crane scheduling problem while considering bays as tasks was initially studied by (Daganzo, 1989). The QCs are assigned to bays for certain time slots in order to minimize the total weighted completion times of ships. As preemption is allowed, more than one QC can perform operations on a bay. Peterkosfsky & Daganzo (1990) further investigate the problem and propose a branch-and-bound (B&B) solution approach. These two early studies do not however consider non-crossing and safety margin constraints. (Lim et al., 2004) and (Kim & Park, 2004) are the first papers that address spatial constraints. Lim et al. (2004) consider non-crossing constraints, safety distance between QCs and non-simultaneous constraints among bays. The objective is to find the QC-to-bay assignment that maximizes the throughput under these constraints. Dynamic programming algorithms, a probabilistic tabu search, and a squeaky wheel optimization heuristic are proposed to solve the problem. This work has been extended by incorporating task processing time in order to find the QC schedule that minimizes the makespan (Zhu & Lim, 2006; Lim et al., 2007). In (Zhu & Lim, 2006), the authors prove that by considering only noncrossing constraints, the problem becomes NP-complete. The problem is solved by a B&B algorithm and a simulated annealing heuristic. Lee et al. (2008) study the same problem, confirm its NP-completeness and solve it using a Genetic Algorithm (GA). In (Lim et al., 2007), the authors prove that for the QCSP with complete bays, there is an optimal schedule among the unidirectional ones; those where all QCs adopt, after their initial repositioning, the same movement direction. This seminal result entails a reduction of the search space to unidirectional schedules when seeking to optimally solve the QCSP with complete bays. Moreover, under the premise of unidirectional QC scheduling, crossing of QCs cannot happen, while the QCSP reduces to a QC-to-bay assignment problem. All the aforementioned studies do not consider QC traveling time constraints. In (Guan et al., 2103), the authors tackle the problem while taking into account non-crossing and QC traveling time constraints and assume that the partition of work units on bays ensures the required safety distance between QCs. A time-space network flow formulation is adopted and used to solve small-sized instances by standard solvers. For mediumsized instances, the authors develop a Lagrangian relaxation approach to obtain a tight lower bound and near-optimal solutions. Furthermore, to face the soaring complexity of large-sized problems, the authors develop two heuristics. In order to assess the quality of the proposed heuristics, a lower bound is determined based on the optimal solution of the preemptive QC schedule solved using a dynamic programming algorithm. Liu et al. (2006) study the QCSP for multiple vessels with different ready times. The objective is to minimize the maximum relative tardiness of vessel departures. The problem is broken down into two levels: the vessel-level and the berth-level. The vessel-level problem is a QCSP and is formulated as a MIP model that anchors the structure of unidirectional schedules. The latter takes into account initial crane positions, non-crossing, safety margin and QC traveling time constraints. The reduction of the research space to unidirectional schedules allows for solving non-trivial instances by a standard solver. Furthermore, the paper compares the preemptive and the non-preemptive QC scheduling approaches and highlights the improvements obtained by adopting the preemptive one. Kim & Park (2004) are the first who address the QCSP for container groups. They propose a model that considers precedence constraints among tasks, QC traveling time, non-crossing constraints and enforces safety distance by setting a non-simultaneity constraint between tasks located in adjacent bays. The objective is to minimize the weighted sum of the makespan and QCs completion times. A branch-and-bound algorithm and a Greedy Randomized Adaptive Search Procedure (GRASP) are proposed. This model has been later refined by Moccia et al.

(2006) in order to account for safety margin constraints in a more stringent way. A branch-andcut (B&C) algorithm is developed to solve medium and large sized instances. The algorithm succeeds in improving solutions for the benchmarks of Kim & Park (2004). To solve the same problem, Sammara et al. (2007) propose a Tabu Search (TS) algorithm where neighbourhoods are defined on the basis of a disjunctive graph. The computing time is significantly reduced at the expense of a slightly weaker solution quality with comparison to the B&C algorithm. (Bierwirth & Meisel, 2009) discloses the weakness of the model proposed in (Moccia et al., 2006) and proposes a model that incorporates a temporal distance between tasks in order to enforce non-crossing and safety margin constraints. Although there is not necessarily a unidirectional schedule among the optimal schedules for the QCSP with container groups, the authors propose a unidirectional search heuristic and show that it clearly outperforms the competing approaches (Kim & Park (2004); Moccia et al. (2006); Sammara et al. (2007)). Meisel (2011) introduces the QCSP with time windows and solves the problem using a tree-searchbased heuristic. In (Meisel & Bierwirth, 2011), the authors propose a unified approach to compare the performance of different model classes based on task granularity (container groups, bays, bay areas). The conducted experiments reveal that for some situations it is worth spending more time on solving a detailed model as this incurs a significant reduction in vessel handling time. Adversely, in other situations, simpler models are preferred. Legato et al. (2012) augment the model proposed in (Bierwirth & Meisel, 2009) by incorporating crane-individual service rates and a unidirectional QC travel mode. They employ a B&B scheme to solve the considered QCSP. Monaco & Sammara (2011) investigated the QCSP while considering time windows on QC availability, QC operational range, and unidirectional QC movement. A MIP formulation is proposed in order to minimize the makespan and QCs’ total delay. A commercial solver is used to solve small and medium size instances and a TS algorithm is designed to solve large instances. In (Chung & Choy, 2012) and (Kaveshgar et al., 2012) a GA is developed to solve the QCSP introduced in (Kim & Park, 2004). Chung & Choy (2012) show that their GA provides better solutions than other existing well-known approaches (TS, B&C, GRASP and B&B), especially for small size instances. For large size instances, even though the performance of the GA deteriorates, it remains attractive as it provides solutions within a short computational time. Kaveshgar et al. (2012) propose a GA using a solution representation that reduces the number of decision variables and employs tighter upper and lower bounds on the number of tasks assigned to a QC and, the set of tasks that could be done by a specific QC. The GA is tested on several benchmark instances proposed in (Meisel & Bierwirth, 2011). The results reveal that the proposed GA provides optimal or near optimal solutions in significantly shorter time with regard to the current best-known solutions. Guo et al. (2013) adapt the idea of generalized extremal optimization (GEO) to heuristically solve the QCSP. The obtained QC schedule is unidirectional and respects non-crossing and safety margin constraints. The heuristic provides optimal or nearoptimal solutions in reasonable time for large-sized problems. (Wang et al., 2013) is the first paper that highlights the importance of considering vessel stability in the QCSP. A GA is proposed to solve the QCSP introduced in (Kim & Park, 2004). The GA discards any QC schedule that violates the QCSP constraints, including those related to vessel stability. In all these works, pre-emption is not allowed. However, it is important to note that by viewing each container as a container group, the model of (Bierwirth & Meisel, 2009) can be used as a formulation of the pre-emptive QCSP. Thus, the pre-emptive QCSP can be assimilated to the QCSP with individual containers. This formulation is proposed for the QCSP allowing vessel bays to be served by more than one QC addressed in (Lu et al., 2012). The problem is solved

using a polynomial heuristic that generates contiguous and unidirectional schedules. The obtained schedules do not take into account QC ready time and initial position. Results demonstrate that the optimality gap is small for practical instances. In (Al-Dhaheri & Diabat, 2015), the authors propose a MIP formulation for the pre-emptive QCSP that consists in assigning QCs to bays over a planning horizon; the objective is to minimize the difference in workload between bays at all time periods. The model does not however take into account QC initial position, QC traveling time and safety margin constraints. It has been shown that optimal solutions can be obtained by a standard solver for small and medium sized problem instances. Some recent papers (Fu & Diabat (2014); Fu et al. (2014)) consider pre-emption while addressing the integrated quay crane assignment and scheduling problem. The problem is formulated as an IP and solved by a Lagrangian approach (Fu & Diabat, 2014) and a GA (Fu et al., 2014). In (Theodorou and Diabat (2014)), the authors use GA to solve an integrated quay crane assignment and scheduling problem, which accounts for crane positioning conditions. In (Theodorou and Diabat (2015)), they use Lagrangian relaxation approach to solve the problem and compare the performance of the new approach with the use of GA. However, all of the proposed models do not however take into account QC traveling time and initial position. The present paper formulates the pre-emptive QCSP for a single vessel, with an objective of minimizing the makespan. The workload of each bay may be performed by multiple QCs, thus the containers in each bay to be assigned to each QC is left as a decision. The model captures non-crossing, safety margin and QC traveling time constraints and takes into account QC initial position and vessel stability restrictions. The objective is to minimize the vessel’s handling time. Hence, the proposed MIP formulation extends the model proposed in (Al-Dhaheri & Diabat, 2015) by including QC initial position, QC traveling time, safety margin and vessel stability constraints. A GA is then proposed to solve medium and large sized instances. A lower bound is determined and used to assess the quality of the obtained solutions.

3. The QCSP: description and model formulation This section provides a description of the considered QCSP and the corresponding mathematical formulation. An example is then presented to illustrate the problem and the solutions obtained while exploring different search spaces. 3.1. Problem description The considered QCSP tackles a single vessel divided into a set of bays, J={1,…, . The bays are indexed sequentially from bay 1 to bay . The workload is expressed in number of work units, where a work unit is a considered fixed number of containers. The workload on the vessel can be processed by a set of identical QCs, I={1,…, . The QCs are indexed in an ascending order along the direction of increasing bay positions. Without loss of generality, bays and QCs are indexed hereafter in ascending order from left to right. Each task refers to container unloading or loading operations to be handled on a bay. Preemption is allowed since a single bay can be assigned to multiple QCs. Nevertheless, at any time, at most one QC can operate on a bay. Each QC has an initial bay position and is ready to commence service at time 0. All QCs are mounted on the same track and can move from one bay to another without crossing. The

traveling time, , between two consecutive bays is constant. The traveling time of QC i from bay j to bay j’, , is calculated as . A safety margin of bays is maintained between each two adjacent QCs i and i+1. In order to enforce non-crossing and safety margin constraints, the minimum temporal distance, formulated in (Bierwirth & Meisel, 2009), is adapted. A temporal distance should be incorporated in the schedule to separate the processing of containers located in bays j and j’ when processed by QCs i and i’ (ii’), respectively. represents the smallest allowed distance between the bay positions of QCs i and i’. The temporal distance to consider for all QCs i and i’ (ii’) when they respectively serve bays j and j’ is calculated as follows.

 if  if

and and

0 otherwise Note that refers to the case in which spatial constraints are respected by QCs i and i’ while serving bays j and j’. Let us denote by the temporal distance between the processing of containers located at bays j and j’ when they are processed by adjacent QCs i and i+1, respectively. As it can be noted, represents the smallest allowed distance between the bay positions of any adjacent QCs. Hence, the temporal distance to set between the processing of adjacent QCs depends only on the bays they serve. Furthermore, in order to respect non-crossing and safety margin constraints, it can be observed that each QC i can serve only bays j such that, +1 and . Let us denote this set for each QC i, by Ji={bi= +1, bi+1, …, ci= . It is also worth noting that all bays can be reached by at least one QC only under the following condition:  (+1). If the aforementioned condition is not satisfied and all bays contain workload then the considered QCSP is infeasible. The minimum number of bays needed to locate the QCs while accounting for safety margin is ( -1) ( +1)+1. For a number of bays =( -1) ( +1)+ with  , only bays 1,.., 1+, 1+(+1),…, 1+(+1)+, …, 1+( -1) (+1), …, 1+( -1).(+1)+ can be processed. The QCSP is hence infeasible unless all vessel operations are located on the aforementioned bays. Moreover, the workload of one bay can be shared between more than one QC only for and such that > ( +1). For  (+1), if the QCSP is feasible, then the QC-to-bay assignment can be easily determined. In this case, each bay can be served only by one QC. Hence, an optimal solution can be found by simply exploring unidirectional schedules. In this paper, the effort will be put into solving challenging QCSP instances where the workload of one bay can be shared between more than one QC. The stability of the vessel is violated if the vessel’s center of gravity (CG) shifts too much toward one side during the loading or unloading process. This shift results from workload distribution along the vessel and the sequence of operations. Henceforth the QCSP should consider vessel stability constraints to allow for obtaining QC schedules that can be used in practice. In particular, the proposed model ensures the longitudinal stability of the vessel

throughout the whole service process. The shift of the vessel’s CG is determined for each time unit based on the distribution of the finished work along the vessel and while assuming that initially the vessel’s CG is located at the point of bay length corresponding to . Nevertheless, the model can be easily accommodated to consider other possible initial positions of the vessel’s CG. 3.2. Model formulation The notation that is used in the proposed formulation is the following: Sets Set of bays, indexed by or where or is the bay number in ascending order from left to right Set of QCs assigned to the vessel, indexed by or where or is the QC number in ascending order from left to right Set of time units, indexed by t or Parameters Workload in number of work units at bay QC service rate in number of work units served per time unit Traveling time between bay and bay ′ Total weight of the vessel before starting operations Permitted shift of CG of the vessel (measured in 1 bay length) Average weight of each work unit Sufficiently large number Number of bays left for safety considerations between any two adjacent QCs 1 for loading operations; -1 otherwise (i.e. -1 for unloading operations) Initial position of QC i Temporal separation of processing bay j and bay j’ by two adjacent QCs (i and i+1) Decision Variables Makespan, latest completion time among bays Unhandled workload at time t in bay 1 If QC is assigned to bay at time segment 0 Otherwise 1 If QC starts moving from bay to bay at time 0 Otherwise 1 If bay is not completed at time 0 Otherwise The MIP formulation of the QCSP can be stated as follows: (1) Subject to:

(2) (3) (4) (5) (6) (7) (8) (9)





,



,



(10)

(11) (12) (13)



(14) (15)



(16) (17) (18)



(19) (20) (21) (22) (23)

The objective is to minimize the makespan of the schedule in order to provide the fastest possible service to the vessel and increase container terminal throughput. The makespan is determined in (2) by the latest completion time among all bays involved in the service. At this level, it is worth noting that the planning horizon must be an upper bound of the makespan. Constraints (3) and (4) guarantee that a given bay j is not completed at time t if and only if there is workload left in j at time t. Constraints (5) ensure that each QC can be assigned to at most one bay in any given

time unit, while constraints (6) guarantee that each bay is handled by at most one QC in any given time unit. Constraints (7) specify the initial workload for each bay. Constraints (8), (9) and (3) define the remaining workload at each bay for each time unit until it reaches zero, by removing the processed work units at the earlier time units. Constraints (9) and (4) guarantee that the remaining workload is set to zero for the remaining time units. Moreover, the combination of constraints (9) with the objective function ensure the completion of operations at every bay. Constraints (10) and (11) enforce a safe QC movement. They impose a temporal distance between adjacent QCs over each time unit. By taking into account all QCs positions at time unit t, they enforce QCs ordering condition and ensure a safety margin of  bays between all QCs. Constraints (12) restrict the assignment of QCs to only the range of bays that allow for the respect of the non-crossing and safety margin conditions. Constraints (13) specify that at any time unit t, QC i can start moving at most once from one bay j to another bay j’. Constraints (14) ensure that QC i can start moving from bay j to bay j’ at most once during any time interval  . Constraints (15) are related to the initial position of the QC. Constraints (16) ensure that a QC cannot be assigned to a bay while moving. Constraints (17) guarantee that a QC can start moving between two bays j and j’ at time unit t only if it is located at bay j at time segment t-1. A QC starts moving at time unit t only if the corresponding QC-to-bay assignment changes between t-1 and t (18). Constraints (19) account for the travel time between bays. Constraints (20) guarantee that all QCs do not start moving after the completion of all bays. The shift of the vessel’s CG is constrained to be within acceptable limits by constraints (21). Constraints (22) and (23) restrict the domains of the decision variables. The minimum time unit that can be used in the proposed QCSP model corresponds to the handling time of one work unit; and in this case, the QC service rate, , must be set to 1. The minimum work unit corresponds to one container. However, it must be noted that with the adoption of new spreader-technologies, QCs are able to simultaneously handle two or three containers (Bierwirth & Meisel, 2015), thus allowing for an increase of the minimum work unit to these aforementioned numbers. For this detailed time discretization, the accuracy of the proposed QCSP model increases at the expense of additional complexity and problem scale. Using a less detailed time discretization (and so a larger service rate) certainly leads to complexity reduction for the model, but it can entail the need to neglect the QC traveling time. The decision variable indicates if QC i is assigned to bay j at time segment t. Let us consider a work unit of one container and a QC handling rate of one container per time unit ( ). It is worth to note that in this case, also indicates if one container at the considered bay j is processed by QC i at time t. According to the classification scheme proposed in (Bierwirth & Meisel, 2010), and without accounting for vessel stability constraints, the QCSP addressed in this paper can be then written as containermove,poscoss,savemax(compl). Furthermore, it is worth to note that the precedence relations in processing pairs of containers at the same bay can be easily derived from the stowage plan and container grouping information. Based on these precedence constraints, the sequence of processing the containers at the same bay can be established. This sequence is then used to specify the container to be processed by QC i at bay j at time t according to the obtained QC schedule. The construction of this sequence is however beyond the scope of this paper.

3.3 Example A small QCSP instance is used to illustrate the proposed MIP solution. The example is then used to present solutions obtained for different options of search space restrictions. Three QCs are assigned to a vessel with 10 bays. As specified in table 1, a given number of containers have to be unloaded from each bay within a planning horizon of 40 time units. The adopted safety margin is equal to one bay. The initial position of QCs 1, 2 and 3 are bays 1, 5 and 10 respectively. The handling rate is set to 1, meaning that 1 time unit is required to process one container. Travel time between two adjacent bays is set to 1 time unit. Permitted shift of the vessel’s CG, total weight of all containers to be unloaded and weight of the vessel before starting operations are assumed to be of bay length 1, 57 and 200, respectively. On the one hand, according to constraints (3), M can be set to . On the other hand and based on constraints (9), M can be set to . Thus, hereafter, the maximum value among these two will be selected for M. In this QCSP instance, M is set to 15. The matrix of temporal distances , is given in table 2. Table 1. Vessel workload in work units Bay index j 1 2 3 4 Initial workload 8 4 0 15 Table 2. Matrix of temporal distances j \ j’ 1 2 3 1 2 1 0 2 3 2 1 3 4 3 2 4 5 4 3 5 6 5 4 6 7 6 5 7 8 7 6 8 9 8 7 9 10 9 8 10 11 10 9

4 0 0 1 2 3 4 5 6 7 8

5 0 0 0 1 2 3 4 5 6 7

5 3

6 0 0 0 0 1 2 3 4 5 6

6 0

7 10

7 0 0 0 0 0 1 2 3 4 5

8 5

9 0

8 0 0 0 0 0 0 1 2 3 4

10 12

9 0 0 0 0 0 0 0 1 2 3

10 0 0 0 0 0 0 0 0 1 2

For the considered problem, J1={b1=1, b1+1, …, c1=6}, J2={b2=3, b2+1, …, c2=8} and J3={b3=5, b3+1, …, c3=10}. Note that only a part of the matrix of temporal distances is needed to ensure safety movement between any adjacent QCs i and i+1, i=1..|I|-1. For example, to ensure safety movement between QC 1 and QC 2, we need to consider j in J1 and j’ in J2. The optimal schedule with Cmax=24, given by GAMS, is depicted in figure 1. Note that the integrality of the decision variables is relaxed while solving the problem. Indeed, given the following conditions: (1) the matrix defining the feasible region of is totally unimodular, (2)

and are binary; the linear relation of with and in constraints (13)-(20) will generate integer (0-1) solutions (Wolsey, 1998). Moreover, in order to strengthen the MIP formulation and speed up the solving time, valid cuts (24), related to the safety margin condition, are added. The latter imply that at any time, at most one QC can operate on any +1 consecutive bays.





(24) The optimal solution is obtained after 9149 seconds. The solved MIP contains 76444 constraints and 14002 variables among which there are 1600 discrete variables.

Figure 1. Optimal schedule obtained by GAMS

An optimal non-preemptive schedule for the same problem instance is presented in figure 2. The latter is obtained after 68 seconds by solving the proposed MIP to which we add the nonpreemption constraints (25). The MIP contains 77644 constraints. (25)

Figure 2. Optimal non-preemptive schedule

The obtained non-preemptive schedule does not represent an optimal solution of the original problem. Our example is enough to prove that better schedules can be obtained by preemptive scheduling even though non-preemption is a common assumption for the QCSP. Obviously, under the assumption of non-preemption, the computational time is significantly reduced.

Moreover, this example points out the interest for including preemptive schedules in the search space of the QCSP. A schedule is referred to as unidirectional if the QCs do not change the moving direction after the initial repositioning and have identical directions of movement either from upper to lower bays (referred to as right-to-left QC movement) or vice versa (referred to as left-to-right QC movement). Since unidirectional QC scheduling is the preferred policy in practice (Legato et al., 2012), two optimal unidirectional schedules are determined and presented in figures 3 and 4. The left-to-right schedule is obtained by solving the proposed MIP to which we add constraints (26), and the right-to-left schedule is obtained by solving the proposed MIP to which we add constraints (27). The two MIPs comprise 81709 constraints and are solved in less than 120 seconds.



(26)



(27)

Figure 3. Optimal left-to-right schedule

Figure 4. Optimal right-to-left schedule

As it can be observed, these two unidirectional schedules represent optimal solutions of the considered problem instance. In this example, reducing the search space to unidirectional schedules allows for finding an optimal solution to the original problem within reasonable computational time. Nevertheless, as we consider here a preemptive QCSP, reducing the search space to unidirectional schedules does not guarantee the generation of an optimal solution (Bierwirth and Meisel, 2009). However, it might represent a good strategy for solving the considered QCSP heuristically.

4. The Proposed Genetic Algorithm based heuristic The Genetic Algorithm (GA) is a search heuristic that simulates the process of natural evolution of living organisms to find near-optimal solutions for combinatorial optimization problems, including sequencing and scheduling. It combines the concept of survival of the fittest with controlled, yet randomized, information exchange to form robust exploration and exploitation of the solution space (Tavakkoli-Moghaddam et al., 2009). It starts by creating an initial population of chromosomes. Each chromosome encodes a solution of the problem and entails a fitness value derived from the value of the objective function for that solution. Then, during each iteration step (referred to as “generation”), genetic operators such as crossover and mutation are employed to breed offspring from two parent chromosomes, that are selected based on their level of fitness. The quality of the population keeps improving along with generations through the selection process of better fitted genes as occurs in the natural evolution of individuals (Goldberg, 1989). The GA is run until a suitable solution is found or a certain number of iterations (generations) is reached. The general structure of the GA is illustrated in figure 5. The steps of the proposed GA are presented with further detail below.

Figure 5. Genetic algorithm process

4.1. Chromosome representation Chromosome representation is the first step in the implementation of a GA. From our model formulation, we can see that a solution specifies: (1) the QC-to-bay assignment and, (2) the sequence of bays to be processed by each QC. The QC-to-bay assignment consists of determining the QCs to serve each bay. Since the workload of a bay is performed by more than

one QC, this decision is also made through the distribution of bays’ workload among QCs (i.e. the assignment of bays’ work units to QCs). Before introducing our chromosome representation, it is important to point out that the proposed GA embraces the strategy of searching a good solution to the considered QCSP within the space of unidirectional schedules. The sequence of bays to be processed by each QC is thus determined based on a unidirectional operating mode. In this case, the assignment of bays’ work units to QCs is the most complex decision (Lim et al., 2007). The GA chromosome is thus designed to generate these assignments. Accordingly, we can represent the chromosome by a matrix with |I|x|J| genes, where each gene (i,j) provides the part of the workload (in number of work units) of bay j assigned to QC i. As each QC i can perform operations only on bays j Ji, the gene (i,j) takes the value 0 for all i I and j J/Ji. Moreover, the sum of the value of genes (i,j) belonging to a same column j must be equal to the workload of bay j. The sum of values of genes (i,j) belonging to the same row i represents the total workload assigned to a QC i. Figure 6 illustrates an example of a chromosome representation based on the problem instance introduced in section 3.

Bay index QC1 QC2 QC3

1 8 0 0

2 4 0 0

3 0 0 0

4 5 10 0

5 2 0 1

6 0 0 0

7 0 6 4

8 0 3 2

9 0 0 0

10 0 0 12

19 19 19

12

57

Total workload to be handled 8

4

0

15

3

5

10

5

0

Figure 6. A chromosome representation

Furthermore, the genes (i,j), i I, j J/Ji can be omitted from the chromosome as their value is fixed to 0. Accordingly, the number of genes in the chromosome is reduced to . This results in the final chromosome representation depicted by figure 7. Bay index QC1 QC2 QC3

1 8

8

2 4

3 0 0

4 5 10

5 6 7 8 9 2 0 0 0 6 3 1 0 4 2 0 Total workload to be handled 4 0 15 3 0 10 5 0

10

12

19 19 19

12

57

Figure 7. Final chromosome representation

After assigning bays’ work units to QCs, the sequence of processing them is determined based on a unidirectional operating mode. The best left-to-right and right-to-left schedules, respecting the travel time and the non-interference constraints, are then derived from the QC-to-bay assignment represented by the chromosome by using the disjunctive graph model described in section 4.3.1. These schedules can indeed entail different makespans. Only one of the obtained schedules is retained for the considered chromosome, as it will be specified in section 4.3.2.

4.2. Generating initial population The initial population can be generated randomly or heuristically. Random generation gives the initial population great diversity. Heuristics generate individuals with better quality but may lead to a local optimal solution due to the limited diversity among the initial population. For the considered QCSP, the random method is used to generate the initial population of chromosomes. The algorithm to generate the initial population is described as follows. Let us denote by Qj the set of QCs that can serve bay j Qj ={i I / j Ji} Qj is a set of successive QC indexes. Let us denote by iminj and imaxj, the lowest and the highest QC index in Qj, respectively. For all j in J do Begin  If  Qj 1 then begin for all i in Qj/{imaxj} do begin randomly generate a discrete number j

between 0 and

 sij -sij

end end 

j

End 4.3. Fitness evaluation The fitness value assigned to a chromosome measures its performance. After the generation of a population, the fitness value of each chromosome must be calculated. The better the chromosome is, the higher its fitness value should be. Hence, chromosomes (individuals) with higher fitness values have more chances to survive. The chromosome used here specifies only the QC-to-bay assignment. Before evaluating its fitness, one schedule must be associated with the chromosome. For that, it is required to find out the best unidirectional schedules that can be derived from the chromosome. The left-to-right and the right-to-left unidirectional schedules that minimize the makespan, while respecting the travel time and the non-interference constraints, are determined using the disjunctive graph model described below. One of these schedules is selected and retained for the chromosome. 4.3.1 Disjunctive graph representation

The basic idea of using the disjunctive graph model is captured from the work of (Samara et al., 2007) and (Bierwirth and Meisel, 2009). In (Samara et al., 2007), the authors base the sequencing and the scheduling of QC tasks on a problem representation using disjunctive graphs. In (Bierwirth and Meisel, 2009), the authors propose a disjunctive graph model to generate unidirectional schedules from task-to-QC assignments that respect QC interference restrictions. Hence, our disjunctive graph model is rather adapted from the one proposed in (Bierwirth and Meisel, 2009). Moreover, it is important here to note that the setting of our QCSP is different of the one tackled in (Bierwirth and Meisel, 2009). In the considered preemptive QCSP, the number of work units that a QC performs on a bay is deemed as a decision. Nevertheless, in the QCSP addressed in (Bierwirth and Meisel, 2009), the group of containers (i.e. the set of work units or containers) that can be processed by a QC on a bay is predefined and is not part of the decision. Regardless of the direction of QC operations, a disjunctive graph model is used to build a unidirectional schedule for the QC-to-bay assignment that minimizes the makespan, incorporates travel time and respects the non-interference constraints. First, we describe how a left-to-right unidirectional schedule can be obtained from a given chromosome. Then, we specify how to obtain a right-to-left unidirectional schedule. In the disjunctive graph, all couples (i,j), i I and j Ji such that QC i processes work on bay j (i.e., the gene value corresponding to (i,j) in the considered chromosome is not null) are represented as nodes. A source node B and a sink node E are added to design the start and the end of the schedule. The sequence of bays that are processed by each QC i, according to a left-toright QC operation, is represented by a set of conjunctive arcs Ai which defines a path from node B to node E. Moreover, a pair of disjunctive arcs is used to link each pair of nodes (i,j) and (i’,j’) that can cause interference between QC i and QC i’. In other words, for each j Ji, for each j’ Ji’ such that , a pair of disjunctive arcs is used to link the pair of nodes (i,j) and (i’,j’). In each chromosome, the gene value corresponding to (i,j) indicates the number of work units of bay j processed by QC i. Let us denote by hij, the time needed to process these work units. The weights of conjunctive arcs are determined as follows:  The weight of the arc joining the source node B to node (i, bi), i I , is equal to  The weight of the arc joining node (i, ci), i I, to the sink node E is equal to  For all other conjunctive arcs in Ai , i I, the weight of the arc joining the nodes (i,j) to (i,j’) is equal to hij+tjj’ The weight of the disjunctive arc linking nodes (i,j) and (i’,j’) is equal to hij+ and, the weight of the disjunctive arc linking nodes (i’,j’) to (i,j) is equal to hi’j’+ . Note that these weights permit to enforce the required temporal distance for a safe QC movement. The obtained disjunctive graph is hereafter denoted by G=( , A, D, W) where  is the set of nodes, A is the set of conjunctive arcs, D is the set of disjunctive arcs and W represents arc weights, as described above.

From the disjunctive graph G, we construct the acyclic graph GL by selecting exactly one arc of each pair of disjunctive arcs. In order to obtain the left-to-right unidirectional schedule, GL retains the disjunctive arcs that are directed from nodes of the upper QC towards nodes of the lower QC. The makespan of the obtained left-to-right unidirectional schedule is determined by the length of the longest path of GL. A right-to-left unidirectional schedule can be obtained by considering the disjunctive graph G’=(, A’, D, W’) where A’ is the set of conjunctive arcs according to a right-to-left operating mode, and W’ represents arc weights, corresponding to the set of conjunctive arcs A’. Then, from the disjunctive graph G’, we construct the acyclic graph GR by selecting the disjunctive arcs that are directed from nodes of the lower QC towards nodes of the upper QC. The makespan of the obtained right-to-left unidirectional schedule is determined by the length of the longest path of GR. To determine the longest path of graphs GL and GR, we can use the Critical Path Method (CPM). Figures 8 and 9 respectively present the acyclic graphs GL and GR corresponding to the GA chromosome reported in figure 7. The makespan of the left-to-right unidirectional schedule is 30 and is given by the longest path of GL, (B, (3,5), (2,4), (2,7), (2,8), E). The makespan of the right-to-left unidirectional schedule is 28 and is given by the longest path of GR, (B, (2,8), (2,7), (2,4), (3,5), E).

Figure 8. Graph GL obtained for the chromosome of figure 7

Figure 9. Graph GR obtained for the chromosome of figure 7

4.3.2 Fitness evaluation

We can derive different schedules from graphs GL and GR since the non-critical nodes have a float time different than 0 and can start at any time between their earliest and latest start times. Hereafter, we consider for a chromosome k the two unidirectional schedules where QCs perform their operations on bays at their earliest start time (i.e. QCs do not stay idle while they can start processing a bay). Let us denote these schedules by SL and SR, respectively. If SL and SR violate vessel stability constraints, the chromosome k is discarded by setting its fitness value to 0. If vessel stability constraints are respected by SL and SR, then the schedule Sk retained for chromosome k is the one that minimizes the makespan. Otherwise, Sk is the feasible schedule among SL and SR. The fitness of a chromosome k can be set as follows:

4.4 Population replacement A new population is composed of old and new chromosomes. The old chromosomes are individuals of the previous population. The new chromosomes are the offspring produced by recombining some old chromosomes (parents). After offspring production, a population replacement strategy is needed to keep the population size constant. The most used strategies consist of removing the least fit chromosomes or the oldest individuals. The present paper employs the first strategy by maintaining in the population the chromosomes with the highest fitness value. 4.5 Stoppage rules We use a stoppage rule that it is commonly used by the GA: the maximum number of elapsed generations (Tavakkoli-Moghaddam et al., 2009). The process of creating a new population will be iterated until the maximum number of elapsed generations is reached. 4.6 Parent selection The roulette wheel is adopted to probabilistically choose the chromosomes that will serve as parents to produce offspring(s). The probability for a chromosome to be selected is given to each chromosome in the population according to its proportion of fitness. More specifically, a chromosome k in the population will have a high probability of being selected if it has a higher fitness, meaning that the corresponding schedule Sk achieves a lower makespan. 4.7 Genetic operators After two parent chromosomes are selected, they are mated to breed an offspring. The latter inherits characteristics from the parent chromosomes through crossover and mutation. In this paper, these genetic operators are designed to ensure that each bay’s workload is completely distributed among QCs. In other words, in any produced chromosome, the sum of the value of genes (i, j) in column j must be equal to the workload of bay j. 4.7.1 Crossover

An arithmetic crossover (AC) is used to produce new chromosomes in which the workload of each bay is fully shared between QCs. The AC produces new offspring as the linear combination of parents 1 and 2 as follows (Tavakkoli -Moghaddam et al., 2009): Offspring = x Parent 1 + (1- ) x Parent 2, where is a randomly generated number within the interval (0.5, 1) and Parent 1 refers to the parent with better fitness. This allows the offspring to inherit more characteristics from the parent with better fitness. Since values in genes (i,j) are integer values (number of work units of bay j processed by QC i), first, the integer part of the obtained gene value is retained. The total workload assigned to each QC is determined. For each bay j, if the assigned workload to QCs in Qj is less than , then (1) the QC i in Qj with the lowest total workload is selected; (2) the value in gene (i,j) is set to: and; (3) the total workload assigned to QC i is updated. If several QCs in Qj are associated with the lowest total workload, then the one with the highest index is selected. The produced offspring from two parent chromosomes of the example presented in section 2 using AC with = 0.8 is shown in figure 10.

Parent 1 Bay index QC1 QC2 QC3

1

2

3

4

8

4

0 0

5 10

8

5

6

Parent 2 7

8

9

2 0 0 0 6 3 1 0 4 2 0 Total workload to be handled 4 0 15 3 0 10 5 0

10

12 12

19 19 19

Bay index QC1 QC2 QC3

1

2

3

4

8

4

0 0

0 15

5

6

7

8

9

1 0 0 0 3 4 2 0 7 1 0 Total workload to be handled 4 0 15 3 0 10 5 0

8 57 Offspring Bay index 1 2 3 4 5 6 7 8 9 10 QC1 8 4 0 4 2 0 QC2 0 11 0 0 5 3 QC3 1 0 5 2 0 12 Total workload to be handled 8 4 0 15 3 0 10 5 0 12 Figure 10. Illustration of the arithmetic crossover

10

12

13 22 22

12

57

18 19 20 57

4.7.2 Mutation Crossover is used to combine parent genes in order to obtain a new chromosome. Mutation is then used to create a new chromosome by causing small perturbations in genes. The objective of mutation is to maintain diversity of the population and to better explore the solution space. This can be achieved by varying the number of consecutive bays served by each QC. The mutation for each chosen chromosome is: Step 1. Randomly choose a QC i and determine the set of bays served by i in the chosen chromosome. Select jmaxi and jmini, the highest and lowest index of bays served by QC i, respectively. Step 2. Randomly select one of the bays amongst jmaxi and jmini. Step 3. If the selected bay is jmini then go Step 4; else go to Step 5.

Step 4. (1) transfer the work units assigned to QC i at this bay to QC i-1; (2) select jmaxi-1, the bay with the highest index served by QC i-1 and; (3) transfer the work units assigned to QC i-1 at this bay to QC i. Step 5. (1) transfer the work units assigned to QC i at this bay to QC i+1; (2) select jmini+1, the bay with the lowest index served by QC i+1 and; (3) transfer the work units assigned to QC i+1 at this bay to QC i. It is worth to note that if in Step 1, QC 1 is selected, then only jmax1 can be considered in Step2. Similarly, if in Step 1, QC |I| is selected, then only jmin|I| can be considered in Step2. Figure 11 shows two examples of chromosomes that can be obtained by the proposed mutation for a chosen chromosome.

Bay index QC1 QC2 QC3

Original chromosome 2 3 4 5 6 7 8 9 4 0 4 2 0 0 11 0 0 5 3 1 0 5 2 0 Total workload to be handled 4 0 15 3 0 10 5 0

1 8

8 Example 1 Bay index QC1 QC2 QC3

1

2

3

4

8

4

0 0

11 4

8

5

6

7

8

9

10

Bay index QC1 QC2 QC3

1

2

3

10

12

18 19 20

12 57 Example 2 4

5

6

7

8

9

2 0 25 8 4 0 15 2 0 10 0 0 0 0 5 3 1 0 5 0 1 0 5 2 0 0 0 5 5 0 12 22 Total workload to be handled Total workload to be handled 4 0 15 3 0 10 5 0 12 57 8 4 0 15 3 0 10 5 0 Figure 11. Two examples of mutation operation

10

12

29 8 20

12

57

4.8 Optimality gap of the GA solution An optimality gap of the GA solution will be determined based on a combined lower bound, LB=max{LB1, LB2}, where LB1 and LB2 represent two lower bounds of the considered QCSP (without including any restriction on the search space). This is particularly needed to assess the quality of the GA solution for problem instances that cannot be optimally solved within reasonable computational time, and for which linear programming relaxation cannot provide a good lower bound. The lower bound LB1, is determined by the dynamic programming algorithm (LBA) presented below. LBA is adapted from the dynamic programming algorithm ALB, proposed in (Guan et al., 2013), by incorporating the range of bays that can be served by each QC. The foundation of

ALB stems from the observation that when preemption is allowed and safety margin and vessel stability constraints are not considered, in an optimal solution, each QC will serve consecutive bays (Guan et al., 2013). Under this scenario, ALB provides a lower bound for the preemptive QCSP and therefore a lower bound for the preemptive QCSP that considers safety margin and vessel stability constraints. By considering in LBA the range of bays that can be served by each QC, the lower bound provided by ALB can be improved as some of the schedules that do not respect the safety margin will be excluded from the search space. Nevertheless, this will not guarantee finding a solution that satisfies all safety margin constraints. The details of LBA can be described as follows: Index the work units in ascending order from left to right as the indices for QCs and bays. Thus we have ek  ek’ if k  k’, where ek represents the bay index of work unit k, k=1..K, with . Let represent a lower bound of the makespan for the first k work units served by the first i QCs. The initial dynamic programming value function can be obtained as follows:

 for k=1..K such that ekc1 , otherwise  for i=1.. I such that bi  e1  ci , otherwise The dynamic programming value function can be calculated recursively as follows:

 



for k=2..K and i=2.. I such that bi  ek  ci , otherwise Value function represents the time required for QC 1 to operate k work units located between bays b1 and c1. As work units located in bays beyond c1 cannot be processed by QC 1, this time is attributed a big value ( . represents the time required to process the first work unit by QC i and the first (i-1) QCs are idle. If QC i cannot process bay e1, this time is attributed a big value ( . represents the minimum time required by the first i QCs to process the k first work units. First, for each p=1..k, such that work unit p+1 is located in a bay within the operating range of QC i, we select the maximum of two parts: (1) the minimum time required for the first (i-1) QCs to serve the first p work units, (2) the minimum time required for QC i to travel and serve the rest of work units (from work unit p+1 to work unit k). The minimum time required by the first i QCs to process the k first work units is set to the minimum value among all the choices p=1..k such that bi  ep+1  ci.   represents the lower bound LB1 of our QCSP.

The idea of lower bound LB2 comes from (Lu et al., 2012). Given that in a feasible schedule at any time there can be at most one QC processing any +1 consecutive bays, then LB2= is another bound of the makespan. Conversely to LB1, it can be  observed that LB2 does not depend on the number of QCs.

5. Computational experiments The numerical experiments in this section are aimed towards evaluating the performance of the proposed MIP formulation and the effectiveness of the developed GA by solving problem instances of different sizes. The MIP models are solved by the commercial software GAMS. The GA and the LBA are coded in C++. All experiments are carried out on a PC Intel ® Core™ [email protected] GHz with the memory of 8.0 GB. In all instances, QCs are assumed to be equally located on the vessel. The traveling time of a QC between adjacent bays is set to one time unit and the safety margin is set to one bay. For each QC, the service rate is set to one work unit per time unit. The permitted shift of the vessel’s CG is assumed to be of one bay length. The average weight of a work unit is set to 1 and the weight of the vessel before starting operations is assumed to be around 4 times the total number of work units on the vessel. Small-sized problem instances are solved by GAMS. More specifically, the commercial software is used to find the optimal objective value of the unidirectional schedule (OVU) and wherever possible, the optimal objective value of the bidirectional one. For GAMS, the maximum computational time (CPU) is set to four hours. The same instances are also solved by the proposed GA. Based on the preliminary tests, the population size, the probability of crossover, the probability of mutation, and the maximum number of generations are set to 50, 0.8, 0.2, and 50, respectively. Each problem instance is solved 10 times and the best objective value (GA_OVU) and the average computational time (GA_CPU) are reported. For each unidirectional schedule, the QC operating direction (OPR_DIR) is flagged. If the optimal objective solution of an instance is not provided by GAMS within the computational time limit, the corresponding best lower bound is determined. This best lower bound is the largest value amongst LB and the lower bound provided by GAMS. The optimality gap of the unidirectional schedule obtained by GAMS is determined based on the optimal objective value or the corresponding best lower bound (OV_LB). Thirty small-sized instances are considered. The number of QCs ranges from 2 to 5 and the number of bays ranges from 8 to 12. For each number of QCs and bays, five problem instances, with different workload at each bay, are generated. The workload of each bay in number of work units is randomly generated from a uniform distribution of U(0,20). Table 3 reports the computational results for small-sized instances obtained from GAMS. An optimal unidirectional schedule is obtained for all the considered small-sized instances. As it can be observed from table 3, for 90% of the considered small-sized instances, GAMS returned the optimal unidirectional schedule in less than 30 minutes. Nonetheless, an optimal schedule has been found only for 19 instances and with an excessive computational burden. Furthermore, the

obtained results confirm the reduced complexity of the considered problem when a larger number of QCs serve a vessel. Similarly, they indicate a reduced complexity of the considered problem when the same number of QCs is allocated to a smaller vessel (with smaller number of bays). An optimal unidirectional schedule bears an optimal solution of the problem on 21 out of the 30 small-sized instances. Only for instances 13 and 30 the optimal unidirectional schedule is not an optimal schedule for the considered QCSP. For the other instances (4, 5, 20, 21, 22, 24 and 25), the optimality of the obtained unidirectional schedule cannot be confirmed and its optimality gap with regard to the optimal one is based on the best lower bound. Overall, the average gap between the optimal unidirectional schedule and the optimal one, is 0.85%. Therefore, the search for a unidirectional schedule seems to be a good strategy to find an optimal or a near-optimal solution to the considered QCSP. Moreover, it is worth mentioning that for 18 out of 30 instances, the proposed lower bound LB either equals the optimal objective value, or outperforms the lower bound provided by GAMS. More accurately, it equals the optimal objective value in 13 out of the 30 small-sized instances. Subsequently, based only on LB, the average optimality gap between the optimal unidirectional schedule and the optimal one, amounts to 2.61%, with a maximum value of 13.79%.

Table 3. Computational results for small-sized instances Instance Nb. Nb. OVU OPR_DIR CPU (s) no. QCs bays 1 36 LTR 17 2 34 LTR 19 3 2 8 37 LTR 11 4 50 LTR 102 5 55 LTR 408 6 26 RTL 3 7 31 LTR & RTL 3 8 3 8 33 LTR 5 9 34 LTR & RTL 4 10 37 LTR 8 11 35 LTR 1166 12 33 RTL 599 13 3 10 33 LTR 142 14 41 LTR & RTL 950 15 46 LTR & RTL 5185 16 31 RTL 31 17 30 LTR & RTL 19 18 4 10 33 LTR 24 19 35 LTR & RTL 57 20 36 LTR & RTL 65 21 32 RTL 1120 22 31 LTR 1385 23 4 12 33 LTR 233 24 36 RTL 2648 25 44 LTR 3253 26 31 LTR 57 5 12 27 30 LTR 62

OV_LB

GAP (%)

36 34 37 50 55 26 31 33 34 37 34 32 33 41 45 31 30 33 35 35 31 30 33 35 43 31 30

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.94 3.12 0.00 0.00 2.22 0.00 0.00 0.00 0.00 2.86 3.22 3.33 0.00 2.85 2.32 0.00 0.00

28 33 LTR 51 29 34 LTR & RTL 66 30 38 LTR & RTL 46 Notes: LTR stands for left-to-right; RTL stands for right-to-left

33 34 37

0.00 0.00 2.70

Table 4 reports the computational results for small-sized instances obtained by the GA and compares them to the ones obtained by GAMS in terms of objective value. GAP_1 measures the gap between the optimal unidirectional schedule provided by GAMS and the one obtained by the GA. GAP assesses the gap between the objective value of the GA solution and the optimal objective value.

Table 4. Computational results of the GA for small-sized instances Instance GA_OVU OPR_DIR GA_CPU (s) GAP _1 (%) no. 1 36 LTR 3.32 0.00 2 34 LTR 2.90 0.00 3 37 LTR 3.15 0.00 4 50 LTR 5.02 0.00 5 55 LTR 5.82 0.00 6 26 RTL 3.38 0.00 7 31 LTR & RTL 3.27 0.00 8 33 LTR 3.62 0.00 9 34 LTR & RTL 4.32 0.00 10 37 LTR 4.86 0.00 11 35 LTR 5.27 0.00 12 33 RTL 4.39 0.00 13 33 LTR 4.81 0.00 14 41 LTR 7.00 0.00 15 46 LTR & RTL 7.52 0.00 16 31 RTL 5.47 0.00 17 30 LTR & RTL 5.30 0.00 18 33 LTR 5.45 0.00 19 35 LTR & RTL 6.60 0.00 20 36 LTR & RTL 6.43 0.00 21 32 RTL 8.17 0.00 22 31 LTR 7.77 0.00 23 33 LTR 7.76 0.00 24 36 RTL 10.59 0.00 25 44 LTR 11.82 0.00 26 31 LTR 8.83 0.00

GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.94 3.12 0.00 0.00 2.22 0.00 0.00 0.00 0.00 2.86 3.22 3.33 0.00 2.85 2.32 0.00

27 28 29 30

30 33 34 38

LTR & RTL LTR RTL LTR

7.87 7.76 9.85 10.26

0.00 0.00 0.00 0.00

0.00 0.00 0.00 2.70

The entries of table 4 show that the GA returns the optimal unidirectional schedule for all smallsized instances. Thus, based on the results of table 3, the GA achieves optimal solutions for 21 out of the 30 small-sized instances. For most of the small-sized instances, the computational time to obtain the GA solution (the average computational time to obtain a GA solution is 6.28 seconds) is markedly lower than the one spent to obtain the GAMS solution (the average computational time to obtain GAMS solution is 591.30 seconds). This indeed highlights the major advantage of the proposed GA. Thirty medium and large size instances are moreover generated and solved by the GA. The number of QCs ranges from 2 to 6 and, the number of bays ranges from 15 to 25, as indicated in table 5. For each number of QCs and bays, five problem instances, with a different workload at each bay, are generated. The workload of each bay in number of work units is randomly generated from a uniform distribution of U(0,100). In these experiments, the genetic parameters are maintained similar to those used for small-sized problem instances, except for the number of generations. The latter is tuned based on problem complexity as revealed by preliminary tests conducted on instances characterized by the same number of QCs and bays. Hence, for problem instances 31-35 and 41-45, the number of generations is set to 50; for problem instances 36-40 and 51-55, the number of generations is set to 100 and; for problem instances 46-50 and 56-60, the number of generations is set to 150 and 200, respectively. Table 5. Computational results of the GA for medium and large size instances Instance Nb. Nb. GA_OVU OPR_DIR GA_CPU (s) LB GAP (%) no. QCs bays 31 367 LTR 240.93 367 0.00 32 406 LTR 339.79 406 0.00 33 2 15 411 LTR 348.62 410 0.24 34 314 LTR 199.67 314 0.00 35 365 LTR 275.61 365 0.00 36 185 RTL 290.25 184 0.54 37 209 RTL 335.44 203 2.95 38 4 15 207 RTL 438.67 205 0.97 39 158 LTR 224.59 156 1.28 40 186 LTR 277.17 183 1.63 41 333 LTR 395.85 327 1.83 42 367 LTR 503.82 358 2.51 43 3 20 366 LTR 505.73 358 2.23 44 312 LTR 366.25 305 2.29 45 337 LTR 428.87 330 2.12 46 166 RTL 827.90 164 1.21 47 186 LTR 846.94 179 3.91 48 6 20 184 RTL 1256.50 179 2.79 49 176 LTR 674.32 173 1.73 50 171 LTR 751.93 165 3.63 51 476 RTL 1852.94 462 3.03 3 25 52 475 LTR 1997.98 462 2.81

53 54 55 56 57 58 59 60

6

25

463 363 431 339 243 232 203 223

LTR RTL LTR LTR & RTL LTR LTR & RTL RTL RTL

1927.99 1112.63 1684.86 2258.25 2735.20 2832.71 1948.35 2181.67

448 351 419 231 232 225 177 210

3.34 3.41 2.86 3.46 4.74 3.11 14.68 6.19

It is worth at this level to note that all lower bounds LB are obtained in a very short time (less than 1 second). Moreover, linear programming (LP) relaxation is used in an attempt to determine a LP lower bound and compare it to LB. Different experiments on small-sized problem instances point out that the best lower bound can be obtained by relaxing the integrality of the decision variables ( ). However, for 40% of the instances, this LP bound cannot be obtained within the computational time limit (4 hours). This relaxation is moreover performed on problem instance 39. Once again, the LP bound is not provided by GAMS within the computational time limit and the best lower bound returned by the solver is 83. This result highlights the significance of LB to evaluate the quality of the GA solution.

Figure 12. GA performance

For the considered medium and large sized problem instances, GAMS is not able to reach a solution within the time limit. However, the GA delivers solutions for these problem instances within reasonable computational time (the average computational time is less than 30 minutes for 22 out of the 30 problem instances). As it can be observed in figure 12, on average, both the computational time and the optimality gap increase when the number of bays increases and/or the number of QCs increases. More computational effort is indeed needed for instances of larger sizes to better explore the unidirectional search space. For the most intractable instances (56-60), the computational time of the GA does not exceed 50 minutes and the quality of the obtained solutions is acceptable. These results can only favor the recommendation of the proposed GA for use, particularly, in medium and large size problems.

6. Conclusion In this paper, we propose a novel MIP model for the QCSP that takes into account several practical features of the problem such as preemption, non-crossing, safety margin, QC traveling time, QC initial position and vessel stability. Moreover, the model can be slightly modified to

search for optimal unidirectional and/or non-preemptive schedules. The MIP formulation provides an optimal unidirectional schedule for all considered small-sized problems, in relatively reasonable time. However, it does not return the optimal solution for all these instances. Moreover, the experimentation of the MIP formulation on small-sized problems and the use of the proposed lower bound (LB), highlight the interest for adopting a search strategy focused on unidirectional schedules. A Genetic Algorithm (GA) is also designed to solve the problem. The GA embraces the unidirectional search strategy. It provides optimal unidirectional schedules for all small-sized instances within significantly lower time than the one required by GAMS. Furthermore, for medium and large size instances, where GAMS fails to solve the problem, the GA returns a nearoptimal solution within a reasonable computational time. This promotes the use of the proposed GA as a solution approach for the considered QCSP, especially for medium and large sized instances. For future work, additional circumstances should be considered, such as variation in QC service rate and in container weight. Also, the MIP model can be extended to account for loading and unloading operations and incorporate ready time per QC. Another prospect of this work would consist in developing other solution methods for the sake of generating better QC schedules.

References Al-Dhaheri, N., & Diabat, A. (2015). The quay crane scheduling problem. Journal of Manufacturing Systems, 36, 87-94. Bierwirth, C., & Meisel, F. (2010). A survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research, 202, 615-627. Bierwirth, C., & Meisel, F. (2015). A follow-up survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research, 244(3), 675-689. Bierwirth, C., & Meisel, F. (2009). A fast heuristic for quay crane scheduling with interference constraints. Journal of Scheduling, 12, 345-360. Carlo, H. J., Vis, I. F. A., & Roodbergen, K. J. (2013). Seaside operations in container terminals: literature overview, trends, and research directions. Flexible Services and Manufacturing Journal, 27(2-3), 224-262. Diabat, A., & Theodorou, E. (2014). An integrated quay crane assignment and scheduling problem. Computers & Industrial Engineering, 73, 115-123. Fu, Y-M., & Diabat, A. (2014). A Lagrangian relaxation approach for solving the integrated quay crane assignment and scheduling problem. Applied Mathematical Modelling, 39(3), 1194-1201. Fu Y-M., Diabat, A., & Tsai, I-T. (2014). A multi-vessel quay crane assignment and scheduling problem: Formulation and heuristic solution approach. Expert Systems with Applications, 41, 6959-6965. Guo, P., Cheng, W., & Wang, Y. (2014). A modified generalized extremal optimization algorithm for the quay crane scheduling problem with interference constraints. Engineering Optimization, 46(10), 1411-1429. Chung, S. H., & Choy, K. L. (2012). A modified genetic algorithm for quay crane scheduling operations. Expert Systems with Applications, 39, 4213-4221. Daganzo, C. F. (1989). The crane scheduling problem. Transportation Research Part B, 23 (3), 159-175. Goldberg, D. E., Genetic algorithms in search optimization and machine learning, 1989, Addison-Wesley Longman Publishing Co., Inc. Boston, MA, USA. Guan, Y., Yang, K-H., & Zhou, Z. (2013). The crane scheduling problem: models and solution approaches. Annals of Operations Research, 203, 119-139.

Kaveshgar, N., Huynh, N., & Rahimian, S. K. (2012). An efficient genetic algorithm for solving the quay crane scheduling problem. Expert Systems with Applications, 39, 13108-13117. Kim, K. H., & Park, Y. M. (2004). A crane scheduling method for port container terminals. European Journal of Operational Research, 156, 752-768. Lee, D. H., Wang, H. Q., & Miao, L. (2008). Quay crane scheduling with non-interference constraints in port container terminals. Transportation Research Part E, 44, 124–135. Legato, P., Trunfio, R., & Meisel, F. (2012). Modeling and solving rich quay crane scheduling problems. Computers & Operations Research, 39, 2063-2078. Lim, A., Rodrigues, B., Xiao, F., & Zhu, Y. (2004). Crane scheduling with spatial constraints. Naval Research Logistics, 51, 386-406. Lim, A., Rodrigues, B., & Xu, Z. (2007). A m-parallel crane scheduling problem with a non-crossing constraint. Naval Research Logistics, 54, 115-127. Liu, J., Wan, Y-w., & Wang, L. (2006). Quay crane scheduling at container terminals to minimize the maximum relative tardiness of vessel departures. Naval Research Logistics, 53, 60-74. Lu, Z., Han, X., Xi, L., & Erera, A. L. (2012). A heuristic for the quay crane scheduling problem based on contiguous bay crane operations. Computers & Operations Research, 39, 2915-2928. Meisel, F., & Bierwirth, C. (2011). A unified approach for the evaluation of quay crane scheduling models and algorithms. Computers & Operations Research, 38, 683-693. Meisel, F. (2011). The quay crane scheduling with time windows. Naval Research Logistics, 58, 619-636. Moccia, L., Cordeau., J-F.,Gaudioso, M., & Laporte, G. (2006). A branch-and-cut algorithm for the quay crane scheduling problem in a container terminal. Naval Research Logistics, 53, 45-59. Monaco, M. F., & Sammara, M. (2011). Quay crane scheduling with time windows, one-way and spatial constraints. International Journal of Shipping and Transport Logistics, 3 (4), 454-474. Peterkosfsky, R. I., & Daganzo, C. F. (1990). A branch and bound solution method for the crane scheduling problem. Transportation Research Part B, 24 (3), 159-172. Sammara, M., Cordeau., J-F., Laporte, G., & Moccia, L. (2007). A tabu search for the quay crane scheduling problem. Journal of Scheduling, 10, 327-336. Tavakkoli-Moghaddam, R., Makui, A., Salahi, S., Bazzazi, M., & Taheri, F. (2009). An efficient algorithm for solving a new mathematical model for a quay crane scheduling problem in container ports. Computers & Industrial Engineering, 56, 241-248. Theodorou, E., & Diabat, A. (2015). A joint quay crane assignment and scheduling problem: formulation, solution algorithm and computational results. Optimization Letters, 9(4), 799-817. UNCTAD, Review of maritime transport 2011. United Nations, New York and Geneva, 2011. Wang, J., Hu, H., & Song, Y. (2013). Optimization of quay crane scheduling constrained by stability of vessels. Transportation Research Record: Journal of the Transportation Research Board, 2330, 47-54. Wolsey, L.A., Integer Programming, 1998, John Wiley and Sons, New York. Zhu, Y., & Lim, A. (2006). Crane scheduling with non-crossing constraint. Journal of the Operational Research Society, 57 (12), 1464-1471.

The Quay Crane Scheduling Problem with Vessel Stability Consideration: Formulation and Heuristic Solution Approach Research Highlights: 1. 2. 3. 4.

Present a novel formulation to the QCSP. We consider vessel’s stability that nobody before has considered. We develop an efficient GA to solve the problem. GA compared with lower bounds derived using dynamic programming.