An adaptive large neighborhood search heuristic for the multi-port continuous berth allocation problem

We present a new MIP formulation for the MCBAP. This formulation is based on the one for the continuous BAP from [25] and the one for the discrete MBAP from [42]:

Sets and parameters:

Decision variables:

Auxiliary variables:

Dependent variables:

The set of external ships $\bar{N}$ is considered fixed. Therefore, the corresponding set of decision variables $x_{i}^{c},y_{i}^{c},h_{i}^{c},r^{i,c}$ for ships $i\in\bar{N}$ are constant and given as input to the problem.

The objective function (2) minimizes the operational costs of the carriers and terminal operators. This is measured as a weighted sum of the waiting time cost, handling time cost, delay cost, and fuel consumption cost. Constraints (3) ensure that each ship berths within the available space. Constraints (4) and (5) define the relative position of each pair of ships in each dimension by enabling the auxiliary variables $\sigma^{c,c^{\prime}}$ and $\delta^{c,c^{\prime}}$. The M value can be limited to the latest finish time of the pair of ships. Constraints (6) ensure that berthing periods do not overlap in time and space. Constraints (7) compute the arrival time to a port based on the sailing speed chosen to travel from the previous port. Constraints (8) and (9) enforce that the berthing starts strictly after arrival at port and after the time window starts, respectively. Constraints (10) compute the delay if the expected finish time is exceeded and constraints (11) define if the last finish time is respected. Constraints (12) compute the handling time for each ship and port visit while constraints (13) and (14) compute the deviation from the preferred berthing position. Finally. constraints (15) ensure that only one speed is chosen to sail between ports, and constraints (16) - (19) define the domain of the decision variables.

The ALNS requires an initial solution to start with. We present a construction heuristic process for this step that aims at finding a good initial solution. Note that the BAP can be seen as a two-dimensional packing problem. However, in the continuous berth setting, the BAP has the increased complexity that the length of the rectangles (i.e., port visits) vary depending on the berthing location in the quay. In the case of the MCBAP, we are solving multiple continuous BAP problems with the additional constraint that some of those berthing times depend on a sailing time. Moreover, the fact that ships follow different routes complicates the problem as greedy approaches become harder to apply. Our construction method prioritizes reducing the delay of ships at ports. We approach this by (i) trying to place port visits early in time and close to their ideal space, therefore reducing the handling time, and (ii) by reducing ”useless” space, or, in other words, placing port visits efficiently not to create empty spots in the decision space that cannot be filled by remaining port visit. Notice that any possible solution is mathematically feasible since we allow it to exceed the latest finish time, and the time horizon is not limited. However, we aim to construct solutions where none of the ships exceed the $LFT$ as those can be perceived as infeasible by the port operators and are also heavily penalized. The method acts as a greedy heuristic, where we schedule one port visit at a time. The port visit to schedule is selected as the most constrained one. To find it, we compute the set of feasible berthing positions and times for each ship and port visit. We consider a finite set of positions and times by dividing the quay into segments of a given length (e.g., 10 meters) and the planning horizon into intervals of 1 hour. For each time instant and segment, we compute if the ship can berth starting at that time and with its left-most side starting at the segment. We do not count berthing times exceeding the latest finish time to measure how constrained a ship’s port visit is. From all unscheduled port visits, we define the one with the fewest possible positions as the most constrained one. We then schedule the port visit in one of the feasible positions. In fact, we do not consider the entire set but only the subset of feasible positions, where the port visit rectangle is directly adjacent to another scheduled port visit or to the limits of the decision space (i.e., the limit of the quay or planning horizon). From this subset of positions, we select the one resulting in the minimal change to the objective function. Besides the handling and delay cost directly computed when scheduling the port visit, we need to compute fuel consumption and waiting time costs. We consider these only if the previous port visit of the ship is scheduled.

Once a port visit is scheduled, we repeat the computation and selection of the most constrained unscheduled port visit and schedule it at the least costly efficient position. The procedure is described in Algorithm 1.

The goal of a removal operator is to select a set of scheduled port visits to be removed from the current solution. All the operators presented in this paper select $K$ number of assignments to be removed, computed as a percentage $\rho$ of the total number of port visits to be scheduled. It should be noted that removing port visits that are totally unrelated does not provide any potential gain. Therefore, a removal operator should aim at removing assignments that are related.

