An Improvement of Branch-and-Price Algorithm for Quay Crane Scheduling Problem in Single Ship Setting

(1)

1

An Improvement of Branch-and-Price Algorithm for

Quay Crane Scheduling Problem in Single Ship Setting

Master’s Thesis SCM MSc Supply Chain Management Faculty of Economics and Business

University of Groningen The Netherlands

Student name: Anisha Maharani Student number: S2540096 email: a.maharani@student.rug.nl

Supervisor: dr. Evrim Ursavas

(2)

2

Abstract

(3)

3

1 Introduction

Due to globalization, an increase in number of container traffic has emerged in recent decades. As a result, the number of port container terminals is growing as well. Nguyen, Zhang, Johnston, and Tan (2013) specified that container terminals play a key role as an interface between land and sea transportation. In order to accommodate this trend, an investment for new equipment and infrastructure of container terminals is needed. On the other hand, container terminals also need to reduce their costs in order to stay competitive. Having said that, those terminals should make use of their resources by developing an effective and efficient operations management.

One of the performance indicators to measure the competitiveness of a port container terminal is the vessel turnaround time, which is the average time of a vessel’s stay in the terminal (Chen, Lee and Goh, 2013). Exposito-Izquierdo, Gonzalez-Velarde, Melian-Batista and Moreno-Vega (2013) stated that the main cause of large turnaround time is the high rate of utilization of a terminal’s infrastructure. Therefore, the most effective way to reduce the turnaround time is to improve the productivity of resources, e.g. berths and handling equipment. One of the important handling equipment of a port container terminal is the Quay Crane (QC). It has the function to handle the containers of berthed vessels. Being one of the main bottlenecks of container terminals, QC performance can largely affect a terminal’s throughput and efficiency (Diabat & Theodorou, 2014). For that reason, this research explicitly focuses on the quay crane scheduling problem (QCSP). According to Nguyen et al. (2013), the objective of QCSP is to determine a sequence of handling operations of a vessel with a given number of QCs in order to minimize the vessel completion time.

From mathematical point of view, QCSP is considered as a complex problem thus resulting to long computational time (Chen et al., 2013). Yet, Choo, Klabjan and Simchi-Levi (2010) emphasized that the decision of QCSP should be obtained in a very short time, because the required inputs for solving the problem get known just before the arrival of the ship. Large-scale optimization techniques are needed to cope the complexity of the QCSP and to provide the optimal solution. There are two alternative for solving QCSP, by using a heuristic procedure or an exact solution algorithm. Recent researches about QCSP opt more to develop heuristic procedures because they perform faster. Several heuristic approaches have been conducted to tackle the QCSP in a limited amount of time, e.g. the LTM algorithm to solve rich-QCSP (Legato, Trunfio and Meisel, 2012), tabu search algorithm (Sammara et al., 2007), simulation (Legato, Mazza and Trunfio, 2008), constraint programming (Unsal and Oguz, 2013), and genetic algorithm (Hakam, Solvang and Hammervoll, 2012; Kaveshgar, Huynh and Rahimian, 2012; Nguyen et al. 2013; Diabat and Theodorou, 2014). Despite their fast computational time, heuristic approaches have disadvantage since they cannot guarantee an optimal solution.

(5)

5

solution and outperform heuristic approaches. Nevertheless, using those exact solution approaches to solve large-scale QCSP will result to longer computational time.

One exact solution algorithm that fits to solve large-scale linear problem as the QCSP is branch-and-price (B&P). B&P algorithm is suitable to handle complex problems because it only considers a subset of variables at each iteration. It is expected to reduce the complexity of the problem whilst still provide an optimal solution (Danna and Pape, 2005). At the same time, several characteristics of B&P can result to high solution times and reduce the efficiency of the algorithm, especially in large scale instances (Dzubur and Langvik, 2012).

The main aim of this research is to develop a novel B&P algorithm by implementing the several improvement techniques in order to expedite the computational time for solving the QCSP. To the best of our knowledge, acceleration techniques for the B&P algorithm have not been developed for the QCSP setting so far. Therefore, the research question of this paper is: “how can the B&P

algorithm be improved in order to solve the QCSP such that global optimal solution can be obtained instantly? “.

The main method is to extend the study of Choo et al. (2010) by developing a new B&P algorithm with several improvement techniques provided in the literature. The result then will be compared to that obtained from other algorithms presented in the research by Choo et al. (2010) in terms of computational time, optimality gap, and other performance measures. Afterwards, a validation is conducted by implementing the advanced B&P algorithm to a new-developed model by Boulaniki (2015). This model formulates an integrated model built upon the model by Choo et at. (201) and solves the QCSP concurrently with the yard storage location assignment.

The theoretical contribution of this research is to improve the efficiency of the B&P algorithm in order to solve the QCSP instantly. Application to an extended model will ensure the robustness of the new B&P algorithm. This has important practical contribution since as mentioned before the terminal only has a little amount of time to make a decision of the QCSP, so a fast algorithm is needed in order to handle this situation.

The remainder of this paper is organized as follows. The second chapter will give an overview of container terminals operations. Related literature about the QCSP is described in chapter three. The fourth chapter will present the initial model to solve the QCSP and also the extended model of integrated approach. Chapter five will present the description of the solution algorithm and a number of new techniques that will be applied to improve the B&P algorithm. Some numerical and computational results will be given in chapter six. The last chapter will conclude this paper and provide some suggestions for future research.

2 Overview of Container Terminal Operations

(6)

6

storage locations. The key function of container terminals is to receive the outbound containers from shippers to be loaded onto vessels and to discharge inbound containers from vessels to be picked up by consigners (Choo et al., 2010). Container terminals are high-capital intensive and several specialized equipment are used to handle the containers. Thereby, an effective and efficient management of port container terminal operations becomes crucial.

Figure 1 Schematic representation of a container terminal system. Adapted from: Lajjam et al. (2013) The process of handling the containers is demonstrated as follows. When a vessel arrives in a seaport, a vessel is assigned to a berth for loading and unloading the containers. The assignment also indicates the required number of QCs to serve that vessel and which ship bay is served by which QC. Sometimes, the previous decision is solved simultaneously with scheduling the arrangement of QCs in each ship (also sometimes called as quay crane scheduling problem/QCSP) because they are highly connected. The main objective of the QCSP according to Nguyen et al. (2013) is to determine the sequence of handling containers in order to minimize the vessel’s stay time given the number of cranes assigned to a vessel. The output of the QCSP is critical since the productivity of the handling equipment highly affects the turnaround time of berthed vessels (Exposito-Izquierdo et al., 2013).

The next process is transporting containers from the quay-side to pre-determined storage locations in the yard. Afterwards, yard cranes (YCs) lifts the containers at the landside from the trucks onto their assigned stacks, the trucks then go back to be used for other jobs. These two last processes are directly affected by the decision of the QCSP, thus it is critical to obtain the optimal solution for it (Chen et al., 2013). Occasionally, the yard storage location assignment for each container is determined at the same time with the QCSP, because the storage location of specific container influences the utilization and the travel distance of the handling equipment (Zhang et al., 2003).

(7)

7

(Diabat & Theodorou, 2014). Terminal productivity and efficiency are measured by the number of containers handled by QCs, so an improvement of QC efficiency will result in improvement of overall terminals throughput, including reduced vessel turnaround time.

Figure 2 Top view of a vessel. Adapted from: Diabat and Theodorou (2014)

