The multiport container ship stowage problem consists in determining the position of the containers on board a ship along its route with the objective of minimizing the number of unproductive moves required in the loading and unloading operations at each port. This paper presents an integer programming model for the problem and proposes several sets of valid constraints that bring its LP-relaxation closer to an integer solution. Moreover, it presents a GRASP algorithm that generates stowage plans with a minimal number of unproductive moves in a high percentage of medium and large-size instances. An extended computational analysis has been performed in which, to the best of the authors’ knowledge, the efficiency of integer programming models for the problem is tested for the first time. With respect to GRASP, the computational results show that it performs well on different sized datasets.

1. Introduction

Seaborne trade accounts for nearly 80% of world trade and is one of the most efficient modes of transportation [1]. About 90% of nonbulk cargo worldwide is transported by container [2]. Since its development, containerization has played an everincreasing role, due to its important advantages, reducing packaging and deterioration of goods and increasing productivity in the intermodal transportation of products. Global containerized trade expanded at a rate of 3.1% in 2016, with volumes attaining an estimated 140 million Twenty-foot Equivalent Units (TEUs) [1].

The high volumes reached have produced an increase in ship size, which still continues. Modern container vessels can handle up to 20,000 TEUs. Currently, the world’s largest container ship, the OOCL Hong Kong, has capacity for 21413 TEUs [5]. For these large vessels to achieve economies of scale, they must be loaded at near full capacity. These vessels sail from port to port loading and unloading thousands of containers. Minimizing the time a vessel spends in port requires an efficient stowage plan, a plan that describes in which position each container should be loaded into the vessel. This plan has to satisfy high-level constraints, related to the stability of the ship, as well as low-level constraints related to the way in which containers have to be stacked, the weight distribution, and the specific conditions regulating containers with dangerous products. The main objective is to load all the assigned containers on board if they feasibly fit into the ship. If this is not possible, then the objective is to load the maximum number of containers. A second very important objective is to minimize the number of unproductive moves in the loading/unloading operations. These moves, also called reshuffles, arise at a port when some containers going to later ports have to be unloaded to gain access to containers that must be unloaded at that port and are situated below them. Once all the containers going to this port have been unloaded, the reshuffled containers are loaded again into the ship, together with all the containers that had to be loaded at this port.

In this paper we study the multiport container ship stowage problem, considering the whole route of the ship and the different sets of containers which must be loaded/unloaded at each port on the route. The problem consists in determining the position of the containers loaded at each port so as to minimize the number of unproductive moves necessary to unload the containers at their port of discharge.

The contributions of our study are basically of two types. Firstly, we have developed a new integer linear model, as an alternative to the existing models by Avriel et al. [3] and by Ding and Chou [4], requiring fewer variables and obtaining better results concerning the number of optimal solutions obtained and the quality of feasible solutions when the time limit is reached without proving optimality. Nevertheless, as the computational study shows the limits of mathematical models for solving large-size instances, in the second part of the study, we have developed a metaheuristic GRASP algorithm, going one step further than the constructive algorithms proposed so far, and including randomization strategies and a local search tailored to the problem. An extensive computational study shows that GRASP obtains much better results than previous procedures, generating stowage plans with a minimal number of reshuffles in a high percentage of medium and large-size instances.

The structure of the paper is as follows. The relevant literature on stowage problems is reviewed in Section 2. Section 3 presents a detailed definition of the problem. In Section 4 we present a new mathematical formulation and propose several families of valid inequalities. The GRASP metaheuristic is presented in Section 5 by describing its randomized constructive algorithm and the local search. The computational experience carried out is summarized in Section 6. Finally, conclusions and future work are highlighted and discussed in Section 7.

2. Previous Work

Academic work on stowage planning problems has increased notably in recent years. The studies can be classified into two main groups: single-port and multiport problems.

