## Berth Allocation Problem with Rectangle Packing

• Post author:
• Post category:Blog

## Berth Allocation Problem (BAP)

The berth allocation problem assigns arriving ships to a section of the berth (where, when & how long) during enough time for the loading and unloading operations to take place.

Photo by Alex Duffy on Unsplash

We make the following assumptions and simplifications

• The berth is divided into sections of constant length, and a number of adjacent sections are simultaneously assigned to each berthed ship.
• Time is discrete and we assume a large enough planning horizon
• An average speed of container unloading is used to compute time needed to unload all the containers in a ship
• We ignore the number of cranes used for unloading containers, their location and their movements.

With these assumptions, the problem becomes a rectangle packing problem, which we will solve with IBM Cplex OPL Studio.

## Rectangle Packing Problem

The rectangle packing is a packing problem where rectangle have to be packed in a larger rectangle without any overlap.

In the berth allocation problem, each ship is be represented by a rectangle where the length is its unloading time and its height is number of segments occupied in the berth.

Ship Unloading time Berthing length
1 3 3
2 5 2
3 5 4
4 4 4
5 4 6

Various objective functions are possible, we will minimize the total time taken to complete all the unloading tasks

### Mathematical formulation

After assigning to each rectangle i a position defined by two floating point variables x_i and y_i, the next problem is to represent the non-overlap constraint. This is done by defining boolean variables that define for each pair (i, j) of rectangles whether i is on the left, right, above or below j. The link between the (x_i, y_i), (x_j, y_j) positions and the \text{left}_{ij}, \text{above}_{ij} variables is created using indicator constraints.

Indicator constraints are logic combinations of linear inequalities, which are automatically linearized by Cplex, but may also be linearized by hand if desired.

\begin{aligned}
& \quad \quad \min \max_i (x_i + W_i) & (1)\\
\\
& \forall ij \quad (l_{ij} = 1) \Leftrightarrow (x_i + W_i \leq x_j) & \quad (2)\\
& \forall ij \quad (a_{ij} = 1) \Leftrightarrow  (y_i + H_i \leq y_j) & \quad (3)\\
& \forall i < j \quad l_{ij} + l_{ji} + a_{ji} + a_{ji} \geq 1 & \quad (4)\\
& \forall i \quad x_i + W_i \leq W & \quad (5)\\
& \forall i \quad y_i + H_i \leq H & \quad (6)\\
\\
& \quad x_i\ y_i \in \mathbb{R}^+ &\\
& \quad l_i\ a_i \in \lbrace 0, 1  \rbrace &\\
\end{aligned}

The objective (1) minimizes the makespan of the assignment (the end of the last task). Constraints (2) and (3) define the boolean variable l_{ij} and a_{ij} as the indicator constraint that measures whether x_i + W_i \leq y_i. Equation (4) is the no-overlapping constraint and forces any par of rectangles to be either one on the left of the other, or one above the other. Constraints (5) and (6) ensure the rectangle remains within the admissible area.

### Linearizaton of the indicator constraint

The indicator constraint may be linearized in various ways

\begin{aligned}
& x_i + W_i + W * l_{ij} \leq x_j + W\\
& x_i + l_{ij} * (W_i + W) \leq x_j + W\\
& x_i + W_i \leq x_j + (1 - l_{ij}) * W\\
\end{aligned}

In all cases when l_{ij} = 1 the inequality reduces to x_i + W_i \leq x_j and when l_{ij} = 0 the equation reduces to something that is always true.

However, the model running time may differ depending on the chosen linearization if any or the usage of an indicator constraint.

### Improvements

Some improvements are easy to integrate in the model

#### Arrival times and waiting time

Arrival times just introduce extra bounds on the position variable x_i \geq R_i, which in turn allows a change in the objective to minimization of the maximum waiting time \min \max_i (x_i - R_i) or a minimization of the tardiness \min \sum_i (x_i + W_i - R_i)

#### Incompatibilities with berth positions

An incompatibility of a ship with some positions in the berth (because of insufficient depth for instance) creates a forbidden range of values for the y_i variable which can be represented by the constraint (y_i \leq L)\ ||\ (y_i \geq U). It is again a logic combination of linear constraints that can either be given directly to Cplex or linearized with the introduction of two boolean variables belowL and aboveU and the constraint b_L + a_U = 1

#### Lexicographic objective

Minimizing the makespan may result in solutions that are not completely "pushed to the left" because moving the yellow ship to an earlier time doesn’t change the end of the green ship.

A possible solution is to optimize lexicographically the makespan and the tardiness

dexpr float makespan = max (s in ships) (x[s] + ship_unloading_time[s]);
dexpr float tardiness = sum (s in ships) (x[s] + ship_unloading_time[s]);

minimize staticLex (makespan, tardiness);

constraints {
...
}

Which results in the following solution (results may vary)

#### Delays and stochastic optimization

Arrivals of ships as well as average loading times are estimations that may be wrong. Errors on average values are typically estimated with a probability distribution fitted to historical data – the most popular one being the normal distribution with parameters average \mu and deviation \sigma.

Stochastic optimization takes into account the fact that constants may deviate from their estimated average, and instead of optimizing over a single scenario, they optimize on expected value over a distribution of scenarios.

#### Cranes

Approximating the time it takes to unload the containers of a ship with an average is a very crude approximation as the number of cranes available on the berth is small, their mobility is limited, and the containers in the ship are located on multiple bays (positions along the berth).