Diabat and Theodorou (2014) indicated in their research two important aspects of the QCSP, which are enforcing the position conditions at all times and considering some spatial constraints, e.g. non-crossing constraints and clearance constraints. In certain cases dealing with multiple vessels, yard congestion constraints also are considered as an important factor to ensure that there will be no traffic in the yard storage location at any point of time. A set of tasks in the QCSP represents handling operations for a vessel given a set of assigned QCs (Bierwirth and Meisel, 2010). In this research, a task is defined as handling a single container. This definition gives more detail schedule and clearly represents the real life situation in such way a QC handles only one container at one specific time.

3 Literature Review

Current studies about the QCSP are trying to find the best algorithm in order to obtain the optimal solution as fast as possible. This is needed since the inputs for solving this problem are only known in a short time, just before the arrival of a ship (Choo et al., 2010). Using a standard commercial MILP solver formerly preferable since it will give the final solution directly. However, the QCSP can be considered as a complex problem hence using a commercial solver can result to long computational time or infeasibility due to the lack of memory.

The QCSP has been studied before in many different settings. The first studies of the QCSP are conducted by Daganzo (1989) and Peterkofsky and Daganzo (1990), where an exact branch-and-bound algorithm is developed to solve mixed integer programming (MIP) for the loading of ships and assigning cranes to bays. The objective is to minimize the aggregate cost of delay incurred on the vessels. However, both studies did not take into account the spatial constraints in their model.

(8)

8

and Laporte (2006). They added clearance constraints to the model and used an exact branch-and-cut (B&C) algorithm to solve the problem. Based on Moccia et al. (2006), Sammara, Cordeau, Laporte, and Monaco (2007) proposed a tabu search approach to solve larger scale problems for practical applications. Nonetheless, there are some drawbacks in the model developed my Sammara et al. (2007) found by Bierwirth and Meisel (2009). The latter research proposed the idea to insert a minimum temporal distance between tasks to fix the drawbacks of the previous model and they used a branch-and-bound algorithm as the core to solve the problem.

Current research streams in the QCSP are more likely to develop a meta-heuristic or heuristic procedure since the result can be obtained in a fast way. One of the heuristic approaches is the LTM algorithm, which developed by Legato, Trunfio and Meisel (2012). This algorithm is applied to solve rich-QCSP. Previously, Legato, Mazza and Trunfio (2008) used simulated annealing and adaptive balance explorative and exploitative search to minimize the total completion time of scheduling. Different approaches called constraint programming (Unsal and Oguz, 2013) and hybrid estimation of distribution algorithm (Exposito-Izquierdo et al., 2013) are developed lately to get good solutions in a little amount of time compared with other heuristics. Recently, solving the QCSP using genetic algorithm (GA) approach is prominent since this approach can efficiently find optimal or near optimal solution. Several researches using genetic algorithm to solve QCSP are conducted by Lee, Wang, and Miao (2008); Hakam, Solvang, and Hammervoll (2012), Chung and Choy (2012); Kaveshgar, Huynh and Rahimian (2012); Nguyen et al. (2013), and Diabat and Theodorou (2014). In the research by Hakam et al. (2012), they proved that a GA can achieve results with a smaller optimality gap and faster computational times compared to the best known solution obtained by other approaches such as tabu search and branch-and-cut (B&C).

As a conclusion, heuristics are widely employed to tackle the complexity of the QCSP (Diabat and Theodorou, 2014). At the same time, the global optimum solution is still more preferable in practice since the output of QCSP will affect the overall terminal operations. In this regards, heuristic approaches are less attractive when compared to exact solution algorithms. In order to obtain the optimal solution, one of the exact solution algorithms such as branch-and-bound, branch-and-cut, and branch-and-price can be used. However, using these exact solution approaches to solve large-scale the QCSP will result to long computational times. Accordingly, there is a need for developing an advanced algorithm which can find the optimum solution for the QCSP in a very short time. This can be done by further development of existing exact algorithm to accelerate the computational times.

One of the exact solution algorithms which are considered to be applicable to solve the QCSP is branch-and-price algorithm (B&P). B&P algorithm in container terminals operations is more developed for other frameworks, such as integration of berth allocation and quay crane allocation problem (Vacca, Salani an Bierlaire, 2013) or integration of berth allocation and yard assignment problem (Robenek et al., 2014). In these studies, it is shown that the B&P algorithm is better than heuristics approaches and commercial solvers in terms of computational time and the optimality gap.

(9)

9

problem with huge number of variables (Barnhart et al., 1998). To the best of our knowledge, only few researches use this algorithm to solve the QCSP because of its complexity. One of them is the research by Choo, Klabjan, and Simchi-Levi (2010). This research proposed a MIP model to solve both single-ship and multi-ship problems with respect to spatial constraints. Choo et al. reformulated the single-ship problem into a set covering problem and used branch-and-price to solve the problem. A combination of Lagrangian relaxation and branch-and-price algorithm is formulated to solve the multi ship problem. It is shown that the result of B&P algorithm is better compared to that of other approaches.

Although B&P is suitable to solve large problem as the QCSP, this algorithm has drawbacks because of its characteristics which will result to long computational time. Thereby, several improvement techniques are needed in order to maximize the utilization of the B&P algorithm for solving the QCSP. Literatures about B&P have already developed certain methods to improve the performance of B&P algorithm. However, those methods have not been applied in the QCSP scheme so far. This research aims to extend the model of Choo et al. (2010) by implementing the improvement techniques of B&P algorithm in order to obtain the exact optimal solution for the QCSP in a fast way. The implementation of the new B&P algorithm to a new-extended QCSP problem to clarify its performance. Chapter five will provide more detailed explanations about the B&P algorithm and the available improvement approaches in the literature.

4 Mathematical Model

In this study, we considered the quay crane scheduling problem in single ship situation. The model and the solution algorithm are built upon the Choo et al. (2010) formulation. The objective of this model is to minimize the total make span to finish all the required jobs. The model was built with respect to the clearance and non-crossing constraints between adjacent QCs for handling containers. There are several assumptions applied in this model as follows.

1. As in for the single ship framework, it is assumed that there is no other vessel berthing/ed during the planning horizon. The entire planning horizon is discretized, and the length of each interval is the time needed to handle the smallest unit of work, i.e. a standard 20-feet container.

2. A QC’s movement time is small compared to the time it takes to handle a container, therefore it can be ignored in the calculation of a vessel’s make span.

3. The number of allocated QCs during the planning horizon is fixed. 4. All QCs are identical and have similar work rates.

5. There is no delay in trucks delivering containers to the QCs, meaning that there is always a job to be handled in each period of time.

6. In this model, we do not differentiate between the process of loading and unloading operations, because it is expected that the processing time to loading and unloading the containers is identical.

(10)

10

Figure 3 Interpretation of variables. Source: Choo et al., 2010.

As mentioned earlier, the proposed solution algorithm will also be implemented to the new model partially built upon the quay crane scheduling model by Choo et al. (2010). To begin with, the model is developed by Boulaniki (2015) considered for a single ship. It is an integrated model of quay crane scheduling problem and yard storage location assignment, respecting to the other spatial constraint: the yard congestion constraint. This constraint limits the number of QCs handling in specific yard storage location at one time period specified with a threshold value. The main objective of the proposed model is to minimize the total travel distance of all QCs for handling containers, thus the storage location of a container is became one of the decision variables. The work of a QC is defined to be moving a container from a particular bay to a specific block in the yard area at any time period. A complete explanation of the exact model is provided in Appendix B.

5 Solution Algorithm