The single-port problem consists in determining the arrangement of containers in the ship at a given port without considering the loading of containers in subsequent ports. The vast majority of papers in the literature related to stowage planning belong to this group and it is usual to refer to this problem as the Master Bay Planning Problem (MBPP). For a detailed description of the MBPP and its constraints, see, for instance, Ambrosino et al. [6]. Even though this problem is difficult, it is not as hard as the multiport problem, so there is a trend of progressively including more realistic constraints. Papers related to the MBPP can be divided into two approaches: those that consider the problem as a whole and those that decompose it into several optimization subproblems. In the first approach, exact methods such as mathematical models [6] or branch and bound algorithms [7] can be found, as well as heuristic algorithms, such as Sciomachen and Tanfani [8], based on an exact method for the 3D-BPP. In the second approach, the problem is divided into two phases. First, the containers are distributed in clusters along the vessel. Then, for each cluster a specific position for each container must be found. Examples of this approach are Ambrosino et al. [9], Pacino et al. [10], Delgado et al. [11], and Parreño et al. [12].

In the multiport problem, the stowage problem is considered as a dynamic loading and unloading problem. According to the classification by Lehnfeld and Knust [13], it is a combined problem, since it is possible to perform two types of moves. Two main approaches have been followed in recent years. On the one hand, some authors use a multiphase approach, proposing a decomposition of the problem into several optimization models that can be solved individually. This decomposition makes it possible to include realistic characteristics of the problem, in the objective function, considering costs related to nonloaded containers, reshuffles, and quay cranes, and also to include constraints concerning the stability of the ship [1418]. On the other hand, there are several studies that follow a single-phase approach. The problem is solved as a whole, looking for loading strategies that minimize the number of reshuffles at each port on the route, thus providing a complementary view of the stowage problem, usually considering routes with many ports and ships with large numbers of containers, but disregarding stability conditions and considering just one type of container. Avriel et al. [3] presented a mathematical model and developed a constructive heuristic, the Suspensory Heuristic (SH), to stow containers of the same size in a ship consisting of a single rectangular bay so as to minimize the total number of unproductive moves in the whole route. The same problem was studied by Dubrovsky et al. [19], who presented a genetic algorithm that was able to incorporate constraints beyond the simplified problem and performed nearly as well as the SH. More recently, Ding and Chou [4] have proposed a new heuristic algorithm that is shown to perform better than the SH in an extensive computational study, and also, in an on-line supplement, an integer linear model. Our study belongs to this class of single-phase approaches to multiport stowage problems, and in addition to a new integer linear model, we develop a more efficient GRASP metaheuristic algorithm.

3. Problem Description

Let us consider a ship traveling along a trade route consisting of ports. All containers are the same size and no stability conditions are considered, so the ship can be seen as having a single bay with rows and columns. The position of each container is described by the pair where , is the bottom row, and , is the first column on the left. We consider sets , , and to simplify notation. Note that here we are following the notation proposed by Avriel et al. [3] and Ding and Chou [4], although in other stowage problems horizontal levels are referred as tiers and vertical levels as rows or stacks.

The loading/unloading operations along the route are defined by an transportation matrix whose input is the number of containers originating at port with destination port . It is a nonnegative upper triangular matrix, so for all . We assume, as has been done in previous studies [3, 4], that all transportation matrices are feasible, so all the containers to be shipped can be stowed at every port according to (1).

The vessel starts its route at port , and sequentially visits ports . At each port , containers with destination are unloaded and containers bound for destinations are loaded. A container whose destination is port is called a -container. At the last port , the vessel is emptied. Therefore, it can be seen as a circular route in which port 1 is the same as port . If there is a -container in position and there is a -container, , in some position with , the first one is said to be a blocking container and will produce a reshuffle in some port , . The objective of the problem is to minimize the number of reshuffles along the route. Note that a non-full loading transportation matrix can always be transformed into a full-loading matrix , adding as many dummy containers from to as necessary. We can consider, without loss of generality, the inequality of (1) as an equality.