After applying a removal operator, the partial solution has $K$ missing port visits that need to be scheduled. They need to be assigned efficiently while respecting the other assignments and ensuring that the solution remains feasible. This is the goal of the insertion operators.

Once a new solution is reconstructed, we either accept it as the new current solution or reject it and reuse the previous one. We use a simulated annealing (SA) based criterion to take this decision. Such an acceptance criterion has been widely used for ALNS studies (see e.g.[36] and [24]). We accept the new solution $x^{\prime}$ over the current one $x$ if it is better ($f(x^{\prime})<f(x)$), or if it is worse with a probability $e^{\frac{-(f(x^{\prime})-f(x))}{T}}$, where $T$ is the current temperature at a particular iteration, and $f(x)$ is the objective function. We define an starting an ending temperature, $T_{start}$ and $T_{end}$ respectively, and the cooling time $t_{cool}$ that defines the duration of going from $T_{start}$ to $T_{end}$. Based on these parameters, we can define the cooling factor $\tau$ ($0<\tau<1$), by isolating it from the formula $T_{end}=T_{start}\tau^{t_{cool}}$. This cooling factor allows us to compute the temperature at any given instant. Given temperature $T$ at iteration $i$, we find the temperature $T^{\prime}$ to be used at iteration $i+1$ by computing $T^{\prime}=T\tau^{t_{it}}$, where $t_{it}$ is the duration of the iteration $i$. Following the strategy used in [24], we compute $T_{start}$ and $T_{end}$ based on the cost of the initial solution $f(x_{0})$ described in Section 4.1, where $\xi$ and $\phi$ define the percentage of the cost used to compute $T_{start}=\xi f(x_{0})$ and $T_{end}=\phi f(x_{0})$.

One of the main differences between the ALNS method and the standard Large Neighborhood Search (LNS) is the adaptive component of the former. The performance of the employed removal and insertion operators is measured at each iteration. These measures are then used to update the weight and, therefore, the probability of choosing the respective methods. The most common way of measuring the performance of a method is to give it a different score depending on the quality of the solution. In our case, we define four reward categories as shown in Table 1.

Let $R$ and $D$ denote the set of insertion and removal operators. Each removal and insertion method has a probability $\pi^{R}_{i},\pi^{D}_{i}$ respectively of being selected at each iteration. Throughout the algorithm run, the probability of selecting these methods gets updated depending on their performance. In our study, we update the probabilities after a $\Delta_{update}$ time interval. During these iterations we accumulate the sum of $\psi^{R}_{i},\psi^{D}_{i}$ rewards for each method, and update the weight $\omega^{R}_{i},\omega^{D}_{i}$ of each method as indicated in Equation 24

where $\lambda$ is a parameter between 0 and 1 that denotes the degree of adaptability of the method. If $\lambda=0$, the weight remains equal to the previous one. This means that each method would have the same probability throughout the entire algorithm run, behaving like an LNS with multiple neighborhoods. If $\lambda=1$, the new operator’s probability solely depends on the score achieved during the last $\Delta_{update}$ and not on previous scores. It is common to use an intermediate value for $\lambda$ strictly between 0 and 1. Once the weights are updated, the probability of each repair method $\pi_{i}^{R}$ and destroy method $\pi_{i}^{D}$ can be computed as indicated in Equation (25).

An extension of the method is implemented where we perform a local search procedure after reconstructing a new solution. This step aims to incrementally improve the solution by testing small adjustments to the port visits.

The procedure is based on the ejection chains strategy used in many routing and network-based problems (see [17, 35, 5]). The idea, in our case, is to perturbate the solution by re-planning a port visit to a better position (i.e., lower operational cost) and iteratively re-plan any port visits that conflict with the change. The chain of perturbations is limited to a maximum number of port visits to re-plan $K_{chain}$, and it terminates if this limit is reached or if a conflict-free solution is achieved. Figure 5 shows an example of this move. Note that the handling time (i.e., the vertical dimension of the port visitsships) is reduced or increased for the ships as their position changes with respect to their ideal position. A pseudo-code of the procedure is described in Algorithm 3.

