A novel mathematical formulation for solving the dynamic and discrete berth allocation problem by using the Bee Colony Optimisation algorithm

Abstract

Berth allocation is one of the crucial points for efficient management of ports. This problem is complex due to all possible combinations for assigning ships to available compatible berths. This paper focuses on solving the Berth Allocation Problem (BAP) by optimising port operations using an innovative model. The problem analysed in this work deals with the Discrete and Dynamic Berth Allocation Problem (DDBAP). We propose a novel mathematical formulation expressed as a Mixed Integer Linear Programming (MILP) for solving the DDBAP. Furthermore, we adapted a metaheuristic solution approach based on the Bee Colony Optimisation (BCO) for solving large-sized combinatorial BAPs. In order to assess the solution performance and efficiency of the proposed model, we introduce a new set of instances based on real data of the Livorno port (Italy), and a comparison between the BCO algorithm and CPLEX in solving the DDBAP is performed. Additionally, the application of the proposed model to a real berth scheduling (Livorno port data) and a comparison with the Ant Colony Optimisation (ACO) metaheuristic are carried out. Results highlight the feasibility of the proposed model and the effectiveness of BCO when compared to both CPLEX and ACO, achieving computation times that ensure a real-time application of the method.

1 Introduction

During the last century, global trade and freight growth impose new challenges and requirements for efficient management of transport processes. Therefore, maritime transport is one of the crucial points of intermodal transportation, concerned with the difficulties of developing more efficient port operations. This paper focuses on solving the Berth Allocation Problem (BAP) by optimising port operations using an innovative model. In the literature, this problem has been studied by several researchers using different approaches. Basically, the BAP is a complex operations research problem based on the assignment of ships to berth areas along a quay [41]. Typically, ships arrive at a port in a specific time window, and port operators indicate their available berths. Models applied to solve the BAP can be divided according to variables related to space and time.

Spatial variables are related to the quay. According to Imai et al. [27] and Bierwirth and Maisel [1], there are three different cases of spatial variables: discrete layout (BAPD), continuous layout (BAPC), and hybrid layout (HBAP). The BAPD is the simplest and the most used model for berth planning. In this model, a quay is divided into several berths with specific lengths, and each berth can serve only one ship at a time. In the BAPC, the quay is not divided into several berths, and each ship can be served in any place without overlapping. In the HBAP, the quay is similar to the BAPD, but ships can occupy a few or more berths according to their length. Recently, Kovač et al. [32] have solved the HBAP, formulated as the Minimum Cost Hybrid Berth Allocation, by using four different metaheuristic approaches, i.e., Evolutionary Algorithms (EAs), Bee Colony Optimisation (BCO), Variable Neighborhood Search (VNS), and the General Variable Neighborhood Search (GVNS).

Temporal variables divide the BAP in two cases related to the ships’ arrival time: Static arrival (SBAP) and Dynamic arrival (DBAP) time [25]. In the SBAP, there are no expected arrival times, so it is assumed that all ships are in the harbor at the beginning of the planning horizon, waiting to be anchored. Thus, the ships’ arrival time puts a single constraint on the berthing times, with the possibility of speeding up the specific operations before the expected arrival time. In the DBAP, arrival times are scheduled and assigned to each ship, in order to avoid overlapping. Recently, Kramer et al. [34] have proposed two novel formulations: a time-indexed formulation and an arc-flow formulation for solving a big-size DBAP considering 250 ships and 20 berths instances in a reasonable computation time. Generally, various models aim at minimising the sum of the ships’ waiting time or the total staying time for completing the activities, increasing the attractiveness and efficiency of container terminal ports.

Another classification of BAP is formulated including the synchronisation between the BAP and the Quay Cranes Assignment Problem (QCAP), and it has been treated by several authors in the literature ([21, 40, 58, 60, 61]).

The problem formulated in this work deals with the Discrete and Dynamic Berth Allocation Problem (DDBAP). Generally, container terminals are highly dynamic systems, and terminal managers have to decide different port operations in short periods [56]. The right decision-making in a short time basis is crucial for port management optimisations. For this reason, an efficient optimisation algorithm that will provide an optimal berth assignment solution in real-time is a current challenge. Umang et al. [55] proposed a real-time formulation for solving the BAP with stochastic arrival and handling time on a rolling planning horizon aiming at minimising the total realised cost of the updated berthing schedule.

Heuristics and metaheuristics are the most used approaches for solving the DDBAP. Metaheuristics are general frameworks used to build heuristic algorithms for NP-hard optimisation problems in reasonable computation time. In the literature, several metaheuristic approaches were proposed, such as Genetic Algorithms [18, 28, 53], Particle Swarm Optimisation [54], Tabu Search [35] and Variable Neighborhood Search [23].