4. Mathematical Formulation

In this section we describe a new integer linear model and several classes of valid inequalities to enhance it.

4.1. An Integer Linear Model

The model involves two types of binary variables describing the configuration of the ship and the containers that are moved unproductively at each port.

According to the problem description, the ship is emptied after the unloading operations at port . Therefore from port to port there are only -containers in the ship, and there will be no reshuffles at port . Moreover, since we assume that the number of containers to be loaded on board is not greater than the number of available locations, the exact position of each container when the ship leaves port is irrelevant. Under these assumptions, neither variables nor variables need to be considered.

The multiport container ship stowage problem can be formulated as follows:

The objective function (3) seeks to minimize the total number of reshuffles. Reshuffles are identified by constraints (4) and (5). Constraints (4) identify changes of destination in a given position. Constraints (5) identify unproductive moves required to unload a container, serving two purposes. On the one hand, to take into account the obvious case in which if , meaning that there is an -container in slot at port , the container above it in slot must also be unloaded at port (making ) or there is a reshuffle and one of the variables must be 1 for some . On the other hand, if there has been a reshuffle in the stack at any tier (and then the variable is equal to 1 for some ) the container above it in slot must be unloaded at port or it will produce a new reshuffle (and some variable must be 1 for some ). In this second case, all blocking containers in the stack will be counted as reshuffles. Finally, constraints (6) and (7) fix the number of containers to be loaded at each port and constraints (8) force each location to have at most one container at each port. At the first port, the ship is empty before loading operations, and after these operations there are -containers on board (constraints (6)). At other ports, the number of -containers loaded at each port , , is equal to the number of -containers on board after loading operations at this port minus the number of -containers on board before loading operations at that port, that is, minus the total number of containers whose destination is port and were loaded at a port before (constraints (7)).

We refer to the example in Figure 1 to better understand the model and detect the different types of reshuffles that can appear. It shows the configurations of a solution on ports 2 and 3. Each square represents a container and the number indicates the destination port. The gray-marked containers are those that are relocated on port 3. As there is a change of destination in position between ports 2 and 3, there has been a reshuffle at port 3 of that container and by constraints (4). We now focus on the 6-container located in slot at port 2. It must have been unloaded at port 3 to unload the 5-container but there is a container with the same destination in slot at port 3, so constraints (4) do not consider this case. It is identified by constraints (5). Since there has been a reshuffle in the stack 1 at tier 3 at port 3, the container above it in slot must have destination or it will produce a new reshuffle as is the case. Constraints (5) also identify unproductive moves required to unload a container at its destination port. Consider now the container located in slot at port 2. There is a 3-container at port 2 in slot , so the container on top of it has to be unloaded unproductively at port 3 because its destination is not port 3.

4.2. New Valid Inequalities

The integer linear formulation described above produces valid solutions for the multiport stowage problem, but its linear relaxation is quite loose, producing very low lower bounds. In most cases, the lower bound at the root node is precisely zero. In order to enhance the integer formulation, we studied the solutions provided by its linear relaxation and identified some valid constraints satisfied by all integer solutions but not by the linear relaxation solutions. Adding these inequalities to the initial model would produce a tighter linear relaxation, improving the quality of the bounds and therefore reducing the search tree required to obtain an optimal solution.

We consider the upper row of the bay at port and the number of containers going to the next port, . If the number of slots in the upper row occupied by -containers at port is less than the minimum number required to ensure direct access to these containers, there will be at least as many reshuffles as necessary at port to reach that number. It is satisfied by adding constraints (9).

Constraints (10) extend (9) to consider not only row , but all rows from to .

We now consider the bottom row. If there are containers on board going to the last port, , either they are placed on top of other -containers or in the bottom row. Otherwise, they will block other containers and that will produce some reshuffles at subsequent ports.