The function $movePortVisit(p,n)$ performs the perturbation for a given port visit (i.e., ship $n$ at port $p$). It should be noted that the direction is given by the first perturbation made. To find the direction of the first perturbation, we compute the cost variation of moving the port visit in three directions: (i) one segment length towards the ideal position along the spatial axis, and (ii) one time instant earlier and (iii) one time instant later along the temporal axis. The direction in the spatial axis is checked if the port visit is not scheduled already at its ideal position. Once the perturbation is performed in the chosen direction, the following port visits in conflict are perturbed in the same direction.

A high value of $K_{chain}$ increases the probability of finding a better solution and the number of operations to compute. The parameter $K_{chain}$ should leverage both solution quality and low computational complexity. Therefore, we define the value of $K_{chain}$ to depend on the number of instances ships and equal to $K_{chain}=2\cdot|N|$. The reason for $K_{chain}>|N|$ is that for some movements, a conflicting port visit may require multiple perturbations to achieve a feasible new position, and selecting a lower $K_{chain}$ value may be too restrictive.

The overview of the solution method is summarized in Algorithm 4. Due to the additional computational effort of the local search procedure, we do not execute it at each iteration. Instead, we only perform it if the reconstructed solution is better than the current one. This reduces the number of times that the local search is performed, allowing the algorithm to perform more iterations while at the same time filtering the times the local search is performed to those where we already have promising solutions.

To the best of our knowledge, [31] is the only study on the MCBAP. The instances presented in the study are rather small and limited. Therefore, we develop a more comprehensive set of benchmark instances. In the absence of real-life data, one could extend current benchmarks instances of the continuous BAP to multiple ports. Instead, we decided to use the public access to port data [30], and we validated it with additional data from an industrial research partner to create an instance generator for the MBAP with a continuous quay.

We consider three different ship types: (i) feeders or small vessels with a length of up to 200 meters, (ii) medium-size vessels with a length between 200 and 300 meters, and (iii) large vessels longer than 300 meters. Each ship type has a different speed-fuel consumption relation. Moreover, we consider three terminals at the three main ports in the north sea: (i) Rotterdam APMT with a quay length of 1600 meters [3], (ii) Bremerhaven NTB, with a quay length of 1800 meters [2], and Hamburg EGH, with a quay measuring 2100 meters [14]. These three ports are relatively close to each other, and large, medium, and small vessels visit them in different sequences as part of their routes.

The duration of the vessel time window is based on the planned port call duration. The planned duration of a vessel’s port call can often be updated the days previous to the arrival time. Therefore, we establish a fixed point in time for each ship two weeks before the actual arrival time and retrieve the planned port call duration as the time difference between the estimated time of arrival (ETA) and the estimated time of departure (ETD). We compute this by averaging the planned port call duration for each port and ship type berthing in a period of three months (January-March 2021). This value is also used to define the minimum handling time $h_{0}$ (see Table 2). Port service times that exceed 48, 72, and 96 hours for small, medium, and large ship types, respectively, are categorized as outliers and removed from the dataset. The reason for this is that such long service times usually involve maintenance or fueling operations that are not usually performed on a regular basis. Thus, they are not part of the problem.

We define six different ship patterns, each with a given route, type of ship, and length. All ships visit two or three ports in different orders. The $N$ ships for a given instance are sampled from the six patterns.

For each ship, we randomize the (i) desired berthing position at each port visited and (ii) the earliest start time $EST_{i}^{c}$, following parameters ensuring that feasible sailing times between ports exist. The estimated finish time $EFT_{i}^{c}$ is computed by adding the corresponding value from Table 2 to $EST_{i}^{c}$ ($EFT_{i}^{c}=EST_{i}^{c}+h_{0}^{i,c}$). In addition, we compute the latest finish time $LFT_{i}^{c}$ by adding half of the difference between the minimum $h_{0}^{i,c}$ and maximum $h_{max}^{i,c}$ handling time that the ship can take at the port to $EFT_{i}^{c}$ ($LFT_{i}^{c}=EFT_{i}^{c}+\frac{h_{max}^{i,c}-h_{0}^{i,c}}{2}$).

