Minimizing the Total Service Time of Discrete Dynamic Berth Allocation Problem by an Iterated Greedy Heuristic

Berth allocation is the forefront operation performed when ships arrive at a port and is a critical task in container port optimization. Minimizing the time ships spend at berths constitutes an important objective of berth allocation problems. This study focuses on the discrete dynamic berth allocation problem (discrete DBAP), which aims to minimize total service time, and proposes an iterated greedy (IG) algorithm to solve it. The proposed IG algorithm is tested on three benchmark problem sets. Experimental results show that the proposed IG algorithm can obtain optimal solutions for all test instances of the first and second problem sets and outperforms the best-known solutions for 35 out of 90 test instances of the third problem set.


Introduction
Containerization has been widely adopted in global freight transportation since the 1950s. Containerization significantly reduces shipping costs and accelerates cargo handling at ports. According to UNCTAD [1], maritime transportation is an important component of the global supply chain, with more than 8.7 billion tons of goods shipped annually. Shippers and carriers benefit from their ships spending as little time harbored in port as possible. Therefore, terminal authorities strive to provide efficient and cost-effective services that maximize terminal efficiency.
The operations of a container terminal include seaside operations, yard operations, and land-side operations [2,3]. One important issue in seaside operations is the assignment of berthing position to a defined set of ships that must be served within a defined planning horizon. This is esteemed as a key process for container terminals. Due to operational correlation, container terminal operators and ocean carrier share the common objective in minimizing the time ships spend in port. The berth allocation problem (BAP) thus arises in dealing with terminal assignments of berths to ships, and one of its major aims is to minimize the total service time, that is, waiting time plus handling time. Various types of aims/objectives pursued by the BAP exist and can be referred to as in the excellent reviews of the literature [3] for the corresponding objective functions.
The BAP can be categorized using temporal and spatial constraints [3]. In terms of temporal constraints, BAP can be static or dynamic, while in terms of spatial constraints BAP can be for discrete, continuous, or hybrid berthing spaces. The static berth allocation problem (SBAP) disregards ship arrival time. That is, ships arrive before berth allocation is planned. The dynamic berth allocation problem (DBAP) assumes ships can arrive at any time with future arrival information being known; ships cannot be berthed before their arrival time. It should be mentioned that the DBAP solved assumes all data (including arrival times) are known in advance and therefore no reoptimization is required.
In the discrete BAP, the quay is divided into a set of berths, each of which can harbor only one ship at a time. In the continuous BAP, the quay is not partitioned into definitive 2 The Scientific World Journal berths, and a vessel can occupy any arbitrary position on the quay. This improves utilization of quay space at the cost of greater computational complexity. In the hybrid BAP, the quay is partitioned into berths, but large ships may need multiple berths, while small ships require only one berth.
The DBAP is known to be NP-hard [4] and thus is generally solved using metaheuristics. The iterated greedy (IG) algorithm is a very effective and efficient metaheuristic algorithm [5]. The major advantages of the IG algorithm are its simplicity and its extension property to be applicable to different problems. IG has exhibited state-of-the-art performance for numerous problems [6][7][8]. Therefore, this study proposes an IG algorithm to solve the discrete DBAP.
The remainder of the paper is organized as follows. Section 2 describes the discrete DBAP and reviews the relevant literature. Section 3 presents the proposed IG algorithm for tackling the discrete DBAP. Section 4 depicts computational experiments and discusses the results. Finally, Section 5 presents concluding remarks.