Instead of considering only the containers going to the last port, , we can consider the containers going to a port and subsequent ports. The argument is similar to the previous constraints, but we now have constraints (12).

If after loading operations at port there is a container in position whose destination is port or a port subsequent to on top of a container whose destination is previous to , the first container will be relocated at some port prior to . Then, constraints (13) hold.

When we introduce these constraints, the linear programming relaxation provides a solution which lies in general not too far away from the integer optimum. More importantly, the lower bound is increased.

We illustrate their effect by an example, in which we consider a container ship with 6 rows and 4 columns and whose route is composed of 8 ports. Although its optimal value is 3 reshuffles, the objective value of the LP-relaxation (without valid constraints), provided by CPLEX at the root node of its internal Branch & Bound process, is zero, which is indeed a loose lower bound. Figure 2 shows the configuration provided by the LP-relaxation at port 2. On the left-hand side, the bay is represented by the large square, divided into smaller squares that represent each slot of the bay, and each one of these squares is split up (using dotted lines) into smaller squares in order to assign the value taken by the variable that represents each destination port.

Looking at the two enlarged slots that represent slots and and considering port as a reference, the variable whose destination is port 8 takes a positive value in the upper slot. Nevertheless, in the lower slot the variable whose destination is port 8 takes value zero, so constraints (13) are not satisfied. Moreover, according to the transportation matrix, there is one container going to the next port, port 3, so in order to avoid reshuffles, it should be stowed in the top row, row 6. However, looking at the top row, at the top right corner of Figure 2, the sum of the values of variables whose destination is port 3 is less than one. This solution could be removed by including constraints (9).

If the valid constraints are included in the initial model and the LP-relaxation of this extended model is solved, its objective value takes value two, so the lower bound is increased and the solution obtained is very near to an integer solution.

5. A GRASP Algorithm

GRASP (Greedy Randomized Adaptive Search Procedure) is a multistart metaheuristic in which each iteration consists basically of two phases: a randomized construction phase and an improvement phase. The first phase builds a feasible solution, and the improvement phase seeks to improve the outcome. Once the two phases have been executed, the solution is stored and the method starts another iteration, keeping the best solution found so far. An extended description of the algorithm can be found in Resende and Ribeiro [20]. GRASP algorithms have been successfully developed for many hard combinatorial problems, in container loading [21], location [22], and routing [23], among other areas.

Concerning stowage problems closely related to that studied here, Parreño et al. [12] also used GRASP for the slot planning problem, but in recent papers other approaches have been successfully applied: particle swarm optimization [24], genetic algorithms [2527], simulated annealing [27], or hill climbing [27].

The following subsections present the specific implementation of the metaheuristic GRASP that we have developed for this stowage problem: the randomizing strategies, and the local search procedure.

5.1. Randomizing the Constructive Algorithms

To the best of our knowledge, there are two constructive algorithms in the scientific literature that deal with the problem studied in this paper: that of Avriel et al. [3] and that of Ding and Chou [4]. The suspensory algorithm proposed by Avriel et al. [3] is basically composed of a hanging step and a filling step. In the hanging step the containers are assigned to columns by increasing order of destination, occupying the top positions (hanging), and in the filling step the columns are completed and containers are put in the usual positions one on top of the other. On the other hand, the algorithm by Ding and Chou [4] is composed of rules. Containers by decreasing order of destination are assigned one by one following the corresponding rule in each case. The computational results provided by Ding and Chou [4] show that their algorithm outperforms that of Avriel et al. [3] and obtains good solutions for the large-size instances tested.

Unlike other implementations of GRASP, the randomization procedure used here does not randomize the selection of the next element to be included in a partial solution, but the constructive process itself. The two constructive algorithms can be randomized, taking into account the way in which they build a feasible solution.