Photo by Freddie on Unsplash

This justifies the interest in the combined berth assignment and crane scheduling problem which we will analyze in a subsequent blog post.

## An OPL model

For simplicity we have chosen the OPL model with indicator variables. The reader may modify it to try any of the linearizations suggested – which in our tests reduce the running times of the model.

The reader is invited to add constraints to the model and experiment with it.

### Model

/*********************************************
* OPL 12.10.0.0 Model
* Author: diegoolivierfernandezpons
* Creation Date: Nov 11, 2020 at 2:58:22 PM
*********************************************/

int number_of_ships = ...;

range ships = 1 .. number_of_ships;
int ship_length [ships] = ...;
int ship_unloading_time [ships] = ...;

int berth_length = ...;
int horizon = sum (s in ships) ship_unloading_time[s];

dvar float+ x[ships];
dvar float+ y[ships];

dexpr int left [i in ships][j in ships] = x[i] + ship_unloading_time[i] <= x[j];
dexpr int above [i in ships][j in ships] = y[i] + ship_length[i] <= y[j];

minimize max (s in ships) (x[s] + ship_unloading_time[s]);

constraints {

// Non overlapping constraint
forall(i, j in ships : i < j)
left[i][j] + left[j][i] + above[i][j] + above[j][i] >= 1;

// Ships within bounds
forall (i in ships) x[i] + ship_unloading_time[i] <= horizon;
forall (i in ships) y[i] + ship_length[i] <= berth_length - 1;
}


### Data

berth_length = 10;
number_of_ships = 5;

ship_length = [3,2,4,4,6];
ship_unloading_time = [3,5,5,4,4];

### Engine log

Version identifier: 12.10.0.0 | 2019-11-26 | 843d4de
CPXPARAM_TimeLimit                               300
Tried aggregator 2 times.
MIP Presolve eliminated 16 rows and 5 columns.
MIP Presolve modified 18 coefficients.
Aggregator did 23 substitutions.
Reduced MIP has 70 rows, 75 columns, and 182 nonzeros.
Reduced MIP has 41 binaries, 0 generals, 0 SOSs, and 25 indicators.
Presolve time = 0.00 sec. (0.17 ticks)
Probing changed sense of 2 constraints.
Probing time = 0.00 sec. (0.07 ticks)
Cover probing fixed 0 vars, tightened 6 bounds.
Tried aggregator 2 times.
MIP Presolve eliminated 2 rows and 2 columns.
Aggregator did 8 substitutions.
Reduced MIP has 60 rows, 65 columns, and 170 nonzeros.
Reduced MIP has 39 binaries, 0 generals, 0 SOSs, and 25 indicators.
Presolve time = 0.00 sec. (0.15 ticks)
Probing fixed 0 vars, tightened 2 bounds.
Probing time = 0.00 sec. (0.06 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 60 rows, 65 columns, and 170 nonzeros.
Reduced MIP has 39 binaries, 0 generals, 0 SOSs, and 25 indicators.
Presolve time = 0.00 sec. (0.11 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Clique table members: 30.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 4 threads.
Root relaxation solution time = 0.00 sec. (0.11 ticks)

Nodes                                         Cuts/
Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                           21.0000        9.0000            57.14%
0     0        9.0000    23       21.0000        9.0000       16   57.14%
0     0        9.0000    23       21.0000       Cuts: 8       22   57.14%
0     0        9.0000    23       21.0000       Cuts: 7       29   57.14%
0     0        9.0000    23       21.0000      Cuts: 22       41   57.14%
*     0+    0                           14.0000        9.0000            35.71%
0     0        9.0000    23       14.0000      Cuts: 12       50   35.71%
*     0+    0                           13.0000        9.0000            30.77%
*     0+    0                           12.0000        9.0000            25.00%
Detecting symmetries...
0     2        9.0000     6       12.0000        9.0000       50   25.00%
Elapsed time = 0.06 sec. (3.57 ticks, tree = 0.02 MB, solutions = 4)
*    10+    2                           11.0000        9.0000            18.18%

Clique cuts applied:  2
Cover cuts applied:  4
Implied bound cuts applied:  43
Mixed integer rounding cuts applied:  8
Gomory fractional cuts applied:  4

Root node processing (before b&c):
Real time             =    0.05 sec. (3.53 ticks)
Parallel b&c, 4 threads:
Real time             =    0.14 sec. (2.86 ticks)
Sync time (average)   =    0.09 sec.
Wait time (average)   =    0.00 sec.
------------
Total (root+branch&cut) =    0.20 sec. (6.39 ticks)

### Result

// solution (optimal) with objective 11
// Quality Incumbent solution:
// MILP objective                                 1.1000000000e+01
// MILP solution norm |x| (Total, Max)            1.04000e+02  1.10000e+01
// MILP solution error (Ax=b) (Total, Max)        0.00000e+00  0.00000e+00
// MILP x bound error (Total, Max)                0.00000e+00  0.00000e+00
// MILP x integrality error (Total, Max)          0.00000e+00  0.00000e+00
// MILP slack bound error (Total, Max)            0.00000e+00  0.00000e+00
// MILP indicator slack bound error (Total, Max)  0.00000e+00  0.00000e+00
//

x = [4 0 6 7 0];
y = [4 7 0 4 1];


## References

[1] Ak, Aykagan. 2011. Berth and Quay Crane Scheduling: Problems, Models and Solution Methods (Book and PhD Thesis, Georgia Institute of Technology, Dec 2008).