Literature Review
Several mixed integer programming (MIP) models for discrete DBAP have been proposed in the literature. Imai et al. [9] were the first to present a discrete DBAP model which was an extension of an SBAP. Decision variables were used to assign ships to berths and scheduled their processing sequence in each assigned berth. Cordeau et al. [4] formulated the problem as a multidepot vehicle routing problem with time windows (MDVRPTW). Their model included different positions along the berth that resulted in different handling times for each ship.
Christensen and Holst [10] modeled the discrete DBAP as a generalized set-partitioning problem (GSPP) which assumed that the time measurements were integers (discrete time periods). In their model a column represented a feasible assignment of a single ship to a specific berth at a specific time. Buhrkal et al. [11] proposed a heterogeneous vehicle routing problem with time window (HVRPTW) model which was a simplified version of MDVRPTW. The HVRPTW model defined the problem on a graph rather than on a multigraph as shown in Cordeau et al. [4] and left the complexity of the problem unchanged. Furthermore, Buhrkal et al. [11] proposed an improved HVRPTW model, denoted as HVRPTW+, by considering break symmetry, variable fixing, and adding valid inequalities to reduce the computation time of HVRPTW. Based on various existing models, the GSPP model optimized solutions of the ship assignment problem [11]. However, the GSPP was still NP-hard. The number of columns became excessive when a large number of berthing time and ship were involved [11]. Furthermore, Vacca et al. [12] developed an exact algorithm for the integrated planning of berth allocation and quay crane assignment. Their exact algorithm could solve the BAP + QCAP (berth allocation problem and quay crane assignment problem) to optimality on the same set of instances of [4] using branch-andprice. Because of the complex nature of the problem, global optimal solutions may be difficult to obtain when the problem involves a large amount of variables. Therefore, researchers have been seeking efficient approximation algorithms that obtain near-optimal solutions reasonably fast.
Lai and Shih [13] developed a heuristic algorithm for solving the DBAP by considering the first-come-first-served rule and evaluated three different berthing policies using simulation experiments. Brown et al. [14,15] explored berth allocation models that allow multiple ships to occupy a single berthing position. Imai et al. [16] formulated an SBAP as a nonlinear integer programming model to minimize the weighted objective of two conflicting criteria: berth performance and berthing satisfaction. Imai et al. [9] introduced a dynamic version of the problem, where each ship had a given arrival time and could only receive service after arrival. The objective was to minimize the accumulative service time of all ships.
Nishimura et al. [17] extended the DBAP to the multiwater depth configuration in a public berth system and proposed a genetic algorithm (GA) to solve the problem. Imai et al. [18] extended the DBAP to assign service priorities to ships and developed a GA-based algorithm to solve the extended problem. Cordeau et al. [4] considered a DBAP with time windows and proposed a Tabu search algorithm to solve it. Monaco and Sammarra [19] derived a more compact formulation than that of Imai et al. [9] and solved it by using a Lagrangian relaxation algorithm and a nonstandard multiplier adjustment method. Imai et al. [20] formulated a BAP that allowed a single berth to simultaneously serve two ships and solved it using GAs. Imai et al. [21] used the Lagrangian relaxation with subgradient optimization and the GA to identify noninferior solutions in a biobjective BAP model which minimized the service and delay time.
Mauri et al. [22] proposed a population training algorithm with linear programming (PTA/LP) to solve discrete DBAP. PTA/LP improved incoming columns in the column generation problem. Hansen et al. [23] presented a minimum cost berth allocation problem (MCBAP) based on an extension of the model previously proposed by Imai et al. [18] and developed a variable neighborhood search algorithm to solve the problem. Imai et al. [24] studied a variant of the DBAP in which an external terminal could be available when the port ran out of berth capacity.
Barros et al. [25] formulated a berth allocation model with tidal time windows, where berths could only be operated when the tide allowed. Their simulation model was based on the dynamic berth allocation model of Imai et al. [9]. De Oliveira et al. [26] proposed an approach based on clustering search (CS) with simulated annealing mechanism. CS was an iterative method which divided the search space into clusters and comprised a metaheuristic for solution generation, a grouping process, and a local search algorithm. Lalla-Ruiz et al. [27] proposed a Tabu search (T 2 S * ) approach and a Tabu search with path relinking (T 2 S * + PR) approach to solve the discrete DBAP. T 2 S * was an improved version of T 2 S [4] that employed different neighborhood structure, and T 2 S * + PR added the path-relink techniques to the T 2 S * . T 2 S, T 2 S * , and T 2 S * + PR were tested using the instances from Cordeau et al. [4] and the newly generated problem set. The experimental The Scientific World Journal results showed that the T 2 S * + PR approach was competitive with GSPP in small-scale problems and outperformed T 2 S and T 2 S * . Lin and Ting [28] proposed two versions of simulated annealing (SA) algorithm to solve the discrete DBAP. The results showed that both versions obtained the same solutions as those of GSPP and outperformed the T 2 S, PTA/LP, and CS approaches. Furthermore, the version of SA algorithm with restarting strategy (SA RS ) outperformed that without restarting strategy (SA WRS ).