In Avriel et al. [3] the positions for the containers are assigned column by column. We randomize the next column to be assigned at each step and the order in which containers are loaded at port 1. Since the algorithm places the containers by increasing destination, the randomization consists of using a geometric distribution of parameter to choose the sequence. We also randomize the number of containers to be loaded in the Hanging Stage (H), and the containers not assigned in this phase are considered leftover containers to be loaded in the Filling Stage (F).

In Ding and Chou [4] the assignment is made according to the destination of the containers, group by group, following a set of rules. We randomize the algorithm using two strategies. On the one hand, we modify Rule 4. In the original algorithm, if the destination of the container to be loaded is the last port (), and the set of empty columns and columns partially filled only with -containers, denoted by , is not empty, Rule 4 assigns -containers to the column with the largest number of empty slots in set until no such loading moves can be carried out. We modify this rule, assigning the -containers to columns selected at random from . On the other hand, we randomize the order of some rules followed in Step 4. If , there are three rules, Rules 5, 6, and 7, that assign containers to columns in such a way that there will not be any blocking containers. These rules are used sequentially. In the randomized algorithm the rules are taken in random order.

5.2. Improvement Phase

In the improvement phase, the solutions of the constructive phase are modified to find better solutions. When defining a move, we have to take into account the specific structure of the multiport problem, in which decisions at one port directly determine the layout at later ports. Looking for a move that was simple to compute and, at the same time, able to modify the solution substantially, we realized that loading operations at port 1 play an important role in the final configuration of the stowage plan. The initial configuration largely determines the position of the empty slots after the unloading operations at the subsequent ports. Therefore, we have developed a local search algorithm based on this idea in which only the configuration at port 1 is subject to changes. At this port, we select two columns and , with . Containers belonging to are grouped according to their destination port and we define as the set formed by all the -containers in column , . Let and be the sets of containers placed in columns and respectively.

A neighboring solution can be obtained by exchanging a homogeneous subset of containers belonging to set with another subset of the same cardinality, but not necessarily homogenous, from set . A subset is said to be homogeneous if it is composed of containers with the same destination. As containers with the same destination tend to appear together in the stowage plans, it seems natural to select homogeneous subsets from column to be moved to other column . However, it would not be easy to find homogeneous subsets of the same cardinality in and therefore we allow any subset in this column to be included in the move. The exchange is carried out according to the following steps:(1)Take containers from any group in column with , . This set is denoted as .(2)Select another set of containers in column with destinations different from . This set does not have to be homogeneous and it is denoted as .(3)Let be the new set of containers to be loaded in column .(4)Let be the new set of containers to be loaded in column .(5)Place containers of in column and of in column in decreasing order of destination. Note that the remaining columns maintain their initial configuration.(6)After the exchange between the two columns, complete the solution with one of the constructive algorithms described above.

Figure 3 shows on the left-hand side the configuration at port 1 of one solution for an instance in which , , and . Each square represents a container and the number indicates the port of destination. The figure also shows the possible exchanges to be considered if , , and two containers from the group of 9-containers of () are chosen.

The local search algorithm starts from a candidate solution and then iteratively moves to a neighboring solution following a first-improving strategy. The search starts with and all possible swaps with are checked until an improved neighboring solution is found. If an improvement is not obtained, it proceeds with and , and continues until all pairs of columns have been considered.

Algorithm 1 shows the pseudo-code of the procedure for exchanging containers between columns and of a solution . The objective function is the total number of reshuffles. If Algorithm 1 returns , it has not been able to obtain a better solution. If it finds a better solution , the search starts again with the new solution as the initial solution and with columns and . In line 6 COMB determines all the ways of selecting containers belonging to column with destinations different from . For each selection, , a new neighboring solution is built.

(1) function  EXCHANGES
(2) for    to    do
(3) if    then
(4) for    to    do
(5) ;
(6) COMB;
(7) for all    in  COMB   do
(8)   SWAP;
(9) if     then
(10) Return  ;
(11) Return  

6. Computational Experiments