In this chapter, the explanation about branch-and-price algorithm will be provided as it will be used to solve the QCSP. The explanation will be followed by the improvement techniques for the B&P algorithm.

5.1 Branch-and-price Algorithm

Branch-and-price (B&P) algorithm was first introduced by Savelsbergh (1997) as an extension of dynamic column generation. In B&P algorithm, the column generation is merged with branch-and-bound (B&B) because the column generation itself does not necessarily generate all columns needed for an MIP solution. Column generation is used to solve large problems because there is an infinite number of a possible sub-problem nodes, so it is impossible to generate the master problem explicitly. In the B&P algorithm only a subset of columns is considered at each iteration. Danna & Pape (2005) stated that the B&P algorithm outperforms heuristic approaches in terms of ensuring the optimal solution. It is designed to achieve a tight bound to find the solution i.e. a lower bound in the minimization model.

(11)

11

repeating the iteration. If such columns cannot be found, the procedure is stopped. In this phase, it means an optimal dual solution has already been established to the original problem, together with the optimal primal solution to the RMP. If the solution of column generation is fractional, the branch-and-bound procedure is used to obtain integer solutions and to improve the approach hence more attractive columns can be found.

Figure 4 General branch-and-price algorithm

In order to implement B&P algorithm in solving the QCSP, the exact model needs to be reformulated as a generalized set-covering problem where each row represents a ship bay and each column serves as a QC position-bay assignment for a single time period. A column has H (number of ship bays) rows, and each having a value of ‘1’ if a QC is positioned at that particular bay. When defining an assignment, we must take into account the clearance between adjacent QCs and the ordering of QCs.

As stated before, in branch-and-price algorithm we decompose the problem into a master and pricing problem (PP). The decomposition is as follows.

Restricted master problem (RMP). In the master problem, a new variable zp for p ϵ P is

introduced as the number of times the QC position-to-bay assignment p is selected. Notation P is defined as the set of all feasible assignments, considering the clearance constraint 𝑗_𝑖+ 𝑟 ≤ 𝑗_𝑖+1 for 𝑖 = 1,2, . . , 𝐶 − 1. The illustration of variable zp can be seen in Figure 5. In the first assignment (z1),

(12)

12 [ 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1] [ 𝑧₁ 𝑧₂ 𝑧3 𝑧₄ ] ≥ [ 3 2 2 4 6 8]

The purpose of the RMP is to provide dual variables to be sent to the sub-problem and to determine if the current set of QC position-to-bay assignment ‘patterns’ provide an optimal solution or not.

(RMP) : 𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ∑𝑝∈𝑃𝑧𝑝 (1)

s.t. ∑𝑝∈𝑃:𝑗∈𝑝𝑧𝑝≥ 𝑓𝑗. 𝑗 = 1, . . , 𝐻 (2)

𝑧_𝑝 𝑛𝑜𝑛𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 𝑖𝑛𝑡𝑒𝑔𝑒𝑟 (3)

The expression of 𝑗 ∈ 𝑝 represents the fact that bay j is assigned to a QC in assignment p, in a sense there exists q such that 𝑗_𝑞 = 𝑗 if 𝑝 = (𝑗₁, 𝑗₂, 𝑗₃, . . , 𝑗_𝐶) . The objective function (1) minimizes the total number of columns (QC position-to-bay assignments) needed to meet work-load requirements which is stated in constraint (2). These constraints impose that for every bay j, we must select at least fj assignments.

Pricing problem (PP). The role of the pricing problem is to provide a column that prices out

profitably or to prove that no column exists. The optimal dual variables πj associated with constraint (2) from the RMP are imported into PP.

(PP) : 1 − 𝑚𝑎𝑥𝑖𝑚𝑖𝑧𝑒 ∑𝐻_𝑗=1𝜋_𝑗𝑥_𝑗 (4)

s.t. ∑𝐻𝑗=1𝑥𝑗 = 𝐶, (5)

𝐶(1 − 𝑥𝑗) ≥ ∑min(𝑗+𝑟,𝐻)𝑙=𝑗+1 𝑥𝑙, 𝑗 = 1, . . , 𝐻 − 1 (6)

𝐶(1 − 𝑥𝑗) ≥ ∑𝑗−1𝑙=max (1.𝑗−𝑟)𝑥𝑙, 𝑗 = 2, . . , 𝐻 (7)

𝑥 𝑏𝑖𝑛𝑎𝑟𝑦 (8)

The decision variable 𝑥𝑗 is ‘1’ if a QC is positioned at bay j and ‘0’otherwise. The objective

function (4) arises from the calculation of the reduced cost of non-basic columns. Constraint (5) states that all QCs must be positioned at a particular bay, which specify a total of C ‘active’ rows in the new assignment. Constraint (6) and (7) have similar purpose as constraint (21) – (24) in the exact model that is to ensure all QCs are assigned considering the clearance and ordering rules.

Column generation and branching scheme. Based on the RMP and PP, column generation is

conducted at every node of the branch-and-bound tree. A valid branching scheme partitions the solution space so that current fractional optimal solution of a node is excluded, integer solutions are obtained and finiteness of the algorithm is ensured (Desrosiers and Lubbecke, 2005). In this research, we branched on the fractional variables z*p with the fractional part being closest to its

(13)

13

integer value. We round up the fractional variable z*p on one branch and round it down on the other branch. In case of ties, we selected the variable with the lowest index number.

There are three available branching schemes: depth-first search, best-bound search, and combination of depth-first and best-bound (Brunner and Stolletz, 2014). Depth-first is used if the selected node has the most components fixed (i.e. one deepest in the search tree), then it is expected that integer solutions are found rapidly. On the other hand, best-bound is used if we select an active partial solution with the best parent bound (in terms of minimization is the lowest lower bound). The latter selects a deepest active partial solution after the branching node, but one with best parent bound after a termination. In this research, at the initial problem we use depth-first search to branch the fractional nodes.

Stopping criteria. The algorithm will keep iterating if the current solution is still fractional and a

new column is produced. In order to avoid large size of search tree, it is possible to discharge nodes which have optimal LP values that are greater than an integer away from the current best solution (𝑦𝑘 _{> 𝑦}

𝐵𝐸𝑆𝑇+ 1) . This rule allows us to disregard non profitable nodes, thus it will significantly

speed up the branch-and-bound procedure.

The reformulation of the integrated model for quay crane scheduling and yard storage location assignment is performed in the same manner which is presented in Appendix C.

5.2 Improvement Techniques

Despite its advantages, the B&P algorithm holds several disadvantages. Gunnerud, Foss, McKinnon & Nygreen (2012) clarified that the B&P algorithm is highly dependent on the problem being solved. Robenek et al. (2014) also stated a similar thing by saying that the convergence of this algorithm is very dependent on how well the decomposition of original problem is carried out. The method can lead to high solution times in cases where search tree becomes large, especially since column generation must be repeated in each node (Dzubur & Langvik, 2012). This phenomenon is often called as long-tail convergence or tailing-off (Alves and de Carvalho, 2006), which can be explained by primal degeneracy and superfluous oscillations of dual variables. Degeneracy happens if two or more columns have the same reduced cost value (column symmetry). This tailing-off phenomenon also occurred at every node of branching tree.

Having known these problems, researchers already tried to find ways to improve the B&P algorithm. By far, the most effective way to encounter tailing-off problem is by stabilizing the value of dual variables (Brunner and Stolletz, 2014). This technique is generally employed to improve the B&P algorithm in many different frameworks, especially when the RMP is shaped as a set-partitioning model (Grønhaug et al., 2010; Gunnerud et al., 2012; Robenek et al., 2014; Brunner and Stolletz, 2014). It is shown in those studies that using dual stabilization technique reduces the computational time for almost 90% without decreasing the effectiveness of the B&P algorithm.