This paper proposes a novel mathematical formulation expressed as a Mixed Integer Linear Programming (MILP) for solving the DDBAP, in order to find an optimal berth assignment for a given schedule. We applied a Swarm Intelligence algorithm, the Bee Colony Optimisation (BCO), since it resulted in being very effective in solving large-sized combinatorial problems [6, 44, 52] and, to the best of our knowledge, it has not been applied to the DDBAP yet. In order to assess the performance of the model, we introduced a new set of instances based on Livorno port real data. We compared the exact solutions obtained for the proposed MILP formulation by using CPLEX with near-optimal solutions obtained by the proposed BCO algorithm. Furthermore, we applied the proposed model to a real case study (Livorno port in Italy), and we compared the results obtained by the BCO with the Ant Colony Optimisation (ACO) metaheuristic, another well-known Swarm Intelligence algorithm.

The paper is organised as follows. Section 2 presents the literature review related to the DDBAP. Section 3 introduces the proposed MILP formulation, while Section 4 describes the proposed metaheuristic based on the BCO algorithm for solving the DDBAP. Section 5 is devoted to a numerical application for comparing BCO solutions with CPLEX exact ones. Section 6 shows the application to a real case study and compares results obtained by both BCO and ACO metaheuristics. Finally, Section 7 gives concluding remarks and proposes some future developments.

2 Literature review on the DDBAP

A brief introduction of BAP has already been presented in Section 1, while in this section we aim to present a brief state-of-the-art focused on the Dynamic variant of the Discrete BAP. The most recent BAP survey paper in the literature is presented by Bierwirth and Maisel [2].

A state-of-the-art on the dynamic variants of the DBAP is necessary for highlighting the contribution of the proposed work. Since DDBAP is an NP-hard problem due to the formulation structure and to the problem size [26, 27, 48], it is necessary to apply heuristic or metaheuristic methods as a resolution approach to obtain near-optimal solutions. A summary of proposed DDBAP mathematical formulations and algorithm approaches is needed. Imai et al. [25] were the first authors who addressed the DDBAP by using a Lagrangian relaxation heuristic as a resolution approach for solving numerical experiments. Hansen and Oğuz [22] implemented a DDBAP equivalent Mixed Integer Program (MIP) correcting and extending the model proposed by Imai et al. [25]. They provided a VNS metaheuristic for solving large instances to near optimality. Thereafter, Hansen et al. [23] extended the mathematical formulation proposed by Hansen and Oğuz [22] in terms of objective function, variables, and constraints introducing the Minimum Cost Berth Allocation Problem (MCBAP) considering dynamic ships’ arrival and discrete quay division with the aim to minimise ships’ waiting and handling times. They proposed a VNS metaheuristic for solving the MCBAP medium and large size instances and demonstrated that VNS outperformed other three metaheuristics, i.e. Multi-Start, Genetic Algorithms (GAs) and Memetic search algorithm (MA). Monaco and Sammarra [45] improved the DDBAP mathematical formulation by Imai et al. [25] introducing a compact and stronger formulation dealing with the same problem. They used a Lagrangian heuristic algorithm for solving large-sized instances reaching near-optimal solutions in a reasonable computation time. Issam et al. [29] extended the Imai’s DDBAP mathematical formulation as a multi-objective MILP including ships’ draft and greenhouse gas emissions evaluation. Xu et al. [59] proposed a MILP formulation for solving the DDBAP and the SDBAP considering water depth and tidal effect with the aim to minimise the ships’ completion time. They developed two heuristic algorithms to solve the dynamic and static DBAP, respectively, and they tested the model performance on numerical experiments. Recently, Lassoued and Elloumi [38] proposed a MILP formulation for optimising the DDBAP in bulk ports. They compared the proposed MILP model with a real case study dataset by using CPLEX solver on small instances.