To test the performance of the mathematical model and the heuristic and metaheuristic algorithms, exact solutions were obtained by CPLEX 12.7 considering a time limit of 3600 seconds and heuristic algorithms were run with a time limit of 300 seconds or 1000 iterations. We conducted all computational experiments on virtual machines with 4 virtual processors and 16 GBytes of RAM running Windows 10 Enterprise 64 bits. Virtual machines are run in an OpenStack virtualization platform supported by 12 blades, each with four 12-core AMD Opteron Abu Dhabi 6344 processors running at 2.6 GHz and 256 GB of RAM, for a total of 576 cores and 3 TBytes of RAM.

6.1. Test Instances

We used two sets of instances in which full transportation matrices were randomly generated using the Authentic Matrices Generation Algorithm proposed by Ding and Chou [4]. For each set of parameters, (ports), (rows), and (columns), five different instances were generated. On the one hand, as the instances generated by Ding and Chou [4] were not available, we generated Set I looking for large-size instances. On the other hand, Set II was generated looking for smaller but more computationally difficult instances. Studying Set I we observed that when the number of columns was increased, the problems became easier, with an optimal solution of zero reshuffles in most cases, making it difficult to assess the relative efficiency of the algorithms. Therefore, we decided to generate Set II with fewer columns but more rows and ports so as to create instances in which the optimal solutions have many reshuffles, in order to obtain more information about the performance of the algorithms.Set I: 625 instances.Set II: 840 instances.

6.2. Performance of the Integer Linear Programming Model

All combinations of valid inequalities were tested in order to decide the best configuration of the integer linear model proposed. In the light of the results, we decided to use the basic model (expressions (3)-(8)), as well as valid inequalities (10) with .

We compared the performance of our mathematical formulation with the two existing models described in Avriel et al. [3] and Ding and Chou [4] by implementing and running all of them with the default CPLEX, with a time limit of 3600 CPU seconds and using 4 threads. Table 1 shows the results obtained by the three models on the second set of instances, Set II, whose main characteristic is that instances are smaller and have a greater number of reshuffles. The instances in Set I are too large for mathematical formulations. Although all models optimally solve instances of four ports, with zero reshuffles, they are not able to obtain an optimal solution when the number of ports increases.

In Table 1 all instances of 4 and 6 ports in Set II are optimally solved by the three models with averages of 0.1 and 0.6 reshuffles. For 8 ports, Avriel’s and Ding and Chou’s models provide the same number of optimal solutions and our model reaches 3 more optimal solutions. For larger numbers of ports, the new model outperforms the other models in terms of both the number of optimal solutions and the average number of reshuffles, with differences increasing as the number of ports increases. However, although the average number of reshuffles in our model is reduced by 74.8% compared to that of Avriel’s model, it is still too high, because when the optimal solution is not found, the average number of reshuffles of the best solution found is sometimes very high. Over the 473 instances optimally solved by the three models, our model achieves a reduction in CPU time of 50.42% compared to Avriel’s model and of 75.9% compared to Ding’s model.

6.3. Selecting the Best Configuration for the GRASP Algorithm

We first test the constructive algorithms developed for the multiport container ship stowage problem and their randomized counterparts developed in this paper. The GRASP presented here applies the best randomized version in the construction phase and the best constructive approach to complete the solutions in the improvement phase. We then study the effect of the local search by itself on the sets studied and we finally establish the stopping criteria used.

6.3.1. Performance of the Constructive and Randomized Algorithms

We reimplemented the two constructive algorithms proposed by Avriel et al. [3], SH, and Ding and Chou [4], DCA, and their randomized versions, in C++. In the SH algorithm, there are two alternative column selection rules, the function rule and the necessary shift rule. The authors tested these two rules and established that the function rule is better than the necessary shift rule if and , and is worse in other cases. Therefore, we apply these two alternatives following that criterion.

