Port optimization: the integrated berth allocation and crane assignment problem

Universitat Polit`ecnica de Catalunya Facultat de Matem`atiques i Estad´ıstica Degree in Mathematics Bachelor’s Degree Thesis Integrated Berth Allocation and Crane Assignment Problem Joan Salv`a Soler Supervised by: Andr´es Puy - Data Science Consultant at Accenture January, 2021 Thanks to Andr´es for his excellent guidance, his great ideas, and the learnings I take from him both professionally and personally. Thanks to Sandra, for the opportunity to work on this thesis as an Innovation project at Accenture, and for her bright ideas and support. Thanks to Jordi for the provided knowledge and expertise in Operations Research. Thanks to my family, friends, and roomies, for the support given. Abstract Port logistics play a key role in today’s global economy, and improving its efficiency has been a subject of research for the past two decades. The Berth Allocation and Crane Assignment problem consists in deciding the location of incoming vessels in a port and defining a crane schedule to serve them. We consider a version of the problem that conceives distinguishable cranes and multiple materials. In this research we give a MILP formulation of the problem and design a Genetic Algorithm that we use for speeding up an exact solver (CPLEX) and for solving large instances that are intractable in an exact way. Keywords Berth allocation, Crane assignment, Port optimization, Genetic Algorithm, Mixed-Integer Programming 1 Contents 1 Introduction 3 1.1 Real-world problem ...................................... 3 1.2 Mathematical problems to solve ............................... 3 1.3 Motivation for this Bachelor’s Thesis ............................. 4 1.4 Process to optimize in the Port of Bayside .......................... 4 1.5 The problem we will focus on in this Bachelor’s Thesis ................... 4 2 The Berth Allocation and Crane Assignment Problem 6 2.1 Versions of the problem .................................... 6 2.2 Literature review ........................................ 6 3 Mathematical formulation 8 3.1 Assumptions and definitions .................................. 8 3.2 The space-time diagram .................................... 9 3.3 Formulation .......................................... 10 3.4 Formulation analysis ...................................... 14 4 Genetic algorithm 15 4.1 Genetic algorithms overview .................................. 15 4.2 Design and implementation .................................. 16 4.2.1 Encoding of the solutions ............................... 16 4.2.2 Crossover ....................................... 17 4.2.3 Mutation ........................................ 18 4.2.4 Selection ........................................ 20 4.2.5 Evaluation ....................................... 21 4.3 Convergence .......................................... 21 4.4 Paremeter tuning ....................................... 22 5 Results 27 5.1 Exact solving with CPLEX .................................. 27 5.2 Mixed Integer Program starts ................................. 28 5.3 First Come First Serve solution ................................ 28 5.4 Solution of different problem instances ............................ 29 5.5 Optimizing the Port of Bayside ................................ 31 5.6 Conclusions .......................................... 34 5.7 More ideas and further research ................................ 34 2 1. Introduction 1.1 Real-world problem Maritime transportation plays a crucial role in the international trade and global economy. Around 80% of the volume of international trade in goods is carried by sea. According to the UN Conference on Trade and Development [12], the cargo volume of maritime transportation has gone through a fast growth. Specifically, the carried cargo of oil tankers, bulk carriers, container ships, and other types of ships have increased an average of 4.3%, 9.8%, 27%, and 11.7% per year respectively, as shown in Figure 1.1. However, this rapid development of the world transported tonnage comes with important challenges on the capacities of ports, leading to frequent vessel delays and heavy port congestion. A recent example is the Port of Los Angeles: arriving vessels were forced to spend up to three weeks before they were allowed to dock and unload, this having massive implications in the global supply chain [16]. It is thus an object of study to improve the operational efficiency of ports. Figure 1: Millions of tonnage carried per year, disaggregated per vessel type 1.2 Mathematical problems to solve Different mathematical problems arise when trying to improve the operational efficiency of ports. A possible solution is to address the port berth allocation problem (BAP) of how to assign arriving vessels to a specific location in a given berth layout so as to minimize the waiting times and thus achieve high productivity. The physical design of some ports sometimes complicates the BAP problem. It’s very common that vessels need to go through a one or two-way channel in order to access the berth area, defining a bottleneck (BAP with Channel Constraints). By looking at the problem from the point of view of the cranes, one gets the Crane Assignment Problem (CAP). It’s the problem of assigning jobs (to unload or load parts of the vessel) to mobile cranes, in order to minimize either the total operations time or the total traveled distance of the cranes. There’s also room for optimization in the container terminal on the land. Containers can be stored in a way that minimizes the total movements of containers, the movements of cranes, or the cost of work of the hoppers that will later transport them to the next locations. 3 Port Optimization 1.3 Motivation for this Bachelor’s Thesis The consulting firm Accenture worked on a logistics optimization project for the company that managed a trading port in Europe some years ago. From now on we will refer to this port as the Port of Bayside, which is a fictional name. Two consulting firms were involved: Accenture was in charge of the simulation part of the project that had to provide insights for the client, and the other team was working on defining a planning optimization problem and solving it. Accenture has an interest now in reproducing the optimization solution that the other company developed, in hopes that the developed asset will be useful in future port optimization projects. 1.4 Process to optimize in the Port of Bayside The process that was wanted to be optimized is the following: cargo vessels arrive at the outside of the port and send a notice of readiness to the port. When there is enough space in the port and the operators are ready, the vessel berths and the unloading and loading operations start. Vessels usually carry either containers or bulk, which are handled by different types of cranes. Some cranes are fixed to the ground and some are mobile. At the time they are being unloaded, bulk materials are loaded onto hoppers (trucks) that bring them to the corresponding warehouses. Containers are stored in a container terminal. These hoppers are also in charge of bringing the cargo that has to be loaded to vessels from the port warehouses to next to the cranes. Summarizing, the main ’bottlenecks’ or constraints of this process are the following: there must be enough space in the berthing area for a vessel to moor, cranes must be available in order to start working on a vessel, cranes must be compatible with the cargo that they are working on, hoppers have to be available in order to process and distribute the cargo that was just unloaded. 1.5 The problem we will focus on in this Bachelor’s Thesis Intuition told us that a mathematical formulation of the global port logistics problem would be too complicated (no literature is found on the BAP + CAP + {Cargo on land management problem}). Considering future applications of the project, we decided to focus on the Integrated Berth Allocation and Crane Assignment Problem (BACAP problem). If we retake the process that we described in the section above, BACAP assumes that the fleet of hoppers has ’infinite’ productivity when it comes to handle and store the materials that the cranes are unloading, and to bring the materials to load next to the cranes that will do it. We are not concerned about how the cargo is handled once unloaded or how it was handled before being loaded. Accenture has some data on the Port of Bayside operations, and more detailed information on how the logistics of the port worked was available through the consultants that worked on that project. But rather than focusing on the particularities of this port, we wanted to define a more general mathematical optimization problem that could be applied to ports with different physical and logistic characteristics. However, one thing we want our model to be is a generalization of Bayside, so all particularities of this port must be included. What were these particularities of Bayside that we want to include in the model? As we previously explained, two different kinds of vessels were served on the port. Some of them carried bulk materials that could only be handled using a special type of cranes, while others carried containers that had to be handled 4 by another set of cranes. Therefore, our model will consider that vessels carry different types of materials, and for each crane and material, we will be given a different productivity for loading and unloading. This productivity could be zero in the case that one crane can’t handle a specific material. In the Port of Bayside, some vessels just unload cargo and leave the port, whereas other unload and load back. Our model will treat separately the two processes, so the amount of work on a vessel will be given disaggregated in to-unload and to-load. The port had some of its cranes fixed to the ground, so if a vessel hold has to be operated by that crane, it must berth right in front of it. Our model will allow cranes to be either fixed or mobile. Mobile cranes in Bayside had wheels (didn’t move on rails), so we will allow cranes to cross each other if they have to work in opposite directions on the berth. 5 Port Optimization Constraint (9.1) says that the departure time of the vessel has to be after all holds have been loaded, and after some departure preparation time (for the holds that had to be loaded). ck≥˜ δk+X t,q zk,it,q·(t+pk,iq) (9.2) for i∈Uk−Lk, 1 ≤k≤n Constraint (9.2) says that the departure time of the vessel has to be after the unload of the holds that don’t need to be loaded back. X t,q zt,q k,i= 1 for i∈Uk, 1 ≤k≤n(10.1) X t,q ˜zk,it,q= 1 for i∈Lk, 1 ≤k≤n(10.2) Constraints (10.1) and (10.2) say that only (and at least) one crane will unload a hold, and only (and at least) one crane will load it back. X t,q zt,q k,i= 0 for i/∈Uk, 1 ≤k≤n(11.1) X t,q ˜zk,it,q= 0 for i/∈Lk, 1 ≤k≤n(11.2) Constraints (11.1) and (11.2) say that no crane will unload a hold that comes empty or load a hold that must be left empty. X k,i zt,q k,i+ ˜zk,it,q≤1 for 1 ≤q≤Q, 1 ≤t≤T(12) Constraint (12) says hat a crane can only be assigned to one place at a time. M·(1 − T X s=t+1 zs,q k,i) + T X s=t+1 s·zs,q k,i≥M·( t X s=1 zs,q l,j−1) + t X s=1 zs,q l,j·s+pq l,j+vq· |bk−bl+i−j|(13.1) for i∈Uk, 1 ≤k≤n,j∈Ul, 1 ≤l≤n for 1 ≤t≤T, 1 ≤q≤Q M·(1− T X s=t+1 ˜zk,is,q)+ T X s=t+1 s·˜zk,is,q≥M·( t X s=1 zs,q l,j−1)+ t X s=1 zs,q l,j·s+pq l,j+vq·|bk−bl+i−j|(13.2) for i∈Lk, 1 ≤k≤n,j∈Ul, 1 ≤l≤n for 1 ≤t≤T, 1 ≤q≤Q 12 M·(1− T X s=t+1 zs,q k,i)+ T X s=t+1 s·zs,q k,i≥M·( t X s=1 ˜zl,js,q−1)+ t X s=1 ˜zl,js,q·s+ ˜pl,jq+vq·|bk−bl+i−j|(13.3) for i∈Uk, 1 ≤k≤n,j∈Ll, 1 ≤l≤n for 1 ≤t≤T, 1 ≤q≤Q M·(1− T X s=t+1 ˜zk,is,q)+ T X s=t+1 s·˜zk,is,q≥M·( t X s=1 ˜zl,js,q−1)+ t X s=1 ˜zl,js,q·s+ ˜pl,jq+vq·|bk−bl+i−j|(13.4) for i∈Lk, 1 ≤k≤n,j∈Ll, 1 ≤l≤n for 1 ≤t≤T, 1 ≤q≤Q Constraints (13.1), (13.2) ,(13.3), (13.4) say that a crane can start working on other places after finishing its previous job. It considers the 4 possible cases (from unload to unload, unload to load...). (bk+i−1) ≤bw+M·(1 −zt,w k,i−˜zk,it,q) (14.1) (bk+i−1) ≥bw−M·(1 −zt,w k,i−˜zk,it,q) (14.1) for w∈G, 1 ≤t≤T, 1 ≤i≤hk, 1 ≤k≤n Constraints (14.1), (14.1) say that if a crane w∈Gis fixed at location bw, then all operations it will do will be at location bq. B·uw k,i·M+ (bk+i−1−bw)·M+X t zt,q k,i≥1 (15.1) for w∈G, 1 ≤k≤n,i∈Uk B·(1 −uw k,i)·M−(bk+i−1−bw)·M+X t ˜zk,it,q≥1 (15.2) for w∈G, 1 ≤k≤n,i∈Lk Constraints (15.1), (15.2) say that if a hold is located in a berth with a fixed crane, then that crane is going to unload and load it. We use the additional variable uto model the fact that bk+i−1−bwcould be either positive or negative. Note that we need the Bfactor multiplying Min order to balance the factor B· |bk+i−1−bw|, which we know is bounded by Bbecause it’s a relative distance inside the berth. 13 Port Optimization 3.4 Formulation analysis The formulation given has a linear objective function. There are common tricks that allow to turn constraints with an absolute value (constraints (13.1) - (13.4)) to linear constraints. Since all other constraints are linear, we just provided a Mixed Integer Linear Programming (MILP) formulation for this version of the BACAP. An exact solver will be able to use Branch and Cut or similar algorithms. For this analysis, we will use has the number of holds per vessel as if it was the same for all of them. The asymptotic analysis works without this assumption if we take h= max hk. The number of variables of the formulation is 2n2+ 2nhQT +nh|G|. The total number of constraints for the problem is n(n−1) 2+ 2(n2−1) + 2n+ 5nh +QTnh + 4nh +TQ +QTh2n2+nh|G|= =−2 + 3n 2+5n2 2+hn(QT +|G|+ 9) + (h2n2+ 1)QT, which grows as O(Q·T·h2·n2) 14 4. Genetic algorithm 4.1 Genetic algorithms overview Metaheuristics are general algorithmic frameworks, often inspired by nature processes, designed to solve optimization problems. Unlike classical heuristics, which are very problem-dependant, the strategy of a metaheuristic can be applied to a variety of problems. Some common metaheuristics are Genetic Algorithms, Ant Colony Optimization, Tabu Search and Simulated Annealing. Metaheuristics fail to provide a lower bound of the objective value of the problems they solve. They only provide a feasible solution, which can be interpreted as an upper bound, so the question is how far this upper bound is from the optimum. Due to the random component of metaheuristics, answers to this question have to be probabilistic. Convergence proofs typically require infinite computation time or memory space [2]. Evolutionary algorithms (EA) are metaheuristics satisfying: they are population based (a set of solutions is considered and improved over iterations) and variation-driven (the way of generating new solutions is introducing variations to the current solutions in the population) Agenetic algorithm is a metaheuristic in the class of EAs inspired by the process of natural selection. It relies on the biologically inspired operators of mutation, crossover and selection. GA works with a set of solutions, which we will refer to as population of individuals. We identify individuals with solutions, but their relation is not necessarily one-to-one. An individual must be able to represent a solution, and this representation is called the chromosome. The chromosome encodes the solution with a set of parameters or features, and each one of this parameters is called gene. Figure 3: Population, chromosome and genes in a genetic algorithm The algorithm keeps a population over a number of iterations, each one of these iterations called generation. In order to create new generations, we have a fitness function that, for every individual, returns the quality of the solution it represents. Then selection is applied based on the fitness of the individuals to select a subgroup of this population called the mating pool. Individuals selected in this group are called parents. Every two parents will generate two new individuals called offspring. By selecting and mating high-quality individuals, there will be higher chances to keep the good properties of the individuals and leave out bad ones, and it won’t allow the bad individuals to generate more bad individuals. In order to widen the explored solutions space and to overcome the limitations of the parents, the algorithm introduces mutations to the offspring genes. The newly generated offspring, along with their parents (to ensure that the we will keep the best individual), replace the current generation. After a fixed number of generations, or until convergence conditions are met, the algorithm stops and 15 Port Optimization returns the best solution of the last population. 4.2 Design and implementation 4.2.1 Encoding of the solutions The encoding of the solutions is a job for the algorithm designer, since it’s very problem-dependant. Common encoding techniques include binary genes, chromosomes containing a permutation (used for routing, scheduling problems...), genes containing numerical or categorical variables... Note that one approach for the individual of our problem could be to just have all decision variables included in the chromosome as genes. zand ˜zwould be binary, whereas b,t,cwould take integer values. We will call it the trivial individual. The trivial individual is O(T·Q·n·h) in size, and the search space would be the same as in the exact formulation. Even though this previous approach could potentially give feasible solutions, the size of the individuals and the equally massive search space are important drawbacks. So let’s discuss the idea of dimensionality reduction in a genetic algorithm and give our design of the individual next. Remember that we are thinking of a solution as an element of the search space (could be feasible or infeasible). Let Sbe the set of solutions and Ithe set of possible individuals, then the encoder is a function e:S→I. Note that the encoding of an individual doesn’t have to be unique. In other words, eis not necessarily injective. This is basically allowing dimensional reduction of the search space, and it’s something that we are going to take advantage of in our design. But dimensionality reduction carries a drawback. Remember that we need a way to evaluate how good an individual is. In other words, given its chromosome information, we will need a function that returns what’s called the fitness of the individual. A natural choice for the fitness function is the negative cost of the solution this individual represents (since we want to minimize cost). For the infeasible individuals, return a minus infinity fitness. Since there might be multiple solutions encoded by the same individual, we would want to return the minus cost of the best solution that is encoded by this individual. To find this best solution, one uses an heuristic or greedy algorithm, that is going to be a function h:I→S. Then the fitness function is just defined as f=−c◦h. Note that one desired property for the heuristic his the following: f(i) = max s∈e−1(i)f(e(s)), ∀i∈I If this condition is not satisfied, then we know for sure that the even in the best case, the algorithm won’t return an optimal solution. Our individuals: We define them with two distinct parts. A first part contains information on the crane assignment. For each job in each hold, there’s a gene with the number of the assigned crane. 16 Figure 4: Crane assignment in the individual In a second part of the chromosome we track the berthing times of the vessels. If we assume that the numbering of the vessels equals the arrival times to the port, then the genes will be a permutation of 1, ... , ngiving the order of berth (as we will see in Section 4.2.5, it won’t necessarily correspond to the berth order in the solution, but simply the processing order of the vessels in our heuristic). Note that we are successfully reducing the dimension of the search space of the problem. Let’s compare our choice with the trivial individual, which belongs to the space {0, 1}n·h·T·Q×[T]3·n, dimtrivial = 2n·h·T·Q·T3·n. On the other hand, our individual has size n X k=1 |Uk|+|Lk|+n, which just grows as O(n·h). It belongs to the space [Q]n·h×Sn, dim = Qn·h·n! Both the reduction in size of the individual and dimension of the space are outstanding. In particular, the small size of the individuals will allow our GA to return solutions in very good times. But as stated before, dimensionality reduction comes with a non-injective encoding function. Therefore we will have to develop an heuristic that tries to recover the best solution encoded in an individual. 4.2.2 Crossover Crossover is the process of taking two solutions (parents) and producing a child solution from them. By recombining portions of good solutions, there’s a chance that we take the best features of each of them to create a better solution than both of their parents. The crossover method has to be chosen to closely match the chromosome’s representation of the solution. In our case, for example, we won’t allow parts of the vessel order permutations to cross with the crane assignment genes. In fact, we won’t crossover the vessel order, and we will just apply the crossover operator for the crane assignment genes. We will use three different methods: one-point, two-point and uniform crossover. •One-point crossover: randomly choose a point in one of the parents and swap information on the left of the point with the other parent. 17 Port Optimization •Two-point crossover: randomly choose two points, cutting the individual into three, and swap parts with the other parent. •Uniform crossover: for each gene in the child, randomly choose from which one of the parents we will take it from. The following figure shows how the one-point crossover would work for our individuals. Figure 5: Crane assignment in the individual 4.2.3 Mutation The purpose of the mutation operator in genetic algorithms is to introduce diversity into the population. Preventing the individuals to become too similar to each other helps the algorithm overcome possible local minima and expands the search space. Let’s discuss the mutation operator used in our genetic algorithm. For the crane assignment part of the individual we use single point mutation. We iterate over the genes and for each one, decide whether we mutate the gene with probability pmutate. The mutation consists in replacing the current crane assigned to the gene’s associated job, by another one with a certain probability distribution. 18 Figure 6: Crane mutations For the order of berth, we similarly iterate over the genes and decide to mutate with probability pswap. Since we want to preserve a permutation, the mutation consists in swapping the selected and contiguous vessels. Since we iterate over all genes, we are actually allowing non-contiguous swaps. Figure 7: Vessel order of berth swap mutations Intelligent mutations: One can always use a uniform probability distribution for the crane mutations. However, in hopes of improving the performance of the GA, we implemented a second mutation strategy for the crane assignment. The idea is that we wanted the cranes that worked the least to be assigned to new jobs, and the algorithm implements that by assigning a new crane to a job with a probability inversely proportional to the total working time of the crane in the encoded solution. In essence, we have pq= 2/Q−wq/W, where wqis the working time of qand W=Pqwq. Motivation for this improvement is backed up with an analysis we did on the total handling times of the cranes. Our hypothesis was that the standard deviation of these times decreased as the solution cost decreased. 19 Port Optimization In the figure below we study these standard deviations of five random problems over a fixed 50 generations of the GA. Although the std function is more unstable, one can clearly identify its decreasing trend pushed by the improvement of the objective function. Figure 8: Analysis of working time distribution over cranes We’ll see in Section 4.4 if intelligent mutations actually have an effect on the performance of the algorithm. 4.2.4 Selection Selection is the stage of a genetic algorithm in which the best individuals are chosen from a population to become parents and generate new individuals using the crossover operation. The best individuals are chosen based on their fitness value, but many strategies to do this selection exist. In our problem, we will use Tournament Selection, which works the following way: from a population of nindividuals we want to select mof them, so we will do mtournament rounds. For each round, choose randomly s<nindividuals, and the one with the best fitness is added to the selected population. sis called the tournament size and affects the selection in the following way: •With a small s, the fittest solutions might not be selected and convergence can be slower. •With a large s, weak candidates have a smaller chance of getting selected causing diversity loss. We will combine Tournament Selection with elitism: the fittest solution will always be part of the next generation in order to prevent the cost to increase over generations. Initial solutions: We will include some feasible solutions in the initial population in order to speed up convergence. These ones are going to be individuals with the trivial berthing order (the arrival order), 20 and with a crane assignment using the following rule: for each hold of each vessel, choose the crane that has the least assigned working time and least duration of the operation in that hold. To do so, we take a relative distance measure of this two parameters with different weights, that are going to give different initial solutions. 4.2.5 Evaluation We discussed in previous sections the need for an evaluation or fitness function. In our case, it’s going to be the cost of a hopefully good solution that the individual represents, since there might be multiple solutions encoded in the same chromosome. Therefore, we need an heuristic strategy to find it, described next. The idea is to process the vessels in the order specified by the individual. For each vessel, we have to first decide its berth allocation, and the genetic algorithm will do that by exploring all the possible berth spots. We’ll do that with depth one, meaning that the best location will be chosen as the one that allows the vessel to leave the earliest (and not taking into account the cost of the complete solution). Once the vessel berths, the cranes assigned to the holds in the vessel do the jobs as soon as they are available. If a crane has to do multiple jobs in one vessel, these rules are applied: do first all the unloading tasks, always going to the closest one. If the same crane has to load a hold that has just been unloaded by it, do it right after. Do then all the loading tasks, always going to the closest one. A pseudo-code would look like this: 4.3 Convergence Since we are dealing with a metaheuristic, when talking about convergence we are not necessarily implying that the point where it converges has any special qualities. We are going to define that the algorithm converges to a point xif xis the best solution in the last Ngenerations. 21 Port Optimization can be accessed after a CPLEX execution. We are going to define the relative gap of any other solution with objective value Sin the same way: |S−LB|/|S|, and it’s going to be useful to evaluate the solutions returned by the GA. 5.2 Mixed Integer Program starts When solving a MIP, one can supply hints to help CPLEX find an initial solution. A MIP start can be a feasible, infeasible, or incomplete solution. If one or more of the MIP starts define a feasible solution, CPLEX sets the best of these solutions as the incumbent. Having an incumbent from the very beginning of B&C allows CPLEX to eliminate portions of the search space and thus may result in smaller trees. Having an incumbent also allows CPLEX to use heuristics that require one [5]. MIP starts are a perfect opportunity to take advantage of the already developed GA. We will feed CPLEX with the best feasible solution found by the GA. The following example helps to get a sense of the amount of time saved time when providing a feasible MIP start. We are working with a small problem (n= 2, B= 7, Q= 3, incompatibilities = 1, fixed cranes = 1), where the genetic algorithm converges to an objective value of 124. We run CPLEX without an MIP start until UB ≤124, and it takes 1 hour and 10 minutes to stop. The returned bounds are UB = 122, LB = 69, 6344. For a fair comparison, we are going to run CPLEX with the 124 cost initial solution until both the UB and LB get better than the ones of the previous execution. After only five minutes, CPLEX has UB = 122 and LB = 76, 1220. The solver takes these five minutes because even though the UB is set directly to 124 after the presolve, the computed LB is only 57, 8236. But as we see here, it’s increased rapidly. In conclusion, for this problem instance, CPLEX with the MIP start takes only five minutes to get better UB and LB than the ones that standard CPLEX found in 1 hour and ten minutes. Note that don’t pretend to generalize the previous result; time reduction is very likely going to be dependant on the problem instance. 5.3 First Come First Serve solution We are going to compare the solutions that our algorithms obtain to a simple heuristic that reproduces the how the berth allocation and crane assignment is done in some ports. In particular, we know that this strategy is used in the Port of Bayside. It’s going to be useful because for the cases where the exact solver is unable to produce lower bounds, we will at least have a way to tell if the genetic algorithm would be valuable for improving the port efficiency. The FCFS approach builds the solutions as follows: •Process the vessels in the order of arrival. •When there’s space available on the berth, do it. The berth allocation uses the same one-depth search as the GA. •Crane allocation works as follows. For each hold, order the cranes from least handling of the cargo in the hold to most. We also keep a list of the assigned working times for the cranes until the process of this hold (ordered from least time to most). Choose the crane whose positions in the two ordered vectors have minimum sum. •If the crane selected for unloading can also load the same hold, assign it there. 28 5.4 Solution of different problem instances In this section we are going to test our algorithms on different problem instances. The problem instances that we will cover are described by: nthe number of vessels, Bthe number of berth sections, Qthe number of cranes, ninc the number of incompatibilities crane-material, nfixed the number of fixed cranes. We will always solve the case of R= 4 materials. The arrival times are randomly generated inside a time window from time 1 to 4·n. The quantity of material to unload or load for each hold is generated as follows: first select two materials to be carried by the ship. For each of these materials, an integer randomly generated from −1 to 5 (negative values are read as zero) is the quantity of the material to load or unload for that hold. The productivity of the cranes, for each material, is generated as a random integer from 1 to 4. If there are incompatibilities, the productivity for a material of a particular crane is set to infinity. For the fixed cranes, their static location is a random berth section. We assume vessels take one time unit from moor to begin of operations and another one from end of operations to departure. We assume there must be a hold of distance between to vessels in port. For each problem instance, we first compute the solution generated by the FCFS heuristic. We then run the genetic mtimes, each of them with different parameters, in order to keep the best solution it converges to. The value of mwill change throughout the different executions. We finally feed this best solution to CPLEX and set two stopping criteria: a maximum of two hours and a half of solving time or an optimality gap of 10% in the best feasible solution. We are going to store the results in a table, whose columns are described next. First five columns are the description of the problem instance. cplex cost is CPLEX’s best feasible solution cost when the stopping criteria is met, using the MIP start given by the GA algorithm. cplex time (min) is the time CPLEX took to solve (not building) the model. cplex total time (min) is the previous column plus the minutes it took to add the constraints in the model building. mip relative gap is r=|S−LB|/|S|, where Sis the feasible solution returned by cplex (cplex cost) and LB is the best lower bound during B&C. ga cost is the cost of the best solution found in the mexecutions of the GA. ga time (min) is the aggregated execution time of the multiple GA runs. ga relative gap is r=|S−LB|/|S|, where Sis the best feasible solution returned by the GA and LB is the best lower bound during B&C. fcst cost is the cost of the feasible solution using the FCFS heuristic approach. ga vs dummy improvement is the improvement in the cost function by using the GA compared to the FCFS solution. cplex vs ga improvement is the improvement in the cost function by feeding the GA solution to CPLEX. 29 Port Optimization Table 2: Multiple executions of FCFS, GA and CPLEX Next we are going to give some remarks on the results of the Table 2. In the first place it’s consistent with our model that for a fixed problem, increasing the number of cranes allow us to find schedules with less cost. For instance, see problems 6-8 or 10-14. The same happens with fixed cranes and crane-material incompatibilities, which worsen solution costs, as we can see in 14-16. The berth length is also a factor that makes makes solutions better, like in problems 33 and 34. Surprisingly, a formulation with more constraints and variables doesn’t necessarily translate to larger solving times. In problems 8-13, adding cranes turns into much shorter solving times. In 20, 21 we can see that fixing the position of one crane also reduces the solve time significantly. The intuitive explanation here is that adding cranes makes it easier to define good schedules, and having fixed cranes reduces the degrees of freedom of the problem. For the problems where we were able to use CPLEX, the optimality gap for the GA is on average 22%, which is fairly good considering the big role that heuristics play in our algorithm (see dimensionality reduction in 4.2.1). CPLEX has been able to improve most of the solutions returned by the GA. On average, with a time 30 limit of two hours and a half CPLEX improves the GA solution by 6%. Taking into account that in most cases it reaches the optimality gap before the time limit, letting CPLEX improve the given GA solution for those problems is the reasonable choice. And even though we don’t reach the desired gap, at least CPLEX will provide a lower bound on the objective function that will be insightful when it comes to validating the solution given by the GA. The entries of the table with blank CPLEX columns values (36 - 42) correspond to those that were intractable by the solver. CPLEX aborted the solve function after some time due to Error 1001: Out of Memory. Even thought we don’t have evidence for this, the performance of the GA on the smaller problems allow us to be at least optimistic about the GA’s performance on these larger and intractable problems. Besides, the GA consistently improves the FCFS solution (24% on average), which makes it relatively valuable as a solution method. When it comes to execution times, GA takes approximately 5 minutes per execution, which is significantly less than the solve time that CPLEX requires to reach the stopping conditions. Besides, the GA uses most of the time to improve the current solution it has, in a small percentage. In other words, convergence speed slows down rapidly after the first generations. In fact, in the first generation, the GA usually outperforms the FCFS solution. We can see an example of this behaviour in Figure 8. But we shall not make the mistake of concluding that the GA is faster than CPLEX in any way, because the solver CPLEX returns more information on the problem than the GA (returns LB and UB). We could think that GA is much faster at returning upper bounds, if that is what we are interested in. Not an extensive effort was put in trying to get over the Out of Memory error because it’s likely that the root of the problem is the complex formulation of it. There are techniques that allow the solver to use more computing power from the PC used, but we decided to not look into that direction in this project. If we wanted to implement the exact solution in a production environment we would already have a notable improvement in the available computing power. It’s difficult to estimate whether CPLEX will be able to handle a particular problem or not. However, it only solved simple instances of problems with n≤6. Increasing the number of holds per vessel naturally makes the problem more difficult. We always have reached the ’limit’ of the problems that CPLEX can solve by increasing the number of holds. In other words, by leaving n,B,Qand the other two parameters fixed and increasing the number of holds per vessel we always ended up with a problem where CPLEX runs out of memory. 5.5 Optimizing the Port of Bayside In this section we are going to test the Genetic Algorithm on the Port of Bayside data. The port had the following particularities: •Vessels served in the port are between 170 and 200 m long. •The berth area is 880 m long. We were told that at most four vessels could fit in the port at the same time. •There were two types of vessels. First type only carries containers, and both unloads and loads containers in the port. Second type carries bulk materials: scrap, grain and coal. Since the productivity of the cranes is independent on the type of bulk material, we will consider it as one material. •There were seven mobile cranes and one fixed crane in the middle of the berth area. Mobile cranes can only handle bulk whereas the fixed cranes only containers. 31 Port Optimization Using this information, we will model model the optimization with the following assumptions. Assume that the smallest vessels have three holds. Fitting four of this vessels in the berth and leaving some safety distance between them gives a B= 15, which comparing to the actual length of the port we have holds of length 59 m. It is a great assumption since this in our model have length 3 ·59 = 167 m. We will have some vessels with 4 holds as well because in the real problem, certain vessels could be up to 220 m long. We will R= 2 materials (bulk and containers), and a set of |Q|= 10 cranes with the following assumptions: cranes 1, ... , 7 will be mobile, with an incompatibility to handle the material 2 (containers). 8, 9, 10 ∈Gwill be fixed in the positions b= 6, 7, 8. These three cranes actually correspond to one in the real port, with the particularity that the real one is big in size and can rotate its lift arm so it can have access to the holds i−1 and i+ 1. We will adjust the productivity of the three cranes so it corresponds to the real one, and add an incompatibility with the material R= 1 (bulk). We consider that cranes move along the berth with speed 1. For instance, they take 30 minutes to move from ito i+ 1, and it makes sense if we consider that we there might be shifts of operators, and that the cranes take some time to prepare before starting the actual operation. We will discretize time in steps of length 30 minutes. We will assume that the mobile cranes can handle one ’amount’ of bulk in thirty minutes. In the real port, each of these ’amounts’ is well defined by a constant tonnage measure. Each hold of a bulk carrier contains the amount of material given by a random integer between 1 and 3. When it comes to containers, we will assume that the fixed crane can handle one ’amount’ of containers in thirty minutes. Recall that we were modelling the big fixed crane of the Port of Bayside as three cranes in our model. So it makes sense that every hold carries the same amount of material, so no different jobs will start or end at different times on different holds of the same container carrier. Due to the confidentiality of the data, we can’t go into the details of how we got to the previous assumptions. We are going to plan the operations to serve n= 10 vessels. By looking at the total time-span of the operations, we would be planning approximately two 48 hours, but in reality there are no operations scheduled at night, so it would correspond to an approximate of 4 days. Below we give the schedule for the schedules returned by the FCFS solution and the genetic algorithm, with a cost of 102 and 72, respectively. 32 Figure 16: Space-time of the FCFS heuristic for the Port of Bayside Figure 17: Space-time diagram of the GA solution for the Port of Bayside At first sight, there does not appear to be an obvious change to the FCFS schedule that improves the cost, but using the GA we are able to provide a more efficient allocation of the vessels, and a feasible schedule that reduces the cost by 29%. In particular, GA introduces changes in the berth order: vessels 8, 9 and 10 are reordered. The crane assignment for the other vessels is slightly modified as well. In this particular example, we reduce the waiting time (time compressed between the vessel arrival and berthing time) from 19 to 13 (32%), and this turns into various cost reductions: emissions, fuel, wages, harbor fees... Accenture owns data about the fuel consumption of cargo vessels when they are in standby (in our case, waiting outside the port), while unloading/loading, or while anchored. A very simple computation gives a sense of the magnitude of the cost reduction; we are reducing the fuel consumption from 22.100 to 17.950 liters (28%) for that particular example. 33 Port Optimization Running on 5 more problems with data generated by some other random seeds, we get an average cost improvement of 22%. According to Data Scientists at Accenture, it’s a decent improvement and it’s likely that clients would be interested in moving on with the project to perfect the algorithm and implement the solution. 5.6 Conclusions In this research we provided a MILP formulation for an extended version of the BACAP problem that considered distinguishable cranes, different crane-material productivity, allowed fixed cranes, and takes into account the time that cranes take to change positions inside the port. We developed and tested a Genetic Algorithm that provides feasible solutions of the problem in relatively short times. This algorithm improves an heuristic First Come First Serve method by an average of 24%, even in medium and large sized problems. We finally programmed CPLEX to solve the MILP formulation. We first checked that providing an initial solution could notably reduce the solve time. For some small-sized problems, CPLEX was able to improve the solution of the GA and provided an optimality gap that we could use to evaluate the performance of the GA: the developed metaheuristic returns solutions that are at least 22% on average from the optimal solution, giving us hope that we could have a similar percentage on larger problems. In conclusion, we developed two solution methods that could be used for different purposes. •We have a MILP formulation that can be solved by commercial solvers like CPLEX. We refined this exact approach by providing CPLEX a start solution returned by the GA. This method is exact in the sense that we also compute a lower bound for the objective function that allows us to have an optimality gap. However, we could only have this method to efficiently work for small instances. •We have a Genetic Algorithm that provides feasible solutions of problems of all sizes in very short times, and we know that it notably improves simple approaches like the First Come First Serve rule. It could be used in a more ’naive’ way, when there are no optimality requirements in the solution and the goal is to improve the current scheduling method used. It could also provide insights in a simulation-like project. For example, running multiple trials of the GA can be useful in order to detect bottlenecks, optimize physical characteristics of the port (optimal number of cranes, optimal size of berth...), or A/B test different logistic strategies. 5.7 More ideas and further research A possible way to deal with the BACAP problem that hasn’t been studied in the literature is the following. Given the arrival times list {a1, ... , an}, we would like to come up with i1, ... , irvessel numbers and consider the subsets of vessels {1, ... , i1−1},{i1, ... , i2−1}, ... , {ir, ... , n}. The main idea of this approach is that we would solve the previous sub-instances independently. Their reduced size would be suitable for the use of an exact solver (for example our GA + CPLEX). Some remarks on this approach. We don’t want rto be too large, since setting a break on vessel i implies that this vessel is going to berth after all vessels j<ihave departed. Naturally, the ideal case for a break is the following: we were able to solve exactly the BACAP for 1, ... , i−1 and max{c1, ... , ci−1}<ai. The thing is that checking if this condition is met for every iis too costly, so some other strategy should be designed in order to provide i1, ... , ir, for example by studying the cargo of the vessels and the gaps 34 between arrival times. Note that for some problem instances, we could ensure that our algorithm solves to the optimum, but in some other cases, the previous strategy turns into an heuristic. Another possible line of research to give continuity to this project could be doing a lower bound analysis, similar to the work of Aykagan Ak [1]. Providing a lower bound of the objective cost in polynomial time would allow us to compute optimality gaps for solutions coming from heuristic methods for large instances, or to implement a custom version of Branch and Cut using our lower bounds instead of solving linear relaxations. An extension of our problem could add the no-crossing constraint, that models the case of ports where cranes go on rails, so they can’t cross each other. Just adding this constraint to our problem would very likely make a lot of problems infeasible because we consider crane-material incompatibilities, and not all sections of the berth would be reachable by all cranes. It’s more natural to consider this constraint on a problem where all cranes are assumed indistinguishable. However, we checked that it would be possible to add the constraint for our problem in a linear way. 35 Port Optimization References [1] Aykagan Ak. Berth and quay crane scheduling: problems, models and solution methods. PhD thesis, Georgia Institute of Technology, 2008. [2] Leonora Bianchi, Marco Dorigo, Luca Maria Gambardella, and Walter J. Gutjahr. A survey on metaheuristics for stochastic combinatorial optimization. Natural Computing, 8:239–287, 2009. [3] Christian Bierwirth and Frank Meisel. A survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research, 202:615–627, 2010. [4] Christian Bierwirth and Frank Meisel. A follow-up survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research, 244:675–689, 2015. [5] IBM ILOG Cplex. V20. 1: User’s manual for cplex. https://www.ibm.com/docs/en/icos/20.1.0. Accessed: 2021-12-30. [6] Giovanni Giallombardo, Luigi Moccia, Matteo Salani, and Ilaria Vacca. Modeling and solving the tactical berth allocation problem. Transportation Research Part B: Methodological, 44:232–245, 2010. [7] Xiao-le Han, Zhi-qiang Lu, and Lifeng Xi. A proactive approach for simultaneous berth and quay crane scheduling problem with stochastic handling time. European Journal of Operational Research, 207:1327–1340, 2010. [8] Akio Imai, Hsieh Chia Chen, Etsuko Nishimura, and Stratos Papadimitriou. The simultaneous berth and quay crane allocation problem. Transportation Research Part E Logistics and Transportation Review, 44:900–920, 2008. [9] Cagatay Iris, Dario Pacino, Stefan Ropke, and Allan Larsen. Integrated berth allocation and quay crane assignment problem: Set partitioning models and computational results. Transportation Research Part E: Logistics and Transportation Review, 81:75–97, 2015. [10] Maciej Machowiak, Ceyda O˘guz, Jacek Blazewicz, and T. C. E. Cheng. Berth allocation as a moldable task scheduling problem. 2004. [11] Frank Meisel and Christian Bierwirth. Heuristics for the integration of crane productivity in the berth allocation problem. Transportation Research Part E: Logistics and Transportation Review, 45:196–209, 2009. [12] United Nations Conference on Trade and Development. Review of the maritime transport. 2018. [13] Young-Man Park and Kap Hwan Kim. A scheduling method for bert and quay cranes. OR Spectrum, 25:1–23, 2003. [14] Birger Raa, Wouter Dullaert, and Rowan Van Schaeren. An enriched model for the integrated berth allocation and quay crane assignment problem. Expert Systems with Applications, 38:14136–14147, 2011. [15] Mario Rodriguez-Molins, Federico Barber, Maria Sierra, Jorge Puenta, and Miguel A. Salido. A genetic algorithm for berth allocation and quay crane assignment. Progress in Artificial Intelligence, 2:177–192, 2014. 36 [16] Augusta Saraiva and Brendan Murray. California’s busiest port is straining the global supply chain. Bloomberg, 2021. [17] St´ephane Zampelli, Yannis Vergados, Rowan Schaeren, Wout Dullaert, and Birger Raa. The berth allocation and quay crane assignment problem using a cp approach. 2013. [18] Canrong Zhang, Li Zheng, Zhi-Hai Zhang, Leyuan Shi, and Aaron Armstrong. The allocation of berths and quay cranes by using a sub-gradient optimization technique. Computers & Industrial Engineering, 58:40–50, 2010. 37