The most applied algorithm approaches are GA-based metaheuristics for solving large and real data instances. Nishimura et al. [47] were the first who applied GAs to address the DDBAP. Imai et al. [26] extended the formulation proposed by Imai et al. [25] including service priority constraints and applying GAs metaheuristic as the solution approach. Han et al. [20] implemented a non-linear mathematical formulation to minimise the total turnaround time and to improve the terminal operation efficiency by using GAs and a hybrid optimisation strategy combining GA and Simulated Annealing (GASA), respectively. Zhou et al. [62] proposed a non-linear programming model considering ships’ total waiting and service times as stochastic parameters. They used a GAs and designed a greedy algorithm for improving ships’ assignment to available berths. Golias et al. [14] proposed a MILP formulation with two objective functions for minimising ships’ late departure costs and for maximising ships’ early and timely departure benefits within a time window. Imai et al. [28] proposed a two-objective mathematical formulation considering ships’ service quality and berths’ utilisation objective functions by minimising delay in ships’ departure and total service time, respectively. They applied GAs and Subgradient Optimisation (SP) approaches for solving the DDBAP. Theofanis et al. [53] implemented a linearised version of the DDBAP formulation proposed in Imai et al. [26], and they compared two heuristics, i.e. GAs and the Optimisation-Based GA (OBGA), on a small dataset. Golias et al. [15] proposed a multi-objective combinatorial optimisation problem considering ships’ priority for minimising the total service time. They developed a GA-based metaheuristic to solve large and real data instances. Golias et al. [19] proposed an additional non-linear MIP formulation for minimising ships’ emissions, waiting time, and delayed departures by using a GA-based metaheuristic to solve real data instances. Saharidis et al. [49] formulated a bi-level optimisation problem introducing a hierarchy among the different objective functions. They used an improved interactive GA-based metaheuristic, i.e. the k-th best algorithm, for solving real data instances. Golias et al. [16] proposed a lamda-optimal based heuristic algorithm, namely the Defined Neighborhood heuristic, for solving berth scheduling sub-problems minimising the total weighted service time. They also presented a GAs metaheuristic as an alternative solver for the cases in which no-efficient objective optimiser is available to solve sub-problems. Golias and Haralambides [17] proposed a non-linear DDBAP mathematical formulation for minimising ships’ late departure total costs and waiting time, and for maximising ships’ early departure total premiums simultaneously. They assumed different variable penalty/premium costs functions based on contractual agreements between the liner shipping company and the terminal operator. Additionally, they used a GA-based metaheuristic algorithm to solve computational examples based on four different berth scheduling policies. Golias et al. [18] proposed a bi-level bi-objective optimisation problem for minimising the average and range of total service times depending on ships’ arrival and handling time uncertainties. They developed a GA-based metaheuristic algorithm for solving big-sized numerical experiments. Simrin and Diabat [51] proposed a non-linear MIP for solving the DDBAP with the aim to minimise the ships’ total handling and waiting times. Moreover, they linearised the previous model and developed a GA metaheuristic as a resolution approach for solving small and large set of instances.

In the literature, other authors dealt with the DDBAP by using different metaheuristic approaches. Cordeau et al. [4] proposed two mathematical formulations considering both discrete and continuous cases with the aim to minimise the ships’ total weighted service time. They developed a Tabu Search (TS) metaheuristic for solving discrete and continuous cases for big-sized instances while they applied CPLEX solver for small-sized instances. Buhrkal et al. [3] improved the mathematical formulation by Cordeau et al. [4] by introducing a Generalized Set Partitioning Problem (GSPP) mathematical formulation. Lalla-Ruiz et al. [34] proposed a hybrid heuristic combining a TS metaheuristic with an artificial intelligence method, namely the Path Relinking algorithm, for solving the DDBAP. They compared exact solutions of the GSPP mathematical formulation proposed by Buhrkal et al. [2] and near-optimal solutions of the TS metaheuristic proposed by Cordeau et al. [4]. Furthermore, de Oliveira et al. [4] proposed an alternative Clustering Search method using the Simulated Annealing (SA) metaheuristic applied on the mathematical formulation by Cordeau et al. [4]. Lin and Ting [42] proposed two versions of the SA algorithm, i.e., with and without restarting strategy, for solving the DDBAP by using the mathematical formulation by Buhrkal et al. [3]. They tested both SA algorithms on several small and medium instances in the literature. Lin et al. [43] proposed the Iterated Greedy (IG) heuristic for solving the DDBAP minimising ships’ total service time. Lee and Jin [39] dealt with the DDBAP considering the optimisation between feeder ships schedule and terminal operations in a container transhipment terminal. They proposed a mixed integer programming model and a MA approach, i.e. the combination between GAs and TS.