As can be seen in Table 2, the average numbers of reshuffles obtained by the deterministic algorithms are very small in all cases in Set I and increase in Set II. The greater the number of ports and the smaller the number of columns, the more difficult the problem. The results shown in the second part of Table 2 underline this statement. It can be seen that the average numbers of reshuffles for groups of 4 to 16 ports in Set II are higher than those for groups with a similar number of ports in Set I. Although the SH algorithm finds more solutions with 0 reshuffles than the DCA algorithm, solving 465 instances compared to 264 in Set I and 177 compared to 146 in Set II, DCA reduces the number of reshuffles by 56.3% and 11.1% in Sets I and II, respectively. Therefore, the constructive algorithm by Ding and Chou [4] outperforms that by Avriel et al. [3] for both sets.

With respect to the randomization strategies, Table 2 also shows the values of the solutions obtained by each randomized algorithm for the instances in both sets. It can be observed that in all cases the average number of reshuffles decreases compared to their deterministic counterparts. The randomized version of the SH algorithm reduces the average number of reshuffles given by its deterministic version by 33.9% and 39.3% on Sets I and II, respectively. On the other hand, randomizing the order of the steps followed in the DCA algorithm reduces the average number of reshuffles by 65.3% on Set I and by 30.83% on Set II. According to these results, we select the algorithm by Ding and Chou [4] and its randomization to be part of the GRASP.

6.3.2. Performance of the Local Search

In order to study the effect of the local search described in Section 5.2 by itself, before including it in the GRASP algorithm, we applied it to the solutions obtained using the DCA algorithm. Once the new configuration for port 1 has been obtained, the solution is completed using the same algorithm. Table 3 shows the results obtained on both sets. Besides the original local search, denoted in the table as LS I, and in order to reduce its running time, we also include a modified version, LS II, in which, given a group in column , instead of considering the exchanges of any number of containers we only consider the exchange of or   -containers.

Two points are noticeable in this table. First, the local search is able to find solutions for all the groups with a significantly lower average number of reshuffles, from 4.9 to 3.31 on Set I, and from 12 to 4.83 on Set II. Second, with the modified version of the local search, the average number of reshuffles is reduced by 24.8% on Set I and slightly increased by 2.1% on Set II. Nevertheless, the running times are reduced by 10.5% and 71.7% on Sets I and II, respectively. We therefore select this LS II version for the final configuration of the GRASP algorithm.

6.4. Results of the Complete GRASP Algorithm

Table 4 presents the aggregated results of GRASP for each group in Set I and Set II and summarizes the number of times that GRASP achieves a solution of zero reshuffles and the number of times that it achieves optimal solutions, either because they match an optimal solution provided by one or more of the IP models or because the solution has zero reshuffles. Differences between columns #Z and #Opt. only appear in Set II, on which the integer models have been tested.

With respect to Set I, GRASP is able to solve the majority of these instances, 603 out of 950, to proven optimality. For the remaining instances, the number of reshuffles in the solutions is greater than zero. Although their optimality cannot be proven, it is possible that some of them are in fact optimal, further increasing the total number of optimal solutions found by GRASP.

The instances from 4 to 14 ports are solved with an average number of reshuffles lower than or equal to 1 in an average time of 49.47 seconds. The instances with 16 and 18 ports are solved with an average lower than 2 and those of 20 and 22 ports with averages of 4.1 and 4.3 reshuffles. The average number of reshuffles in Set I obtained with GRASP is just 1.3, while that obtained with the DCA algorithm, the best according to Section 6.3.1, is 4.9. The results obtained with GRASP are also better than that obtained with the randomized version of DCA, 1.3 compared to 1.7 reshuffles, and with the local search technique LSII, 1.3 compared to 2.49. As expected, the stopping criterion that prevails when the number of columns and ports is so high is the time limit of 300 seconds.