Development of the Proposed Iterated Greedy Heuristic
A generic IG algorithm usually starts from an initial solution 0 and then generates a sequence of solutions by iterating the greedy method through the destruction and construction phases [5]. The destruction phase obtains a partial candidate solution by removing a fixed number ( ) of elements from the current candidate solution . In the subsequent construction phase, a greedy constructive approach is used to sequentially insert the removed elements into the partial solution ( ) until a full solution is reconstructed. Once a complete solution is reconstructed, an acceptance criterion is applied to determine whether the new solution should replace the existing one. The process iterates through the destruction and construction phases until reaching the termination conditions. Additionally, another local search method may be applied before both the main loop and acceptance test to improve the initial solution and the reconstructed solution.
Based on the framework of the generic IG algorithm, the following subsection further discusses the solution representation, the objective function calculation, and the main steps of the proposed IG algorithm.

Solution Representation and Objective Function Calculation.
A solution can be represented by a numerical sequence that consists of a permutation of ships and − 1 zeros, where denotes the number of berths. That is, the numerical sequence contains segments separated by "zeros, " where each segment corresponds to the service sequence of certain ships on an assigned berth. Figure 1 represents an example of a solution, which is explained as follows: 15 ships are to be processed on three berths, and the service sequences of the ships on berths 1, 2, and 3 are 3-8-10-9-7, 2-4-15-12-11-5, and 13-14-6-1, respectively.
The completion time of each ship on an assigned berth is calculated according to its arrival time, the sequence in the berth, and the availability of the berth. The service time of each ship is obtained by subtracting its arrival time from its completion time. Finally, the total service time can be calculated by summing up the service times of all ships.

Main
Steps of the Proposed IG Algorithm. The main steps of the proposed IG algorithm are as follows.
Step 1 (generate the initial solution). The initial solution is generated using the first-come-first-served rule, which is usually implemented in real world operations. That is, the ships are sorted in ascending order of their arrival times. Each ship is then sequentially assigned to the berth with the earliest completion time of the ship that has been assigned to it. In the event of a tie, the berth with the shortest waiting time of the assigned ship is chosen. If the tie persists, the berth with the smallest number will be selected. The obtained initial solution is set as the incumbent solution * and the best solution * best .
Step 2 (destruction and construction phases). Consider the following.
(a) Randomly choose distinct ships from the ships of * . The value of is selected randomly between min and max , where min and max are the minimal and maximal numbers of unrepeated ships to be removed, respectively. Subtract these numbers from * and add them to * in the order in which they are chosen, where * is a permutation list of the removed ships.
(b) Sequentially reinsert the ships of * into * until a complete solution * new is obtained, where * is the partial sequence of * obtained after removing ships. When inserting a ship into * , all possible positions in all berths of the incumbent partial solution are considered. The best position is then chosen and recorded as the incumbent partial solution.
, and is the temperature.
Step 3 (stopping criteria). If the computational time exceeds a specified threshold, stop the algorithm. In Step 1, an initial solution is generated according to the first-come-first-served rule. Steps 2(a) and 2(b) are the destruction and construction phases, which comprise a perturbation mechanism. In Step 2(c), the Boltzmann function that is commonly used in the annealing process of SA algorithms is applied to enable the proposed IG algorithm to escape the local minimum. This is achieved by generating a random number ∈ [0, 1] and replacing the incumbent , the probability of replacing * with * new is set to one. Subsequently, the incumbent solution * can be improved by the destruction and construction phases of Step 2 iteratively until the computational time reaching a predetermined limit. To clearly illustrate the process, Tables 1 and 2 together give a small discrete DBAP instance with 15-ship and 3-berth. The start time of berth ( ) and finish time of berths ( ) are listed in Table 1. The arrival time of ship ( ), the end of the service time window on ship ( ), and the handling time of ship at each berth ( 1 , 2 , 3 ) are given in Table 2. Figure 2 presents an iteration of the proposed IG algorithm.
The computational complexity of the proposed IG algorithm is as follows. In Step 1, the time complexity needed for sorting the ships in ascending order of their arrival times is ( log 2 ). When trying to assign each ship to the berth, there are at most berths to be considered. In each iteration of Step 2, distinct ships are removed from the ships of current solution in the destruction phase. The computational complexity is linear. In addition, in each iteration of Step 2, ships are needed to be reinserted to the partial solution in the construction phase. When trying to find the best position to reinsert the first one of removed ships, there are ( + − ) possible positions to be tested. Therefore, there are ( + − ) times to calculate objective function value. In a similar fashion, when trying to find the best position to reinsert the th ship of removed ships, there are ( + − 1 − + ) possible positions to be tested. Therefore, the total number of calculating objective function is ∑ =1 ( + − 1 − + ). It should be noted that recalculation of the objective function is localized. That is, only the ships whose orders are affected by     1  71  300  20  20  40  2  90  300  44  44  88  3  39  300  22  22  44  4  17  300  34  34  68  5  12  300  12  12  24  6  1 1 7  3 0 0  3 0  3 0  6 0  7  9 4  3 0 0  2 8  2 8  5 6  8  2 9  3 0 0  6  6  1 2  9  43  300  26  26  52  10  79  300  22  22  44  11  2  300  20  20  40  12  129  300  16  16  32  13  123  300  26  26  52  14  43  300  14  14  28  15  5  300  18  18  36 the inserted ship on the same berth are taken into account for recalculating the objective function. The proposed IG algorithm is thus adaptive and efficient.