Karafa et al. [31] proposed a modified version of the DDBAP model proposed by Golias et al. [16] considering ships’ handling time as a stochastic parameter. They formulated a bi-objective optimisation problem with the aim to minimise the berth schedule risk and to maximise the berth throughput by using an Evolutionary EA-based heuristic and a simulation-based Pareto front pruning algorithm. Dulebenets [9] proposed a non-linear MIP model for solving the DDBAP with the aim to minimise the ships’ total weighted service time. He developed a modified version of MA with a Deterministic Parameter Control (MA-DPC) algorithm and compared it with other four metaheuristics, i.e., MA, EA, SA, and VNS algorithms by testing it on numerical instances. Furthermore, Dulebenets et al. [13] proposed a MILP model for solving the DDBAP minimising the ships’ total service cost, including the total carbon dioxide emission cost. They developed a Hybrid Evolutionary Algorithm (HEA) for solving a set of numerical instances. Dulebenets [10] proposed a MILP for solving the DDBAP including service priority with the aim to minimise the ships’ total weighted service cost. He developed an Adaptive Evolutionary (AEA) metaheuristic with a parameter control strategy to test the model effectiveness through extensive numerical experiments. Dulebenets et al. [12] proposed a MILP model aimed at minimising the ships’ total weighted turnaround time and the ships’ total weighted late departure. They improved the previous AEA algorithm proposed by Dulebenets [12] with a Self-Adaptive Evolutionary Algorithm (SAEA). Additionally, Dulebenets [11] proposed a MILP model aiming at minimising the ships’ total weighted service cost. He developed an Adaptive Island Evolutionary Algorithm (AIEA) which can execute simultaneously multiple EAs in parallel on its islands based on an adaptive mechanism. Ting et al. [54] proposed a MIP model for solving the DDBAP as a vehicle routing type problem with the aim to minimise the ships’ total service time. They developed a Particle Swarm Optimisation (PSO) metaheuristic algorithm to test and validate their model through two different size sets of benchmark instances in the literature. Lalla-Ruiz and Voß [36] and Lalla-Ruiz and Voß [37] proposed a Partial Optimisation Metaheuristic Under Special Intensification Conditions (POPMUSIC) metaheuristic for solving the DDBAP applying the mathematical formulation proposed by Cordeau et al. [4]. They proposed two variants of POPMUSIC metaheuristic, i.e. POPMUSIC and POPMUSIC-G, applying different initial solution methods, i.e. Random-Greedy method and First-Come First-Served Greedy, respectively. Recently, Nishi et al. [46] proposed a Dynamic Programming-Based Matheuristic (DP-Math) based on POPMUSIC metaheuristic to solve the DDBAP proposed by Cordeau et al. [4]. They compared DP-Math metaheuristic with GSPP and dynasearch algorithms on large scale instances. Issam et al. [30] proposed a DDBAP mathematical formulation with the aim to minimise the ships’ total staying time at a container terminal. They developed a Modified Sailfish Optimiser (MSFO) metaheuristic and compared MSFO with other heuristic-based approaches in the literature on a set of instances proposed by Cordeau et al. [4]. Wang et al. [57] proposed a DDBAP mathematical formulation minimising the ships’ total service time considering tidal effects. They developed a Levy-Flight (LF) metaheuristic to solve the proposed model and compared LF with other state-of-the-art exact and heuristic methods on large scale instances.

In Table 1, we report a summary of the state-of-the-art as an extension of the survey presented by Kovač [32], in which main differences in models and metaheuristic algorithms proposed for solving the DDBAP are highlighted.

Table 1 Problems and metaheuristic algorithms for the DDBAP

According to Table 1, the proposed mathematical formulation can be classified as discrete (disc), and draft ships dependent (draft) as a spatial attribute; dynamic (dyn), and due time-dependent (due) as a temporal attribute; and position-dependent handling time (pos) as a handling time attribute. Among the above-mentioned authors, Zhou et al. [62], Han et al. [20], Xu et al. [59], and Issam et al. [29] have considered draft constraints in their mathematical formulations. Accordingly, as highlighted by Sheikholeslami et al. [50], we considered the compatibility between ships’ draft and berths’ draft since it is relevant, especially for shallow ports where tide effects and low-draft access (e.g., channel port access) could influence the berth planning negatively.

Thus, considering the existing literature, the contribution of this paper is summarised as follows:

• 1. A novel mathematical formulation expressed as a MILP for solving the DDBAP was proposed. Additionally, we introduced a new parameter that considers the berth approaching/leaving time, to the best of our knowledge, not considered in the literature yet.

• 2. A metaheuristic for solving the DDBAP, the BCO, was proposed since there are few Swarm Intelligence algorithms applied for solving the DDBAP in the literature. Furthermore, a new set of benchmark instances, based on Livorno port real data, for the DDBAP was introduced to test the feasibility of the proposed optimisation model by using an exact solver and the BCO.

• 3. An application to a real case study (Livorno port) considering a long-time horizon (from 1 up to 4 weeks) was introduced. We compared the real monthly berth planning with the solution obtained by the proposed BCO and another well-known Swarm Intelligence method, the ACO, to assess the effectiveness of the proposed approach.

3 Problem description and formulation

In this paper, we propose a novel mathematical formulation for solving the DDBAP, formulated as a MILP. The formulation was obtained starting from the Non-Linear Mixed Integer Programming (NLMIP) model reported in Appendix 1 in order to provide an equivalent linear model for algorithms comparison. In both formulations, the objective function minimises the total ships’ staying time. Our model considers the following assumptions: the port, at its initial state, is empty, and each ship corresponds at least to one compatible berth due to the discrete distribution model of quays. Furthermore, compatibility between ships and berths is related to geometric constraints, such as length and draft. Ships arrive at the port in the harbor area according to scheduled arrival times and arranged in chronological order. Arrival, service, and total ships’ staying times are related to a fixed time reference, while berth approaching/leaving, handling, and waiting times are fixed time intervals. Berth approaching/leaving, and handling times are calculated as the average of observed values, respectively. The warehouse capacity for the storage of goods and/or containers is unlimited. The number and performance of cranes are not considered. Finally, ships can move freely within the port: priority and any additional waiting times due to internal congestion have been excluded. The list of notations used for the MILP formulation is reported in Table 2 including notations used in the NLMIP (Appendix A).