(14)

14

number of runs of the associated dual LP, the gain obtained in the number of iterations in column generation does not take much time to compute better dual values. There are several dual stabilization methodologies, i.e. interior point stabilization (Rousseau, Gendreau and Feillet, 2007), boxstep method (Marsten, Hogan and Blankenship, 1975), polyhedral penalty terms (du Merle, Villeneuve, Desrosiers and Hansen, 1999), bundle methods (Briant, Lemarechal, Meurdesoif, Michel, Perrot and Vanderbeck, 2008), and convex combination with previous dual solutions (Silva, 2005). The first four methods are similar in the concept of algorithm; the difference is on the definition of stabilization center and penalty for dual values. The polyhedral penalty is considered as the most complete and flexible methodology since it takes into account both parameters to set the dual values (Lubbecke, 2010). Thus, we apply this methodology alongside with the convex combination methodology which will be explained in the following section.

Polyhedral penalty methodology. In this methodology, we defined a trust region in the dual space

to limit the value of dual variables (Figure 6). It also can be considered to be a piece-wise linear stabilization function. With the trust region, the dual value of next iteration will be restricted according to some pre-determined parameters. This methodology is considered as flexible because it allows the dual variables to lie outside the trust region, but incurs a finite penalty for deviation variables in the objective function. Using this technique, the objective function of the RMP is modified because several additional parameters and decision variables are included as follows.

Additional parameters

𝜀_𝑗−_{, 𝜀}

𝑗+ : Upper bound for deviation variables

𝛿_𝑗−_{, 𝛿}

𝑗+ : Objective function penalties for deviation variables

Additional decision variables

𝑦_𝑗−_{, 𝑦}

𝑗+ : Deviation (surplus and slack) variables for constraint (2)

SRMP (Stabilized RMP) Associated dual variables 𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ∑ 𝑧_𝑝− ∑ 𝛿_𝑗−_𝑦 𝑗∓ 𝑗 ∑ 𝛿𝑗 𝑗+𝑦𝑗+ 𝑝∈𝑃 (9) s.t. ∑ 𝑧_𝑝− 𝑦_𝑗−_{+ 𝑦} 𝑗+≥ 𝑓𝑗. ∀𝑗 𝑝∈𝑃:𝑗∈𝑝 (10) 𝜋𝑗 0 ≤ 𝑦_𝑗− _{≤ 𝜀} 𝑗−, ∀𝑗 (11) 𝜔𝑗− 0 ≤ 𝑦_𝑗+ ≤ 𝜀_𝑗+, ∀𝑗 (12) 𝑤_𝑗+

To understand the idea of stabilization, the corresponding dual program is provided as follows.

Stabilized Dual RMP (SDRMP): 𝑚𝑎𝑥𝑖𝑚𝑖𝑧𝑒 ∑ (𝑓𝑗 𝑗𝜋𝑗− 𝜀𝑗+𝜔𝑗−− 𝜀𝑗+𝜔𝑗+) (13) s.t. 𝜋_𝑗 ≤ 1, ∀𝑗 (14) −𝜋𝑗− 𝜔𝑗−≤ −𝛿𝑗−, ∀𝑗 (15) 𝜋_𝑗− 𝜔_𝑗+_{≤ 𝛿} 𝑗+, ∀𝑗 (16)

(15)

15

𝜋_𝑗 is bounded to the interval[𝛿𝑗−− 𝜔𝑗−, 𝛿𝑗++ 𝜔𝑗+], that is, deviation of 𝜋𝑗 from the soft interval

[𝛿_𝑗−, 𝛿_𝑗+] is allowed but penalized by per unit amount of 𝜀_𝑗−, 𝜀_𝑗+ respectively.

Figure 6 Penalty policy of dual value. Adapted from: Desrosiers and Lubbecke (2005)

To update the value of parameters during iterations, we use dynamic strategy as presented by Oukil, Ben Amor, Desrosiers and Gueddari (2007). The value of stabilization centers of the next iteration (t) are the optimal solution of the stabilized dual RMP (𝜋_𝑗(𝑡−1)) of the previous iteration

(t-1). The strategy to update the size of trust region and penalty is describe in the following.