Computational Results and Discussion
This section discusses the computational tests used to evaluate the performance of the proposed IG algorithm. The details The Scientific World Journal 5   [4] provided two sets (I2 and I3) of instances that were randomly generated based on data from the port of Gioia Tauro (Italy). The I2 set includes five instance sizes: 25 ships with 5, 7, and 10 berths; 35 ships with 7 and 10 berths; and 10 instances generated for each size. The I3 set includes 30 instances with 60 ships and 13 berths. Lalla-Ruiz et al. [27] provided a new set of instances that was generated according to Cordeau et al. [4] with longer time horizon, higher traffic, and fewer available berths. New instances can be found at https://sites.google.com/site/ gciports/berth-allocation-problem. This site provided nine instance sizes: 30 ships with 3 and 5 berths; 40 ships with 5, 7, and 10 berths; 55 ships with 5, 7, and 10 berths; and 60 ships with 5 and 7 berths. Ten instances are generated for each problem size.

Parameter Selection.
The proposed IG algorithm was implemented using the C language on the Windows XP operating system and run on a personal computer with an Intel Core 2 2.66 GHz CPU and 2 G RAM. Parameter selection may influence the quality of the results. One instance was randomly selected from each size in the I2 problem set and the new problem set, and three instances were randomly selected in the I3 problem set for preliminary testing. The following combinations of parameters were tested on these instances:

Results and Discussion
. Tables 3-6 list the computational results for the discrete DBAP. The optimal solution was provided by the GSPP model using CPLEX 11 [10]. CPLEX was able to get optimal solutions for all small-sized problem instances of the three benchmark problem sets. However, for part of medium-sized and all large-sized problem instances, CPLEX is terminated with memory depletion and no solution was found [27]. Tables 3, 4, and 5 list the required computation times for the GSPP model. The proposed IG algorithm is compared with SA RS for the I2 problem set and results are listed in Table 3. Table 4 lists the results of the PTA/LP, CS, SA RS , and IG algorithms for the I3 problem set. Furthermore, Table 5 compares the proposed IG algorithm with T 2 S * , T 2 S * + PR, and SA RS for the new problem set with known optimal solutions, while Table 6 lists the results of T 2 S * , T 2 S * + PR, SA RS , and IG for the new problem set with unknown optimal solutions. In Table 3, column one represents the name of the instance, while columns two and three are the optimal solutions and the computational time required by the GSPP, respectively. Furthermore, columns four to five display the best solutions obtained by SA RS in 10 trials and the computational times required for SA RS to obtain the optimal solutions, respectively. Columns six and seven show similar information for IG. As shown in Table 3, both the IG and SA RS can obtain the optimal solutions for all instances of the I2 problem set. Table 4 lists the computational results for the discrete case of the I3 problem set. Besides the above columns in Table 3, Table 4 lists information for PTA/LP and CS. The CS, SA RS , and IG obtain all optimal solutions, whereas PTA/LP cannot reach the optimal solution for all instances. In such case, PTA/LP is only 1 time unit away from optimality. Notably, if the optimal solution cannot be obtained within a certain number of trials, the computational time is de facto the maximum computational time (Max = 0.2 × seconds) for these trials. This table indicates that the proposed IG algorithm is as effective as CS and SA RS in optimally solving small-scale problems.
In Table 5, columns one and two represent the problem size and instance number. Columns three and four show the optimal solutions and the computational time using the GSPP. Furthermore, columns five to seven display the best solutions obtained by T 2 S * in 30 trials, the required time, and The Scientific World Journal 7    1.20 0.00 the gaps relative to the optimal solutions for T 2 S * , respectively. The gap is calculated as (Sol ℎ − Sol GSPP )/Sol GSPP × 100%, where Sol ℎ is the solution obtained by algorithm ℎ and Sol GSPP is the optimal solution obtained by GPSS. Columns 8-10 listed similar information (best solutions obtained in 30 trails, the required time, and the gaps relative to the optimal solution) for T 2 S * + PR. Columns 11 to 14 show the best solutions, the average solutions, the relative gaps between best solutions to the optimal solutions, the time required to obtain the optimal solutions, and the maximal computation time for SA RS in 30 trials. Similar information for the proposed IG algorithm is listed in Columns 15 to 18. The table shows that 20, 27, 30, and 30 out of 30 optimal solutions are obtained by T 2 S * , T 2 S * + PR, SA RS , and IG, respectively. The SA RS and the proposed IG algorithm perform best for new problem sets with known optimal solutions. Table 6 lists similar information to Table 5, except that the optimal solution is unknown and replaced by the BKS, which exhibits the best solution among the T 2 S * , T 2 S * + PR, SA RS , and IG approaches. For 60 problems, all solutions obtained by the proposed IG algorithm are equal to BKS, while the 13, 28, and 47 solutions obtained by T 2 S * , T 2 S * + PR, and SA RS heuristic are equal to BKS.
The running time of the IG algorithm depends on various factors, including CPU, operating system, compiler, computer program, and the precision used during execution. Therefore, the relative efficiency of the algorithms is hard to determine. Tables 3-6 show the time required for GSPP, PTA/LP, CS, T 2 S * , T 2 S * + PR, SA RS , and IG. The GSPP formulation was implemented by generating all columns a priori using a Java program and solving the resulting integer program using CPLEX 11 (32-bit version) on a PC with an Intel Xeon 5430 (2.66 GHz) processor. T 2 S was implemented in ANSI C, and computational experiments were performed on a Sun workstation (900 MHz). Both the PTA/LP and CS were implemented using C++ and run on a PC with an AMD Athlon 64 3500 with a 2.2 GHz processor and 1 GB of RAM. T 2 S * and T 2 S * + PR were coded in Ansi C and the experiments were performed on a PC with a T4300 processor at 2.10 GHz. The computational time and results for GSPP, PTA/LP, CS, T 2 S * , and T 2 S * + PR are cited from their original papers, while the SA RS is recalculated by setting the maximal computational time to be the same as that of the proposed IG algorithm. As shown in Tables 3-6, the computational time is within the acceptable range for the proposed IG algorithm.
To verify the effectiveness of the proposed IG algorithm, the proposed algorithm is compared with PTA/LP, CS, SA RS , T 2 S * , T 2 S * + PR, and IG by conducting paired -tests on the RER (relative error rate). The RER is calculated as (TST ℎ − BKS * )/BKS * × 100, where TST ℎ denotes the total service time from algorithm ℎ. The BKS * is obtained from all the compared algorithms, including those of Cordeau et al. [4], Mauri et al. [22], Buhrkal et al. [11], de Oliveira et al. [26], Lin and Ting [28], and the proposed IG approach.  At a confidence level of 95%, Tables 7 and 8 show that the proposed IG algorithm outperforms the PTA/LP in terms of the best objective value obtained for discrete DBAP (PTA/LP for the I3 set). However, IG is not statistically better than SA RS and CS on the best objective value obtained, probably because these state-of-the-art algorithms are equally able to obtain the best solutions as IG. For the new data set, the proposed IG algorithm outperforms the T 2 S * , T 2 S * + PR, and SA RS approaches, as shown in Table 9. This outperformance demonstrates the superiority of the proposed IG algorithm.

Conclusion
This paper studies the berth allocation problem with dynamic arrival time. Because the berth allocation problem is NPhard, exact solution approaches cannot optimally solve realistic large-scale problems while maintaining acceptable computational complexity. An IG algorithm is proposed as an alternative method to the problem. The proposed IG algorithm is tested using three benchmark problem sets and compared with the optimal solutions (or best known solutions) from the literature. Computational results indicate that the proposed IG algorithm is effective. The proposed IG algorithm obtains all the optimal solutions of the discrete DBAP instances for the first and the second problem sets, as well as exhibiting best-known solutions for 35 out of 90 test instances in the third problem set. Future research can further examine the integration of the berth allocation and quay crane assignment problems.