Table 2 List of notations used in problem formulations

All definitions, as depicted in Fig. 1, are given as follows:

• Arrival time: the arrival time of a ship to a given port in the harbor area.

• Waiting time: the time spent by a ship, assigned to a berth, in waiting for the previous ship to leave the same berth.

• Berth approaching/leaving time: the total time necessary for a ship to approach and leave its assigned berth. This time includes all mooring manoeuvres starting and finishing at the harbor area.

• Handling time: the time spent by port operators for loading or unloading goods of a ship, assigned to a berth.

• Service time: the sum of waiting, berth approaching/leaving, and handling times of a ship assigned to a berth.

• Leaving time: the port leaving time is the sum of arrival and service time of a ship assigned to a berth.

• Total staying time: the cumulative sum of arrival and service times of all ships assigned to their respective berths.

Fig.1

A space-time representation of the proposed DDBAP model

3.1 The MILP mathematical formulation

In order to find the global optimal solution of the DDBAP using an exact solver, we modelled the problem as a MILP whose mathematical formulation is given as follows:

$$ {\displaystyle \begin{array}{c}\min \sum \limits_{i=1}^{n_s}{t}_i^t\\ {}\mathrm{s}.\mathrm{t}.\end{array}} $$

(1)

$$ \sum \limits_{j=1}^{n_b}{x}_{ij}=1,\forall j\in B $$

(2)

$$ {s}_i={t}_i^w+{t}_i^a,\forall i\in S $$

(3)

$$ {s}_i\ge {t}_i^a,\forall i\in S $$

(4)

$$ {s}_i+M\cdot \left(1-{x}_{ij}\right)\ge {s}_k+{t}_k^h+{t}_{kj}^b-M\cdot \left(1-{x}_{kj}\right),\forall i\in S,\forall k\in S,\forall j\in B,k>i $$

(5)

$$ {l}_i^s\cdot {x}_{ij}\le {l}_j^b,\forall i\in S,\forall j\in B $$

(6)

$$ {d}_i^s\cdot {x}_{ij}\le {d}_j^b,\forall i\in S,\forall j\in B $$

(7)

$$ {t}_i^B\ge {t}_{ij}^b\cdot {x}_{ij},\forall i\in S,\forall j\in B $$

(8)

$$ {t}_i^t={t}_i^a+{t}_i^h+{t}_i^B+{t}_i^w,\forall i\in S $$

(9)

$$ {x}_{ij}\in \left\{0,1\right\},\forall i\in S,\forall j\in B $$

(10)

The proposed MILP formulation (1)–(10) aims at minimising the total staying time by the objective function (1). Constraints (2) ensure the uniqueness of the assignment. We define the leaving time s_i in constraints (3), while constraints (4) and (5) ensure that s_i must be consistent with all previous arrival times. We introduce a large positive number M in constraints (5) in order to avoid unfeasibility and to faster reach the optimal solution. Constraints (6) and (7) ensure the compatibility between the length and the draft of the ships with respect to those of the berths. Constraints (8) evaluate the berth approaching/leaving time auxiliary variables ${t}_i^B$related to the assignment of the i-th ship to the j-th berth, allowed by the right-side term ${t}_{ij}^b\cdot {x}_{ij}$. Finally, constraints (9) and (10) define the total staying time decision variables ${t}_i^t$ and the binary nature of the decision variables x_ij, respectively.

4 The Bee Colony Optimisation approach

In order to face the complexity of the problem, we applied the BCO metaheuristic proposed by Teodorović et al. [52] to solve the DDBAP.

Generally, the BCO is based on bees’ natural collaborative behaviour for solving the complex combinatorial optimisation problem. We used, in this case, artificial bees for performing the search process. Firstly, bees are located in the hive and then they start to interact when the search process is undertaken. Each artificial bee undertakes a set of moves through stages determining a partial solution of the problem. These partial solutions remain until some other suitable ones have been found. At each stage, bees’ movements are called forward and backward step. During the forward step, they generate partial solutions through previous individual or collective experiences, while during the backward step, they return to the hive for participating in the decision-making process. At the end of the stages, we completed an iteration generating one or more feasible solutions. The search process continues until a fixed maximum number of iterations is reached. When applied to the DDBAP, each bee (b) represents a solution, i.e., a feasible assignment of ships to available berths. We can represent an assignment matrix using a three-dimensional structure, as shown in Fig. 2a.

Fig. 2

a Three-dimensional assignment matrix structure; b flowchart of the proposed method

The complete list of all parameters used for BCO is reported in Table 3.

Table 3 List of Notations used in BCO algorithm