As input to the instance generator, we define the number of external ships at each port $Nout$ that are considered fixed. For each external ship, we randomly define: (i) the berthing position and time, (ii) the length (comprised between 180 and 330 meters), and (iii) the handling time, adapted to be proportional to the length of the vessel.

We assume the entire quay is available for berthing unless an external ship is occupying it. We also consider 10 different speed levels, ranging uniformly between 17-21.5 knots. The ALNS heuristic we present can handle any continuous value of the speed but not the MIP formulation we present. To ensure a fair comparison of the methods we employ a discretized set of speed levels in both the formulation and the solution method. Furthermore, the distance between ports is computed based on the actual sea distance of the routes.

The ALNS algorithm has a total of 18 algorithm parameters. These include parameters used to calibrate the operators, the selection and acceptance criteria, and the weight of the scores for new solutions (see Table 4). To select the best value setting for the algorithm parameters we conducted a parameter tuning. We selected a subset of 12 instances that are representative of the entire set. For the tuning, we run the automatic algorithm configurator Pydgga [1] for the 30 generations with a time limit of 5 minutes for each ALNS run. The configurator allows to parallelize the process, and the entire parameter tuning lasted 8 hours. In Table 4, we define the domain of each algorithm parameter used for tuning and the found setting. The models and solution methods are written in Julia and run in a 2.90 GHz Intel Xeon Gold 6226R using one thread and 16 GB of RAM.

To measure the quality of the method, we compare the presented method with its variants and with a baseline commercial solvers; CPLEX v12.10.

Table 5 compares the results between CPLEX and the ALNS method. We compute the objective gap for each instance run as $\frac{z_{obj}-z_{best}}{z_{best}}$ where $z_{obj}$ is the objective value of the best solution of the run, and $z_{best}$ is the best-known solution across all experiments. We have grouped all instances per number of ships, external ships, and distance between berthing positions. Each group contains 10 instances and is named X-Y-Z according to their common characteristics. X is the number of ships, Y is the number of external ships per port, and Z is the distance between consecutive positions considered in meters (i.e., segment length). Therefore, each row in the table corresponds to the average value across instances with different seed values. The ALNS is tested by running each instance 10 times and computing the average, best and worst run.

We observe that CPLEX scales poorly, especially in instances with more than 30 ships. The gap is better for the smallest instances but quickly worsens. On average, the ALNS method outperforms the commercial solver by achieving tighter gaps in most instances. The gap is an indicator of the method performance relative to each other but does not provide an optimality guarantee. The lower bounds obtained with CPLEX indicate a high optimality gap. This could be due to a low-quality solution or a poor lower bound. [32] indicated that the relaxation of the MIP formulation for the MBAP with a discrete quay could be poor and showed that a network-flow reformulation could tighten the relaxation significantly. As indicated in [31], network-flow formulations for the MCBAP can suffer from scalability but show that the relaxation is stronger. We compare the results from the branch-and-price method presented in [31] with the ALNS method. The formulation presented in [31] is slightly different from the one addressed in this paper. The formulation from [31] defines the latest finish time for each ship berthing at a port that must be satisfied. We have adapted the method from [31] to the formulation of this study. The branch-and-price method is based on a graph representation, and therefore, we need to establish the latest possible berthing time. This is set to 50 % more than the latest finish time. This allows the method to exceed the latest finish time while maintaining the graph at a reasonable size. This is a generous bound, and our empirical studies show that this bound does not affect the optimal solution. For that, we have also generated a new set of instances of similar size to the ones presented in [31] using the instance generator defined in Section 5.1.

The input parameters for the instance set are defined in Table 6. The entire set comprises 720 instances, one for each combination of input parameters. The results are shown in Table 7, where we have grouped the instances in batches of 60 according to the number of ships. We compare the branch-and-price method with the ALNS and MIP formulation presented in this study. For both the branch-and-price and CPLEX we compute their optimality gap, where we observe that the branch-and-price method achieves a better gap due to the tighter lower bound in most cases. We also compute the gap to the best-known solution for all three methods. In this case, we observe that CPLEX provides the best performance, showing that despite its poorer relaxation, the upper bounds found are near-optimal. The branch-and-price method still shows a robust performance but for short computational times and larger instances, the ALNS method is able to provide better solutions.