The situation with Set II is somewhat different, because this set is composed of instances with larger numbers of reshuffles. The average number of reshuffles is more than twice that of Set I and the percentage of solutions with zero reshuffles decreases from 63.47% to 40.35%. While the average number of reshuffles for this set with the DCA algorithm is 12, and 8.3 with its randomized version, the average with GRASP is 4.2 reshuffles. With up to 10 ports, the new model provides slightly better solutions than GRASP, but for instances with 12 ports, GRASP outperforms the model with an average of 5 reshuffles instead of 7.4. The difference increases as the number of ports increases; in instances with 16 ports GRASP obtains an average of 11.6 reshuffles while for the model it is 68. Although the time limit was set at 300 seconds, the average time for all groups of instances in Set II is much shorter. In fact, GRASP finishes after 300 seconds only in 6 instances of 8 ports, 45 of 10 ports, 51 of 12 ports, 55 of 14 ports, and 83 of 16 ports. In most cases, the algorithm stops after 1000 iterations.

The results show the effectiveness of the GRASP we propose in this study. It obtains extremely good solutions in very short computing times, less than 5 minutes for the larger instances. These solutions are clearly better than those provided by other heuristic methods, which could be used to obtain acceptable solutions if extremely fast answers were needed. Integer models can only be used for small instances, and even in the cases in which the models can prove optimality, the solutions obtained by GRASP match or are very close to the optimal solutions. For medium and large-size instances, the integer models fail to obtain feasible solutions in many cases, and when they do obtain feasible solutions these are much worse than those provided by GRASP.

7. Conclusions and Future Work

Efficient stowage plans for increasingly larger container ships have a critical impact on the performance of operations at container terminals. The problem is being studied from various different perspectives, using a range of decomposition strategies and algorithmic approaches. In this article we have chosen to study the complete multiport problem, considering loading and unloading operations along the route, and have studied it from two complementary points of view. First, we have developed a new integer model and compared it with existing proposals in an extensive computational study. The conclusion is that, though our model performs consistently better than previous models, none of them is able to produce acceptable solutions for large-size instances. Therefore, in a second part we have developed a GRASP algorithm with new features, such as a randomizing strategy that randomizes the steps of the existing constructive algorithms. A new local search, tailored to the problem, improves the solutions of the randomized constructive phase and allows the GRASP to obtain many optimal solutions and extremely good solutions for instances for which optimality cannot be proven, clearly outperforming previous heuristic algorithms. Our conclusion is that the GRASP framework is very well suited to this problem because using a constructive algorithm is a good way of studying the sequence of operations in which the empty spaces left when unloading containers are then used in the loading operations. By adding randomization and a local search phase we achieve a very efficient algorithm for this complex problem.

The version of the problem considered in this study can be extended in several ways, progressively adding characteristics of the real problem. The two most important characteristics are the container type and weight. Considering the joint stowage of 20-foot and 40-foot containers introduces new constraints in the way they can be stacked. Reefer containers have only a limited set of available slots, those with electric plugs, high-cube containers reduce the number of tiers in a stack, and containers with dangerous goods bring new strict stowage conditions. In real problems these container types have to be considered and their characteristics taken into account by proposed solution procedures. Container weight is also of utmost importance, because it directly affects stability conditions that have to be respected for safety reasons. Although these features can be included in the integer linear model, the results obtained clearly indicate that the extended model will not be able to solve real-world instances. Extending the heuristic and metaheuristic algorithms seems a more promising approach, because their flexibility and modest computing time requirements make them able to cope with more complex versions of the problem.

Data Availability

The datasets Set I and Set II used to support the findings of this study are available from the corresponding author ([email protected]) upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This study has been partially supported by the Spanish Ministry of Education, Culture, and Sport, FPU Grant A-2015-12849, by the Spanish Ministry of Economy and Competitiveness, DPI2014-53665-P, cofinanced by FEDER funds, and by the Generalitat Valenciana, PROMETEO/2013/049, Grant CPI-13-351.