In the initialisation, we define the maximum number of iterations (n_k) and the number of bees (n_bee). Each solution is constructed through stages following a given time schedule according to the dynamic nature of the problem. Thus, the number of stages (G) is equal to the number of ships (n_s) in the schedule, and a stage represents a single ship-to-berth assignment. As stated before, each stage is composed of a forward and backward step. During a forward step, each bee finds partial assignments of the ships to the berths. All the partial solutions are identified observing the constraints of the optimisation problem, and the associated fitness value is given by the objective function (1). During a backward step, each bee could improve its partial solution or move to a better one. At the end of each iteration, a bee constructs a feasible solution. Then, a new iteration starts searching for new solutions until the maximum number of iterations is reached. The flowchart of the proposed method is shown in Fig. 2b. In Algorithm 1, we reported the proposed BCO pseudo-code to better understand the forward/backward process related to each iteration k of the algorithm.

Algorithm 1. Pseudo-code applied to solve the DDBAP

As a result, a path of the created artificial network corresponds to a berth assignment found by a bee. At the end of each iteration, all the possible solutions found are evaluated, referring to the associated fitness value, and the best assignment is saved. Then, a new iteration starts searching for new solutions until the maximum number of iterations n_k is reached [44]. We can represent the decision space as an artificial network composed of nodes and stages, as shown in Fig. 3a and b.

Fig. 3

a Forward step and b Backward step

The roulette wheel selection [24], applied to the BCO algorithm, is a common algorithm used to select an item proportional to its probability. The main inspiration for its development derives from the example of roulette used in gambling. Each solution can be chosen, and the probability of selection (the size of a roulette slot) depends on its fitness value, which is one value of the objective function. In practice, the size of the circular sector of the roulette associated with each solution is determined by the ratio between the normalised value of the corresponding objective function and the sum of the normalised values of the objective function for all the solutions found. On the one hand, a solution with a better value than the objective function is more likely to be selected; on the other, there is a probability that it will be eliminated following a subsequent search process. The probability that a partial or complete solution of b can be chosen by any unengaged bee is equal to the following formulation:

$$ {p}_b=\frac{O_b}{\sum \limits_{r=1}^{n_r}{O}_r},\forall b\in R $$

(11)

The probability that the b-th bee, at the beginning of the next forward step, remains loyal to its own total or partial solution found in previous steps is represented as follows:

$$ {p}_b^{u+1}={e}^{-\left(\frac{O_{\mathrm{max}}-{O}_b}{u}\right)},\forall b\in B,\forall u\in G $$

(12)

The normalisation of O_b, related to the minimisation of the objective function (1), is calculated by the following equation:

$$ {O}_b=\frac{Q_{\mathrm{max}}-{Q}_b}{Q_{\mathrm{max}}-{Q}_{\mathrm{min}}},\forall b\in B $$

(13)

The normalised value O_b represents the fitness value of a partial or complete solution. According to Eq. (11)–(13) and a random number generator, each artificial bee decides whether to become an engaged follower or to continue exploring its own solution. If the random number chosen is less than the calculated probability, the bee remains faithful to its solution. Otherwise, if the random number is greater than the probability, the bee becomes unengaged. A higher O_b value corresponds to a better-generated solution and to a greater probability that the bee remains faithful to the solution discovered previously. The higher the index value of the step u the higher the probability that the bee remains faithful to its previously found solution, according to Eq. (12).

5 Numerical application

In order to check the feasibility of the proposed optimisation model, we compared the BCO with an exact solver for the DDBAP as formulated above. We tested both algorithms on a PC equipped with an Intel i7 2.4GHz CPU and 4 Gb of RAM. The BCO algorithm was developed in C++, while CPLEX was used as an exact solver for the MILP. We introduced a new set of instances based on a Livorno port scheduling real dataset considering 25, 50, 75, 100 ships and 5, 10, 15, 20 berths. Generated instances are available at the following URL: https://bit.ly/37avUC7. The proposed benchmark instances are obtained as follows: the arrival time is considered as a fixed value for all instances, while all the other parameters are assigned randomly from the real dataset. Moreover, according to problem constraints, at least one berth must be compatible with all ships in terms of length and draft.

Concerning the BCO algorithm, we set n_k = 1000, n_bee = 10, n_c = 1, n_r = n_bee, while computation times (t) are related to the average CPU time over 10 runs. Regarding CPLEX, we set the time limit for finding the global optimum equal to 7200 s. The results comparison is reported in Table 4.

Table 4 Results comparison between CPLEX and BCO

We can observe that for small-sized instances, the exact solution is achieved by both approaches. Furthermore, for medium and large-sized instances, we can observe that results in terms of the objective function gap give low values (0.41% in the worst case). In contrast, the comparison in terms of CPU time is relevant. The proposed BCO approach assures near-optimal results in low computation time, thus, it could be used for real-time applications in port management operations. Accordingly, the proposed model could avoid data uncertainties by updating the berth scheduling process in real-time.

6 Application to a real case and results