1. If the current dual solution (𝜋_𝑗(𝑡)) lies inside the trust region: 𝑖𝑓 𝜋_𝑗(𝑡)𝜖]𝛿𝑗𝑖, 𝛿𝑗+[ 𝑡ℎ𝑒𝑛 {

∆_𝑗±(𝑡+1)= ∆_𝑗𝑡/2 𝜀_𝑗±(𝑡+1) = 𝜀_𝑗±(𝑡)× 2

2. If the current dual solution (𝜋_𝑗(𝑡)) lies outside the trust region: 𝑖𝑓 𝜋_𝑗(𝑡)≥ 𝛿_𝑗+_{, 𝑡ℎ𝑒𝑛 ∆}

𝑗+(𝑡+1)= ∆𝑗+𝑡 × 2 𝑎𝑛𝑑 𝜀𝑗+(𝑡+1)= 𝜀𝑗+(𝑡)/2

𝑖𝑓 𝜋_𝑗(𝑡)≤ 𝛿_𝑗−, 𝑡ℎ𝑒𝑛 ∆_𝑗−(𝑡+1)= ∆_𝑗−𝑡 × 2 𝑎𝑛𝑑 𝜀_𝑗−(𝑡+1)= 𝜀_𝑗−(𝑡)/2

With both strategies, the updated widths of ∆_𝑗± is set to a maximum 2,0 and minimum 0.1. Both values of 𝜀_𝑗± are set to maximum 1.0 and minimum 10-4_{. The stopping criteria for this stabilized} procedure are: (1) there is no other profitable column generated, and (2) 𝑦_𝑗−_{= 𝑦}

𝑗+= 0, ∀𝑗. The

(16)

16 Solve RMP Generate Pricing Problem(PP) based on stabilized dual value

Solve PP Any favourable column found?

Dual bound<Best

bound? Update the dual

bound

Solution integral Branch on fractional

variables

Add such column to RMP Solution is found YES NO YES NO YES

Update the trust region Deviation variables = 0 ?

YES

Figure 7 Stabilization scheme of polyhedral penalty

Convex combination methodology. Another stabilization technique that can be applied for

(17)

17

𝜋_𝑗𝑆𝑇(𝑡)= 𝛼𝜋_𝑗(𝑡)+ (1 − 𝛼)𝜋_𝑗𝑆𝑇(𝑡−1), 𝑡 > 1 (17)

A column is generated when it has negative reduced cost with respect to 𝜋_𝑗(𝑡). At each iteration, the dual bound is updated whenever 𝐿 (𝜋_𝑗𝑆𝑇(𝑡)) > 𝐿 (𝜋_𝑗𝑆𝑇(𝑡−1)). The stabilization scheme of this technique is depicted in Figure 8.

Initiate α value and π*=0 Solve RMP Update dual variables value Generate Pricing Problem(PP) based on stabilized dual value

Solve PP Any favourable _{column found?}

Dual bound<Best

bound? Update the dual

bound

Solution integral

Branch on fractional variables

Add such column to RMP Solution is found YES NO YES NO YES

Figure 8 Stabilization scheme of convex combination methodology

Potential enhancement in branch-and-bound phase. Another approach to reduce the

(18)

18

6 Computational Experiments

Experiments are performed in several different scenarios ranging from small problems to large problems (realistic size) by changing the input parameters. These inputs are obtained from those presented in Choo et al. (2010) and provided in Appendix D. There are two types of data sets created: small problem with 12 data instances and large problem with three data instances. We used the exact same data sets as provided from the article to ensure the validity of our improvements. We compared the result of the enhancement techniques to that obtained from research by Choo et al. (2010). Different measures are used for experimental tests, those are: (1) computational times, (2) the number of column generation iterations, (3) the number of nodes in the search tree, and (4) the optimality gap. The optimality gap is also measured in order to compare the effectiveness of the approaches. Experiments were implemented using Mosel XpressIVE 7.5 software (64-bit full version) on Intel® Core™ 2 Duo CPU 3.00 GHz personal computer with 4.00 GB of RAM and Windows 7 Enterprise as the operating system.

6.1 Parameter Settings

Several parameters that affect the performance of the computation must be set for testing. We chose the initial parameter values based on what has been done in the literature. Further, after doing trial-and-error on the initial set of testing, we can obtain the best parameter values for each problem size.

Polyhedral penalty methodology. In the first dual stabilization technique, we set the initial value

of stabilization center (𝜋𝑗′) to be zero and the initial interval for the dual value [𝛿𝑗−, 𝛿𝑗+] are set

differently to be [-1, 1], [-0.5, 0.5], or [-2, 2]. The interval of trust region can be calculated as: ∆𝑗= 𝛿𝑗+− 𝜋𝑗′= 𝜋𝑗′− 𝛿𝑗−. As an initial setting, the penalties are selected in the interval [10-4, 1].

Convex combination methodology. The only parameter that needs to be defined in this

methodology is α. We employed several tests to find the best value of α for each problem instance in a range of [0.1, 0.9] with the increment value of 0.1.

6.2 Results

Small problems. We used the computational times as the indicator to give initial insight about the

(19)

19

Table 1 Computational times of small problems

Problem Instances CPLEX (secs) Heuristic (secs) B&P from article (secs) B&P with depth first (secs) B&P depth-first & convex combination (secs) B&P depth-first & polyhedral penalty (secs) B&P with best bound (secs) B&P best-bound & convex combination (secs) B&P best-bound & polyhedral penalty (secs) SS1-1 0.87 0.09 0.15 1.2 0.799 0.1 0.252 0.08 0.1 SS1-2 0.98 0.08 0.19 1.02 0.13 0.77 0.301 0.144 0.301

SS1-3 1.1 0.08 0.14 0.76 root node root node 0.24 0.046 0.219

SS1-5 1.13 0.01 0.02 0.41 0.09 0.75 0.219 0.09 0.219

SS2-1 0.84 0.31 1.09 2.543 0.12 0.14 1.1 0.09 0.45

SS2-2 11.58 0.5 0.51 1.74 0.12 0.13 0.858 0.07 0.27

SS2-3 9.67 0.36 0.64 1.82 0.13 0.14 1.02 0.06 0.59

SS3-1 6.54 1.04 1.05 5.243 n/a n/a 4.268 0.574 2.786

SS4-2 368.7 1.82 20.66 4.26 0.1 0.141 4.16 0.07 5.19

SS5-1 n/a 2.99 105.29 7.676 n/a n./a 5.126 1.421 5.84

*n/a: unknown – takes long to compute model

The combination of depth-first search and both dual stabilization methods is not attractive in SS3-1 and SS5-1, because it needs a longer time to find the optimal solutions than the other branch-and-price schemes. Table 2 presents the output at the root node of both data instances. We can observe that the performances after applying dual stabilization are better than the ‘naïve’ B&P algorithm with depth-first search. It can be concluded that for these instances, dual stabilization techniques could improve the performance of column generation phase. On the other hand, it results to a huge number of iterations in branch-and-bound phase in order to find the optimal integer solutions.

Table 2 Results for SS3-1 and SS5-1 at the root node

Objective value Computation time (secs) Number of iterations

(20)

20

solution at the end of column generation iteration. The number of nodes in the branching tree for each algorithm is provided in Table 3.

Table 3 The number of nodes for small problems

Problem Instances B&P from the article B&P with depth first B&P depth-first & convex combination B&P depth-first & polyhedral penalty B&P with best bound B&P best-bound & convex combination B&P best-bound & polyhedral penalty SS1-1 13 24 1 1 8 1 1 SS1-2 ₉ 12 ₁ ₁ ₈ ₁ ₇ SS1-3 7 1 1 1 7 1 5 SS1-4 1 1 1 1 1 1 1 SS1-5 3 15 1 1 8 1 8 SS2-1 ₁₅ 72 ₁ ₁ ₁₆ ₁ ₉ SS2-2 1 31 1 1 15 1 1 SS2-3 4 39 1 1 20 1 13 SS3-1 1 113 n/a n/a 150 7 63 SS4-1 9 1 1 1 1 1 1 SS4-2 43 100 1 1 140 1 172 SS5-1 99 205 n/a n/a 132 13 138

As the next step, we applied the dual stabilization techniques to both of the original B&P algorithms. The implementation of dual stabilization results to fewer numbers of nodes being explored in the branching tree. When implementing convex combination approach to stabilize dual values, the optimal integer solution was directly found at the root node; meaning no branching was necessary. The polyhedral penalty method also gives the same result when implemented to the B&P algorithm with depth-first search in the branch-and-bound tree. When we implement this stabilization technique to B&P algorithm with best-bound search, it produced fractional solutions thus branching scheme was needed. Accordingly, the required nodes are fewer than the original B&P algorithm. This relates to the reduction of computational times of the algorithm.

(21)

21

Table 4 The number of column generation iterations of small problems

Problem Instances B&P with depth first B&P depth-first & convex combination B&P depth-first & polyhedral penalty B&P with best bound B&P best-bound & convex combination B&P best-bound & polyhedral penalty SS1-1 39 2 2 16 2 2 SS1-2 25 5 4 20 5 17 SS1-3 2 2 2 15 2 16 SS1-4 1 1 1 1 1 1 SS1-5 30 2 2 17 2 17 SS2-1 135 5 5 79 5 29 SS2-2 80 4 4 62 4 4 SS2-3 86 5 3 71 3 3 SS3-1 235 n/a n/a 270 40 140 SS4-1 16 16 16 270 40 15 SS4-2 231 4 3 250 3 260 SS5-1 392 n/a n/a 243 73 257

The efficiency of improvement techniques should be aligned with their effectiveness. We evaluate the effectiveness based on the value of optimal make span of each iterations. the optimality of those results.

Table 5 shows the comparisons of optimal value. It is shown that the optimal values of each algorithm that we developed are not exactly the same as those of prior research. One of the reasons of this phenomenon could be the differences of distributions of the workloads in each bay. It also affected by the type of generated columns. The optimal make span values obtained with the B&P algorithm with best-bound search are larger than those with other procedures. Furthermore, the strangest result is emerged in data instance SS5-1 when the optimal value obtained is drastically larger than those with other algorithms. The status report from XpressIVE shows that these values are not the global optimum solutions, but the optimal value of simplex dual. It means the optimal make span values obtained using branch-and-price with best-bound search are the optimal value of corresponding dual problem. Therefore, we cannot guarantee the optimality of those results.

Table 5 The optimal make span of small problems

Problem Instances B&P of Choo et al. (2010) B&P with depth first B&P depth-first & convex

(22)

22 Problem Instances B&P of Choo et al. (2010) B&P with depth first B&P depth-first & convex

combination B&P depth-first & polyhedral penalty B&P with best bound B&P best-bound & convex combination B&P best-bound & polyhedral penalty SS2-2 38 33 28 27 31 28 27 SS2-3 36 40 39 39 63 39 37 SS3-1 56 56 n/a n/a 618 48 180 SS4-1 99 101 101 101 101 106 101 SS4-2 66 80 84 79 564 84 96 SS5-1 83 98 n/a n/a 510 87 566

Large problems. In this group of problems, the complexity of problem is increased based on the

combination of the number of bays, the number of assigned QCs and the clearance distance. Although several instances in small problems have larger number of bays, but the combination of the number of QCs and clearance distance is quite simple. It is confirmed by the previous research (Choo et al., 2010), the input combination of larger problems makes them harder to be solved and it result to long computational times. It is shown that for the two largest problems, computation with exact model cannot be done since the problems are considered to be hard MIP. An attempt to solve larger instances not only aimed for the practical contribution, but also to validate the performance of each new algorithm.

The comparison of computational times of these large problems is provided in Appendix E. It is shown that our proposed algorithms reach the optimality faster than previous algorithms to a great extent. Mostly, the advanced algorithms only need less than 10s to solve the problem. Moreover, the combination of convex dual stabilization and best-bound search could solve all the instances less than 1s. On the other hand, the combination of depth-first search and dual stabilization techniques is not applicable to the last two instances. It needs significant more time than the original B&P algorithm, thus we disregarded this combination from the performance comparison to give more focus on the other algorithms.

In the Table 6 we only present the result from prior research with that from the B&P algorithm using convex dual stabilization and best-bound search since the latter is confirmed as the fastest approach to solve large problem instances.

Table 6 Performance of original B&P algorithm and the new B&P algorithm

Computational times

(secs) Number of nodes Objective value

(23)

23

Performing with the same manner as in the small problems, the B&P algorithm with convex dual stabilization and best-bound search could reach the optimality directly at the root node and they only need two iterations to find the optimal solutions. Although for SSP-2 instances the new algorithm still produce fractional solutions at the root node, the number of nodes needed still less than the original B&P algorithm. However, the optimal values still differ from the result of prior algorithm although both of them are confirmed to be the global optimum solution.

6.3 Solution Algorithm Evaluation

Based on the experiments, if we compare the performance of two branching strategies, depth-first search produced a larger search tree (more number of nodes) than best-bound search. This can happen since in the depth-first search, the active partial solution that will be chosen in the next iteration is the children of current node. It will take a longer time to find the optimal value. On the other hand, the best-bound search tries to branch the node based on the value of its parent bound. In their study, Brunner and Stolletz (2014) stated that finding the best parent bound quickly will crucially improve the efficiency of branch-and-price algorithm. We also tested the combination of depth-first and best-bound search, yet it did not give a better performance than best-bound search.

In dual stabilization approaches, we did a sensitivity analysis to find the most suitable value of each parameter for each data instances. Based on our experiments, different α values result to insignificant changes to the performance of dual stabilization. However, for SS3-1 and SS5-1 instances smaller α value gives better result than larger α in terms of total iterations required to find optimality. It means that giving a higher weight to current dual values is better than to previous one. The α value itself gives more influence to the total column generation iterations than to the number of nodes in branching tree. We present the example for SS3-1 instances in Figure 9.

Figure 9 The effect of α value in small instances

We also provide the evaluation of the adjustment for α value to the performance of the new B&P algorithm in large problems (Figure 10). As we examined in the small problems, smaller α value produced less number of columns than the larger one. Thus, we conjecture for this type of linear problem using smaller α value to stabilize the dual values is more preferable.

0,00 10,00 20,00 30,00 40,00 0 200 400 600 800 1000 1200 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 To tal tim e ( sec s) α

The effect of α value (SS3-1)

Total Iteration Number of nodes

(24)

24

Figure 10 The performance for different α value in large instances

Determining the value of parameters used in polyhedral penalty approach is more complicated since the adjustment in one parameter will affect the others. However, we already tested in the experiment setting larger size of trust region and smaller penalty value at the initial node is the best combination to solve the problem. Clearly, the parameter values will affect the performance of each enhancement. Nevertheless, we cannot come with exact conclusion about the best value of each parameter since the performance itself also related to the characteristic of the problems being solved.

6.4 Implementation to the Extended Model

As to validate the performance of the proposed solution algorithm, we implemented of the best combination of improvement techniques to the new-extended model by Boulaniki (2015). This model is more complex than the prior model by Choo et al. (2010), because it considers a new variable in the calculation. We used the same data sets as in Choo et al. (2010) to did the experiments which shown in Table 7. Additional required data for the calculation are the number of blocks in the yard, the capacity of each blocks, and the distance between bays and blocks.

Table 7 Input data for the new model

Data Sets H C R Z SS1-2 10 3 1 2 SS1-3 10 4 1 2 SS1-5 10 2 3 2 SS2-2 20 4 4 2 SSP-1 25 3 8 2

We confronted the outcomes of three different experiments: (1) the compact formulation by Mosel XpressIVE, (2) the branch-and-price algorithm with best-bound search, and (3) the advanced branch-and-price algorithm with best-bound search and convex combination of dual stabilization technique. As a realization, the time to solve the problem using the original B&P algorithm is faster than the direct exact model (Table 8). The novel B&P algorithm could find the optimal solution at the first iteration of column generation, thus it was not necessary to do the branch-and-bound procedure. This yields to the increase of the efficiency of the computation to a great extent.

0 20 40 60 80 100 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 α

The effect of α value (SSP-2)

Total Iteration

Number of nodes

(25)

25

Table 8 The output for integrated model by Boulaniki (2015)

Data Sets

Exact Model B&P with best-bound search B&P with best-bound and convex combination Optimal value Runtime (secs) Optimal value Runtime (secs) No. iterations No. of Nodes Optimal value Runtime (secs) No. iterations No. of Nodes SS1-2 3050 0,6 3770 0,281 5 2 3770 0,078 1 1 SS1-3 3050 1,4 3210 0,717 5 2 3210 0,078 1 1 SS1-5 3050 0,6 3130 0,172 3 1 3130 0,073 1 1 SS2-2 5300 14,3 6020 0,171 2 1 6020 0,124 1 1 SSP-1 25950 17,1 30670 0,546 5 1 30670 0,125 1 1

However, the optimal value of the algorithm is considerably higher. We noticed that this arises because of the differences on the definition of the optimal value. Using the B&P algorithm, the optimal value was calculated based on the total distance of all active assignments. However for each assignment there will be wastage, meaning the number of bay-block assignments is higher than needed. An assignment itself is a product of combination among several active number of bay-block pairs. As a result, further computation should be conducted in order to attain the real value of total distance (Table 9).

Table 9 The calculation of exact optimal value of the B&P algorithm

block _{No. of} containers handled No. of containers should be handled wastag e block 1 2 1 2 ba y 1 5 5 5 0 bay 1 5 2 5 5 5 0 2 5 3 7 7 7 0 3 7 4 7 7 7 0 4 7 5 7 7 7 0 5 7 6 7 7 7 0 6 7 7 7 7 7 0 7 7 8 6 6 6 0 8 6 9 5 5 10 5 5 9 5 10 5 4 9 5 4 10 5 total distance 3770 total distance 3050

(26)

26

7 Conclusion and Future Research

In this paper, we focus on improving the branch-and-price algorithm for solving the quay crane scheduling problem (QCSP) in single ship setting. We use the research of Choo et al. (2010) as the main reference in terms of developing a mathematical model and a solution algorithm to solve the problem. The objective of the model itself is to minimize the optimal make span to finish all the tasks for one singles ship. We have implemented several enhancement techniques related to the column generation procedure and branch-and-bound iteration. There are two type of techniques that have been implemented to improve the performance of the column generation procedure, namely polyhedral penalty and convex combination methodology. These two strategies have improved the performance of column generation phase in terms of the number of iterations needed to find the optimal solution. In order to improve the branching scheme, we have implemented three types of branching strategy with the aim is to reduce the size of the search tree. Those strategies are depth-first search, best-bound search, and the combination of both. Further, we also combined the enhancements for column generation with those of the branch-and-bound procedure.

The performance of the algorithms is compared with several criteria, namely: the computational times, the number of nodes in search tree, the objective value, and the number of column generation iterations. Based on the experiments, the combination of best-bound search and convex dual stabilization has the best performance of all. Additionally, we validated the performance of the leading B&P algorithm by implementing it to an integrated model built upon the same research by Choo et at. (2010). This new-developed model by Boulaniki (2015) solves the quay crane scheduling and yard storage location assignment problem concurrently. As expected, the performance of the advance branch-and-price algorithm when applied to the new model is better than the ‘naïve’ one and also the exact model formulation.

The major theoretical contribution is enhancing the branch-and-price algorithm to solve the quay crane scheduling problem. This paper fills the gap in the literature since we propose a new algorithm to solve QCSP to optimality in a fast way. In this research, it has been shown that the performance of the new branch-and-price algorithm in terms of computational times outperforms that of the heuristic approach. We also have tested that the novel branch-and-price scheme can also be implemented to a more complex MIP problem. This output gives significant impact to the practical contribution since the decision of the QCSP is critical and must be obtained in a very short time.

(27)

27

References

Alp, O.N. & Barach, H. 2009. Yard crane scheduling in container terminals – a literature review. Industrial Engineering Management, Yildiz: Turkey.

Alves, C. and de Carvalho, V. (2006). A stabilized branch-and-price-and-cut algorithm for the multiple length cutting stock problem. Computer & Operations Research (35): 1315-1328

Barnhart, C., Johnson, E. L., Nemhauser, G., Savelsbergh, M. W. P. & Vance, P. (1998), `Branch and Price: Column Generation for Solving Huge Integer Programs', Operations Research 46: 316-329.

Bierwirth, C., Meisel, F., 2009. A fast heuristic for quay crane scheduling with interference constraints. Journal of Scheduling 12, 345–360.

Bierwirth, C., Meisel, F., 2010. A survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research 202 (3), 615–627.

Boulaniki, S. (2015).Crane sequencing and storage location assignment problem in container terminals: An integrated approach. Faculty of Economics and Business University of Groningen, Groningen, NL.

Briant, O., Lemarechal, C., Meurdesoif, Ph., Michel, S., Perrot, N. & Vanderbeck, F. (2008). Comparison of bundle and classical column generation. Math. Program., Ser. A, 113:299-344.

Brunner, J.O. & Stolletz R. (2014). Stabilized branch and price with dynamic parameter updating for discontinuous tour scheduling. Computer & Operations Research (44): 137-145.

Chen J.H., Lee D.H. and Goh M. (2013). An effective mathematical formulation for the unidirectional cluster-based quay crane scheduling problem. European Journal of Operational Research 232: 198-208.

Choo, S., Klabjan, D., Simchi-Levi, D., (2010). Multiship crane sequencing with yard congestion constraints. Transportation Science 44 (1), 98–115.

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 – Methodological, 23(3), 159–175.

Danna, E. & Pape, C. L. (2005), Column Generation, Springer.

(28)

28

Du Merle, O., Villeneuve, D., Desrosiers, J. & Hansen, P. (1999). Stabilized column generation. Discrete Mathematics, 194.

Dzubur L. & Langvik A. S. (2012). Optimization of Oil Production – Applied to the Marlim Field. Norwegian University of Science and Technology.

Exposito-Izquierdo, C., Gonzalez-Velarde, J.L., Melian-Batista, B. & Moreno-Vega, M. (2013). Hybrid estimation of distribution algorithm for the quay crane scheduling problem. Applied Soft Computing, 13, 4063-4076.

Granata, D., Steeger, G. & Rebennack S. (2013). Network interdiction via a critical disruption path: branch-and-price algorithms. Computer & Operations Research (40): 2689-2702.

Grønhaug, R., Christiansen, M., Desaulniers, G. & Desrosiers, J. (2010) A Branch-and-Price Method for a Liquefied Natural Gas Inventory Routing Problem. Transportation Science 44(3):400-415. http://dx.doi.org/10.1287/trsc.1100.0317.

Gunnerud, V., Foss, B.A., McKinnon, K.I.M. & Nygreen, B. (2012). Oil production optimization solved by piecewise linearization in a Branch & Price framework. Computers & Operations Research 39, 2469-2477.

Hakam, M.H., Solvang, W.D., & Hammervoll, T. (2012) A genetic algorithm approach for quay crane scheduling with non-interference constraints at Narvik container terminal, International Journal of Logistics Research and Applications: A Leading Journal of Supply Chain Management, 15:4, 269-281, DOI: 10.1080/13675567.2012.741220.

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 Operation Research 156, 752–768.

Lajjam, A., El Merouani, M., Medouri, A. & Tabaa, Y. (2013). Mathematical moel for quay crane scheduling problem with spatial constraints. International Journal of Innovation and Applied Studies 4(3), 547-551.

Lee, D. H., Wang, H. Q., & Miao, L. X. (2008). Quay crane scheduling with noninterference constraints in port container terminals. Transportation Research Part E – Logistics and Transportation Review, 44(1), 124–135.

Legato P, Mazza R, Trunfio R. (2008). Simulation-based optimization for discharge/ loading operations at a maritime container terminal. OR Spectrum 2010;32: 543–67.

Legato, P., Trunfio, R., & Meisel, F. (2012). Modeling and solving rich quay crane scheduling problems. Computers and Operations Research, 39(9), 2063–2078.

(29)

29

Mestry, S., Damodaran, P. & Chen, C. (2011). A branch and price solution approach for order acceptance and capacity planning in make-to-order operations. European Journal of Operational Research 211: 480–495.

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(1), 45–59.

Nguyen, S., Zhang, M., Johnston, M., & Tan, K.C. (2013). Hybrid evolutionary computation methods for quay crane scheduling problems. Computer & Operations Research, 40:2083-2093.

Peterkofsky, R. I., & Daganzo, C. F. (1990). A branch and bound solution method for the crane scheduling problem. Transportation Research Part B – Methodological, 24(3), 159–172.

Robenek T., Umang N., Bierlaire M. & Ropke S. (2014). A branch-and-price algorithm to solve the integrated berth allocation and yard assignment problem in bulk ports. European Journal of Operation Research 235: 399-411.

Rousseau, L.M., Gendreau, M., & Feillet, D. (2007). Interior point stabilization for column generation. Operations Research Letters 35:660-668

Sammara, M., Cordeau, J., Laporte, G., & Monaco, M. (2007). A tabu search heuristic for the quay crane scheduling problem. Journal of Scheduling, 10(4), 327–336.

Savelsbergh, M. W. P. 2001. Branch-and-price: Integer programming with column generation. C. Floudas, P. Pardalos, eds. Encyclopedia of Optimization, Volume 1, Kluwer Academic Publishers, Boston, MA. 218–221.ilva E. F. (2004). Improving branch-and-price algorithms and applying them to stochastic programs. Naval Postgraduate School, California.

Silva, E.F. (2005). Improving branch-and-price algorithms and applying them to stochastic programs. Naval Postgraduate School, California.

Souza, F.S.H., Gendreau, M. & Mateus, G.R. (2014). Branch-and-price algorithm for the resilient multi-level hop-constrained network design. European Journal of Operational Research 233, 84-93.

Unsal, O. & Oguz, C. (2013). Constraint programming approach to quay crane scheduling problem. Transportation Research Part E, 59, 108-122.

Vacca, I., Salani, M., & Bierlaire, M. (2013) An Exact Algorithm for the Integrated Planning of Berth Allocation and Quay Crane Assignment. Transportation Science 47(2):148-161. http://dx.doi.org/10.1287/trsc.1120.0428.

Zhang, C., Liu, J., Wan, Y. W., Murty, K. G., & Linn, R. J. (2003). Storage space allocation in container terminals. Transportation Research Part B: Methodological, 37(10), 883-903.

(30)

30

Appendices

Appendix A: Model Formulation of Choo et al. (2010)

Notation:

Indices:

j bay number, in increasing order of their location on vessel, i.e. left to right; k QC number, in increasing order of their relative position, i.e. left to right; t time period index, denoting the interval of time from t-1 to t.

Parameters:

C the number of allocated WCs; H number of bays in the vessel;

fj number of containers to be handled in bay j;

T total number of time periods in the planning horizon, which can be set to ∑𝐻𝑗=1𝑓𝑗

r QC clearance value, in terms of the numbers of bays Decision variables:

𝑥_𝑗𝑘(𝑡) 1 if QC k is positioned at bay j at time period t, and 0 otherwise;

𝛿_𝑗𝑘(𝑡) 1 if QC k is handling a container at bay j at time period t, and 0 otherwise;

𝛾(t) work completion flag; 0 if all container jobs have not yet been completed at period t, and 1 otherwise;

The objective of this model is to minimize the total make span needed to finish all the jobs. The single ship model for crane sequencing problem can be written as follows.

𝑚𝑎𝑥𝑖𝑚𝑖𝑧𝑒 ∑𝑇_𝑖=1𝛾(𝑡) (18) Subject to: ∑𝐻𝑗=1𝑥𝑗𝑘(𝑡) = 1, 𝑘 = 1, . . , 𝐶, 𝑡 = 1, . . , 𝑇 (19) ∑𝐶𝑘=1𝑥𝑗𝑘(𝑡)≤ 1, 𝑗 = 1, . . , 𝐻, 𝑡 = 1, . . , 𝑇 (20) 1 − 𝑥_𝑗𝑘(𝑡) ≥ ∑𝑗−1_{𝑙=max(1,𝑗−𝑟)}𝑥_𝑙𝑚(𝑡), 𝑗 = 2, . . , 𝐻, 𝑘 = 1, . . , 𝐶, 𝑚 = 1, . . , 𝐶, 𝑡 = 1, . . , 𝑇 (21) 1 − 𝑥_𝑗𝑘(𝑡) ≥ ∑min(𝑗+𝑟,𝐻)_𝑙=j+1 𝑥_𝑙𝑚(𝑡), 𝑗 = 1, . . , 𝐻 − 1, 𝑘 = 1, . . , 𝐶, 𝑚 = 1, . . , 𝐶, 𝑡 = 1, . . , 𝑇 (22) 𝑥_𝑗𝑘(𝑡) ≤ ∑𝐻 𝑥_𝑙,𝑘+1(𝑡) 𝑙=𝑗+1 , 𝑗 = 1, . . , 𝐻 − 1, 𝑘 = 1, . . , 𝐶 − 1, 𝑡 = 1, . . , 𝑇 (23) 𝑥_𝑗𝑘(𝑡) ≤ ∑𝑗−1_𝑙=1𝑥_{𝑙,𝑘−1}(𝑡), 𝑗 = 2, . . , 𝐻, 𝑘 = 2, . . , 𝐶, 𝑡 = 1, . . , 𝑇 (24) ∑ ∑𝑇 𝛿_𝑗𝑘(𝑡) 𝑡=1 = 𝑓𝑗, 𝑗 = 1, . . , 𝐻 𝐶 𝑘=1 (25) 𝛿_𝑗𝑘(𝑡) ≤ 𝑥𝑗𝑘(𝑡), 𝑗 = 1, . . , 𝐻, 𝑘 = 1, . . , 𝐶, 𝑡 = 1, . . , 𝑇 (26) 𝛾(𝑡) =∑𝐶𝑘=1∑𝑡𝑙=1𝛿𝑗𝑘(𝑙) 𝑓𝑗 , 𝑗 = 1, . . , 𝐻, 𝑡 = 1, . . , 𝑇 (27) 𝑥, 𝛿, 𝛾 𝑏𝑖𝑛𝑎𝑟𝑦 (28)

(31)

31

are restricted from being positioned r bays to the left and right-respectively (clearance constraints). Constraint (23) and (24) describe the QC “ordering” conditions, where “higher-numbered” QC should be positioned to the right of “lower-numbered” QC (ordering or non-crossing constraints). Constraint (25) states that all required container jobs must be completed within the planning horizon. Constraint (26) ensures that a QC must be positioned at a bay if it is working there. The work completion flag defined as the number of work done for each bay. The objective function (17) ensures that when the value of the right-hand side of equation (27) is 1, ɣ(t) will be 1. Otherwise, if the right-hand side is less than 1, ɣ(t) will be zero. Constraint (28) states that all variables are binary.

Appendix B: Integrated model formulation by Boulaniki (2015)

Notation

Indices

j bay number, in increasing order of their location on the vessel and from left to right k QC number, in increasing order of their relative position and from left to right

b Yard storage location number, in increasing order of their location in the yard and from left

to right

t Time period index, denoting the interval of time from t-1 to t. Parameters

C Number of allocated QCs

H Number of bays in the vessel

𝑓_𝑗 Number of containers to be discharged in bay j 𝑑𝑗𝑘𝑏 Distance between block b and bay j if served by QC k

B Number of blocks in the yard

T Total number of time periods in the planning horizon

𝑤𝑏 Yard activity threshold, in terms of number of QCs allowed to work on containers

headed for the same storage location b at any time

r QC clearance value, in terms of the number of bays

𝑍_𝑏 Yard block capacity, in number of containers (twenty foot equivalent units)

The aim of the model is to determine the QC sequence and block storage location that minimizes the total travelling distance of QCs and therefore the discharging handling time of all containers, while respecting QC and yard congestion related constraints. The model will use the following decision variables:

Decision Variables

𝑥𝑗𝑘𝑏(𝑡) 1 if QC k allocated to block i is positioned at bay j at time period t and 0 otherwise

𝛿_𝑗𝑘𝑏(𝑡) 1 if QC k is handling a container at bay j at time period t and 0 otherwise

Both 𝑥 and 𝛿 variable are used to differentiate operating QCs from QCs that are just positioned at a bay but their status is idle. The mathematical formulation for the model is as follows:

Minimize ∑𝐻_𝑗=1 ∑𝐶_𝑘=1 ∑𝐵_𝑏=1 ∑𝑇_𝑡=1 𝛿𝑗𝑘𝑏(𝑡) ∗ 𝑑𝑗𝑘𝑏 (29)

An Improvement of Branch-and-Price Algorithm for Quay Crane Scheduling Problem in Single Ship Setting