We applied the proposed method to the real case of the Livorno port by monitoring and acquiring data referred to the container terminals of Darsena Toscana (5 berths) and Sintermar (3 berths), as shown in Fig. 4. We used a time schedule of 28 days referred to the berth planning of February 2016 and considered a total of 140 container ships. Acquired data included arrival times, staying times, ships’ length and draft, berths’ destination, and status referred to a daily berthing planning and extended to a monthly time schedule. In order to better evaluate the effectiveness of the proposed method, we compared the obtained results with the real schedule of the Livorno port and the Ant Colony Optimisation (ACO) metaheuristic. The comparisons were carried out executing 10 runs of the BCO algorithm. The results of the comparison with the Livorno port assignment are reported in Table 5 in terms of objective function value, service time, and waiting time. We can observe that the percentage of improvement in terms of the total service time is 9% over 4 weeks. This represents a significant result, considering the wide schedule time interval observed. Moreover, we obtained a saving in terms of total waiting time of about 10 cumulated days when compared to the real schedule.

Fig. 4

Container terminal quays considered in the Livorno port (Italy)

Table 5 Comparison between observed and computed times considering 1000 iterations and 100 bees over 10 runs – monthly berth plan

The optimisation of the total staying time improves the performance of the container terminal and, consequently, of the whole port for better ships’ management. The lower the objective function, the greater the number of ships that the port can serve for the same time interval, thus increasing the volume of the handled goods.

Figure 5a shows a comparison between the observed and computed total staying times. We can observe significant improvements in terms of service time (Fig. 5b) and waiting time (Fig. 6a). The latter resulted to be about 80% lower than the observed times. In Fig. 6b, we can highlight a resulting homogeneous distribution of the ships over the berths.

Fig. 5

Comparison between observed and computed (a) ships’ total staying time (obj. function), (b) service time – monthly berth plan

Fig. 6

Comparison between observed and computed (a) ships’ waiting time and (b) the percentage of berths’ usage – monthly berth plan

In the second comparison, we adapted the Ant Colony Optimisation (ACO) to the DDBAP by developing it in C++ as the proposed BCO. The ACO is one of the most used and known Swarm Intelligence method metaheuristics, firstly proposed by Dorigo [7]. It is an algorithm based on the behaviour of ants seeking a path between their colony and a source of food, suitable to search for an optimal path in a graph. Further details can be found in Dorigo and Gambardella [8]. We applied the ACO to the Livorno port schedule and executed 10 runs as the previous comparison. The results of the comparison are reported in Table 6. We can observe that the proposed BCO found better results than ACO in terms of minimum and average objective function value expressed in hours. In particular, we obtained a saving of 3.5 h in terms of minimum value and about 3 h on average. Standard deviation values are almost the same, which means that both metaheuristics are characterised by high robustness.

Table 6 Comparison between ACO and BCO considering 1000 iterations and 100 ants/bees over 10 runs – monthly berth plan

Finally, the computation performance of the BCO in terms of CPU time (seconds) are shown in Table 7 and Fig. 7. We can assert that the obtained computation times, on average 4 min in the worst case, are compatible with real-time applications for container ports.

Table 7 Comparison of the BCO average CPU time (s) over 10 runs varying the number of iterations and bees

Fig. 7

Computational time values varying the number of iterations and artificial bees

7 Conclusions

In this paper, we proposed a novel mathematical formulation to solve the Discrete and Dynamic Berth Allocation Problem (DDBAP). Starting from a Non-Linear Mixed Integer Programming (NLMIP) we formulated the DDBAP as a Mixed Integer Linear Programming (MILP) to find the global optimum with an exact solver. The objective function minimises the total ships’ staying time, introducing a new parameter related to the berth approaching/leaving time. Moreover, a solution method based on the Bee Colony Optimisation (BCO) was introduced for its capabilities in solving complex problems. In order to check the feasibility of the proposed optimisation model, we introduced a new set of instances based on Livorno port (Italy) real data. We compared the optimal solutions of the proposed MILP, obtained using the exact solver CPLEX, with near-optimal solutions obtained by the proposed BCO. Results highlighted the effectiveness of the BCO in terms of near-optimal solutions (average gap of 0.12%) in a reasonable low computation time (on average 10.3 s). Furthermore, we applied the proposed model to the real case of the Livorno port in Italy over a time horizon of 4 weeks (February 2016). We compared the proposed BCO with another well-known Swarm Intelligence metaheuristic, the Ant Colony Optimisation. Results highlighted the effectiveness of the proposed approach in terms of average improvement obtained over 10 runs. We evaluated improvements in terms of objective function value, service, and waiting times. As a result, we obtained a saving of 10 cumulative days when compared to the real berth assignment of the Livorno port. Moreover, BCO resulted in being more effective than ACO, obtaining an improvement of about 3 h on average. Finally, computation performance highlighted the possibility of a real-time application of the proposed model to support port operators in the berth assignment task.

From a practical point of view, several benefits can be achieved. It is possible to serve a greater number of ships in a certain time interval by minimising the total staying time. Thus, the port performance increases proportionally in terms of volume of goods, economic benefits due to anchorage fees, port taxes, possible freight surcharges, cost reduction of port operators, and resource management. Other benefits can also be referred to shipping companies. The reduction of waiting time can result in different cost benefits, such as a reduction in personnel costs, possible anchored surcharges, maintenance, and fuel costs. Furthermore, non-monetizable costs such as water noise and air pollution can be also reduced.

Further developments will concern the integration of the proposed DDBAP model with the Quay Crane Assignment Problem (QCAP) for container terminals and to include data uncertainty, budget limit, and priority, in the mathematical formulation. To this aim, the NLMIP formulation, reported in Appendix 1, can serve as a basis to develop extended non-linear models that embed other optimisation problems in port management like the QCAP. Additionally, the comparison with other metaheuristic approaches would be considered in future research.

Appendix 1. The NLMIP mathematical formulation

According to the notations reported in Table 2 and the assumptions presented in Section 3.1, the NLMIP formulation aiming at minimising the total staying time is given as follows:

$$ {\displaystyle \begin{array}{l}\min \sum \limits_{i=1}^{n_s}\sum \limits_{j=1}^{n_b}{t}_{ij}^t\cdot {x}_{ij}\\ {}\mathrm{s}.\mathrm{t}.\end{array}} $$

(14)

$$ \sum \limits_{j=1}^{n_b}{x}_{ij}\le 1,\forall i\in S $$

(15)

$$ \sum \limits_{j=1}^{n_b}{l}_i^s\cdot {x}_{ij}\le {l}_j^b,\forall i\in S $$

(16)

$$ \sum \limits_{j=1}^{n_b}{d}_i^s\cdot {x}_{ij}\le {d}_j^b,\forall i\in S $$

(17)

$$ \sum \limits_{j=1}^{n_b}{l}_i^s\cdot {x}_{ij}+{d}_i^s\cdot {x}_{ij}\le {l}_j^b\cdot {x}_{ij}+{d}_j^b,\forall i\in S $$

(18)

$$ {t}_{ij}^W=\max \left\{0,-{t}_i^a+\sum \limits_{k=i-1}^1{t}_{kj}^t\cdot {x}_{kj}\right\},\forall i\in S,\forall j\in B,{t}_0^a=0,{t}_0^h=0,{t}_{0j}^b=0 $$

(19)

$$ {t}_{ij}^s={t}_{ij}^W+{t}_{ij}^b+{t}_i^h $$

(20)

$$ {t}_{ij}^t={t}_i^a+{t}_{ij}^s $$

(21)

$$ {t}_{ij}^w\ge {t}_{ij}^W,\forall i\in S,\forall j\in B $$

(22)

$$ {t}_{ij}^w\le M\cdot \left(1-\left|{x}_{ij}-{x}_{\left(i-1\right)j}\right|\right),\forall i\in S,\forall j\in B $$

(23)

$$ {x}_{0j}=0,\forall j\in B $$

(24)

$$ {x}_{ij}\in \left\{0,1\right\},\forall i\in S,\forall j\in B $$

(25)

The objective function (14) minimises the total staying time. Constraints (15) ensure the uniqueness of the assignment between the i-th ship and the j-th berth, i.e., a ship can only be assigned to a single berth. Constraints (16)–(18) are related to the compatibility between the length and the draft of the ships with respect to those of the berths. The waiting times ${t}_{ij}^W$ are calculated by constraints (19) and they are calculated as the maximum value between zero and the difference between the total staying time of the (i-1)-th ship and the arrival time of the i-th ship assigned to the same j-th berth. Constraints (20) and (21) define the service time ${t}_{ij}^s$ and the total staying time, respectively. Constraints (22) and (23) define the assignment of the waiting time auxiliary variable introducing a large positive number M. Accordingly, the waiting time auxiliary variables ${t}_{ij}^w$ are the minimum value between the waiting times ${t}_{ij}^W$and M. Consequently, constraints (24) and (25) introduce three possible cases, referred to the value of the decision variables x_ij:

• 1 x_ij = 1, x_{(i − 1)j} = 1: ${t}_{ij}^w$is equal to the waiting time ${t}_{ij}^W$ according to the constraints (22);

• 2 x_ij = 0, x_{(i − 1)j} = 0: ${t}_{ij}^w$is equal to zero according to constraints (19);

• 3 x_ij = 1, x_{(i − 1)j} = 0, or vice versa: ${t}_{ij}^w$ is equal to zero, according to the constraints (23).

Constraints (24) ensure that the decision variables in constraints (19) and (23) with the subscript i-1 are equal to zero when i = 1. Finally, constraints (25) define the binary nature of the decision variables.