A Fuzzy Optimization Model for the Berth Allocation Problem and Quay Crane Allocation Problem (BAP + QCAP) with n Quays

Abstract

This work introduces a fuzzy optimization model, which solves in an integrated way the berth allocation problem (BAP) and the quay crane allocation problem (QCAP). The problem is solved for multiple quays, considering vessels’ imprecise arrival times. The model optimizes the use of the quays. The BAP + QCAP, is a NP-hard (Non-deterministic polynomial-time hardness) combinatorial optimization problem, where the decision to assign available quays for each vessel adds more complexity. The imprecise vessel arrival times and the decision variables—berth and departure times—are represented by triangular fuzzy numbers. The model obtains a robust berthing plan that supports early and late arrivals and also assigns cranes to each berth vessel. The model was implemented in the CPLEX solver (IBM ILOG CPLEX Optimization Studio); obtaining in a short time an optimal solution for very small instances. For medium instances, an undefined behavior was found, where a solution (optimal or not) may be found. For large instances, no solutions were found during the assigned processing time (60 min). Although the model was applied for n = 2 quays, it can be adapted to “n” quays. For medium and large instances, the model must be solved with metaheuristics.

1. Introduction

Maritime container terminals (MCTs) are vital elements in the global supply chain. The essential objective of an MCT is to provide the resources and organization to the transport of containers between the landside and maritime mediums. In this process, the MCT must guarantee the best conditions of speed, efficiency, and safety, in accordance with the environment [1].

Due to the current rise in the global maritime trade, many ports have suffered from resource constraints, such as space, time, quay cranes, etc. The problems that exist in an MCT are different. This work addresses the berth allocation problem (BAP) and the quay crane allocation problem (QCAP). The BAP is a NP-hard (Non-deterministic polynomial-time hardness) combinatorial optimization problem [2], which involves assigning each incoming vessel a berth position and arrival time at the quay. On the other hand, the QCAP tries to assign a number of quay cranes (QC) to each berth vessel, aiming to perform all the necessary unload or load movements of the containers in the vehicles. The QC are giant cranes that are mounted on rails, therefore, they cannot break through to each other.

The vessel arrival times are very uncertain, i.e., they can arrive earlier or later than their expected arrival time. This imprecision depends on several factors, such as: technical problems, weather conditions (winds, storms), additional terminal visits, or other reasons [3,4]. This has effects over other MCT activities, such as load and unload operations, negatively impacting the services required by customers. This work assumes the imprecise arrival times of the vessels, where the arrival is known but the exact time is uncertain, as this kind of uncertainty (ambiguity) can be modeled with fuzzy logic [5].

In [6], an exhaustive literature review which show a very limited number of reports, the authors deal with the imprecise times in the BAP + QCAP.

The authors of [7] use discrete event simulation to show that the collaboration between two MCTs (sharing resources such as berths, quay cranes, and the container yard) can help reduce uncertainty in arrival time and the number of containers brought in by vessel, although the authors do not present any mathematical models.

A reactive strategy for the BAP + QCAP with discrete berths is proposed in [8], a mixed integer linear programming (MILP) model with practical constraints, is formulated to obtain a basic planning when uncertainty problems appear (deviation of vessels’ arrival times, deviation of vessels’ loading and unloading operation times, unscheduled vessel calls, quay breakdowns, etc.), and a moving horizon heuristic is used to obtain good feasible solutions. Another reactive strategy for the BAP + QCAP is proposed in [9], where basic planning is obtained through the use of a multi-objective optimization model that penalizes: the deviations of the vessels from their preferred berthing positions, the delay in the berthing time in comparison to estimated arrival times, and the delay in departure times compared to estimated departure times; then, with movable time windows the berthing plans are updated.

Given that the model we propose is proactive, a review of works on this approach is made below.

Regarding discrete and dynamic BAP, where a single berth quay is available with some kind of imprecision present, a fuzzy mixed integer linear programming (MILP) model was proposed in [10], the imprecise arrival times were represented by triangular fuzzy numbers; however, this model does not addressed continuous BAP. According to [6], the design of a continuous model would have a more complicated berthing plan than a discrete one. Nevertheless, a continuous model has the advantage of better using the quay space. In [11], a fully fuzzy linear programing (FFLP) model was developed for dynamic and continuous BAP. The model obtained a robust berthing plan, which supports incidents in the vessel arrivals. In [12], a MILP model and a genetic algorithm (GA) were created for dynamic and continuous BAP+QCAP. In the modeling methodology, spaces of times were added in the vessel departure, which mitigates the effects of imprecision and strengthens the model’s accuracy.

In the case of BAPs, where vessels can berth in two quays; [13] proposed a MILP model and a GA. However, imprecision was not included in the problem parameters. On the other hand, a FFLP model for BAP which addressed the imprecision in the vessel arrivals was presented in [14]. A MILP model for integrated laycan and berth allocation and quay crane assignment problem (LBACAP) that considers multiple quays is proposed in [15].

As far as we know, no models have been developed for the BAP + QCAP which considers “n” quays and the uncertainty in the arrival times.

In this work, we present a fuzzy optimization model for the BAP + QCAP with two quays, continuous and dynamic, which considers the vessels’ imprecise arrival times, meaning they can arrive early or late of an allowed time. The vessels’ imprecise arrival times are represented by triangular fuzzy numbers. The model optimizes the use of the quays. In order to analyze the behavior and efficiency, the model is applied to a small, medium, and large instances respectively.

2. Materials and Methods

2.1. Fuzzy Arithmetic

Definition 1. 

Let X be the discourse universe, and a fuzzy set$\tilde{\mathrm{A}}$in X is a set of ordered pairs.

$$\tilde{\mathrm{A}}=\left\{\left(\mathrm{x},{\mathsf{\mu }}_{\tilde{\mathrm{A}}}\left(\mathrm{x}\right)\right), \mathrm{x}\in \mathrm{X}\right\}$$

(1)

where${\mathsf{\mu }}_{\tilde{\mathrm{A}}}:\mathrm{X}\to \left[0,1\right]$is the membership function, which represents the belonging degree of x with the set$\tilde{\mathrm{A}}.$

In this work, fuzzy sets are defined over the real numbers $ℝ$. A membership function can be triangular, trapezoidal, sigmoidal, quadratic, etc.

Definition 2. 

A fuzzy singleton is a fuzzy set with a membership function.

$${\mathsf{\mu }}_{\tilde{\mathrm{A}}}\left(x\right)=\{\begin{array}{c}1, x \in \left[a-\epsilon ,a+\epsilon \right]\\ 0, x \notin \left[a-\epsilon ,a+\epsilon \right]\end{array}$$

(2)

where ε is a margin of error. The fuzzy singleton allows the expression of a real number as a fuzzy set.

Definition 3. 

A fuzzy number is a normal and convex fuzzy set in$ℝ.$

Definition 4. 

A triangular fuzzy number is represented by$\tilde{\mathrm{a}}=\left(\mathrm{a}1,\mathrm{a}2,\mathrm{a}3\right).$

Definition 5. 

A triangular fuzzy number$\tilde{a}=\left(\mathrm{a}1,\mathrm{a}2,\mathrm{a}3\right)$is positive if and only if$\mathrm{a}1>0$.

Definition 6. 

A real number$a$can be represented by a fuzzy singleton in the form of a triangular fuzzy number$\tilde{a}=\left(\mathrm{a}-\mathsf{\epsilon },\mathrm{a},\mathrm{a} + \mathsf{\epsilon }\right)$, where ε is a margin of error.

Definition 7. 

For two non-negative triangular fuzzy numbers$\tilde{\mathrm{a}}=\left(\mathrm{a}1,\mathrm{a}2,\mathrm{a}3\right)$and$\tilde{\mathrm{b}}=\left(\mathrm{b}1,\mathrm{b}2,\mathrm{b}3\right)$, the sum difference and multiplication operations are defined as follows:

$$\begin{gathered} \tilde{a}+\tilde{b}=\left(\mathrm{a}1+ \mathrm{b}1,\mathrm{a}2+ \mathrm{b}2,\mathrm{a}3+ \mathrm{b}3\right) \\ \tilde{a}-\tilde{b}=\left(\mathrm{a}1-\mathrm{b}3,\mathrm{a}2-\mathrm{b}2,\mathrm{a}3-\mathrm{b}1\right) \\ \tilde{a}.\tilde{b}=\left(\mathrm{a}1.\mathrm{b}1,\mathrm{a}2.\mathrm{b}2,\mathrm{a}3.\mathrm{b}3\right). \end{gathered}$$

(3)

Ordering methods allows us to decide the greater number between two fuzzy numbers $\tilde{a}$ and $\tilde{b}$. Fuzzy numbers do not always provide an ordered set as found in real numbers. Ordering methods in fuzzy numbers have advantages and disadvantages. There are different ordering methods, depending on the representation: preference, rationality, and robustness, etc. [16].

In this work, the objective of the fuzzy model optimization is to achieve an ordering (planning) for the berthing of vessels, therefore, any ordering method that orders fuzzy numbers could be used. Ordering methods that use intervals to order fuzzy numbers are not recommended, as they could add more imprecision to the model.

This work utilizes the Yagger First Index [17] ordering method, which is defined below.

Definition 8. 

Given two triangular fuzzy numbers$\tilde{a}$and$\tilde{b}$alongside the ordering function,

$$\Re \left(\tilde{a}\right)=\frac{\mathrm{a}1+ \mathrm{a}2+ \mathrm{a}3}{3}$$

(4)

$\tilde{a}\le \widetilde{b }$when, $\Re \left(\tilde{a}\right)\le \Re \left(\tilde{b}\right)$, Meaning,$\tilde{a}\le \widetilde{b }$when$\mathrm{a}1+ \mathrm{a}2+ \mathrm{a}3\le \mathrm{b}1+ \mathrm{b}2+ \mathrm{b}3$.

2.2. Fully Fuzzy Linear Programming (FFLP)

There are different approaches in fuzzy mathematical programming, the classical fuzzy linear programming methods are used when the parameters are fuzzy; in this work, some of the decision variables are fuzzy, for this case, the FFLP approach is the most appropriate.

In the FFLP approach, the decision parameters and variables are fuzzy and linear, respectively. Several solution methods can be applied to the FFLP model [18]. Most of them convert the fuzzy model into a classic linear programming one. In this work, we use the method proposed by Nasseri [19].

Given the FFLP problem:

$$max{\displaystyle \sum }_{j=1}^n{\tilde{C}}_j{\tilde{X}}_j$$

(5)

Subject to:

$${\displaystyle \sum }_{j=1}^n{\tilde{a}}_{ij}{\tilde{x}}_j\le {\tilde{b}}_i , \forall i=1\dots m$$

(6)

where the parameters ${\tilde{C}}_j , {\tilde{a}}_{ij} , {\tilde{b}}_j$ and the decision variable ${\tilde{x}}_j$ are non-negative triangular fuzzy numbers, ${\tilde{C}}_j=\left(c{1}_j,c{2}_j,c{3}_j\right)$, ${\tilde{a}}_{ij}=\left(a{1}_{ij},a{2}_{ij},a{3}_{ij}\right)$, ${\tilde{b}}_j=\left(b{1}_i,b{2}_i,b{3}_i\right)$ y ${\tilde{x}}_j=\left(x{1}_j,x{2}_j,x{3}_j\right)$.

The Nasseri method transforms the previous model into a classic linear programming problem,

$$\max \Re ({\displaystyle \sum }_{j=1}^n\left(c{1}_j,c{2}_j,c{3}_j) (x{1}_j,x{2}_j,x{3}_j\right))$$

(7)

Subject to:

$${\displaystyle \sum }_{j=1}^{n}a{1}_{ij}x{1}_j\le b{1}_i , \forall i=1\dots m$$

(8)

$${\displaystyle \sum }_{j=1}^{n}a{2}_{ij}x{2}_j\le b{2}_i , \forall i=1\dots m$$

(9)

$${\displaystyle \sum }_{j=1}^{n}a{3}_{ij}x{3}_j\le b{3}_i , \forall i=1\dots m$$

(10)

where $\Re$ is an ordering function (see Definition 8).

2.3. Problem Description

According to [20], an MCT is a composition of several subsystems integrated into a single one. This system has physical and information connections with the transport networks (landside and maritime). The subsystems (see Figure 1) are:

• The loading–unloading container subsystem, which is responsible for resolving the maritime interface.
• The storage container subsystem, which occupies most of the MCT surface.
• The landside reception and delivery subsystem, which act as gates in the landside for trucks and/or railways.
• The internal connection subsystem. To the previous three subsystems, which address the basic terminal functions, a fourth subsystem must be added, this ensures the horizontal transport of containers between the previous subsystems.

The BAP + QCAP discussed in this work occur in subsystem 1, however any incident will affect the other subsystems. For example, if a vessel berths outside their scheduled time, it will cause problems to the assigned cranes, affecting the prepared warehouses for the containers and the trucks waiting to receive them.

The existence of n quays in the port is assumed, where vessels can berth in any of them. Arrival time is imprecise, i.e., vessels are allowed to be early or late until a predetermined time. The BAP focuses on choosing a quay (if any) and the arrival time and position where each arriving vessel at the MCT must berth. The other problem is the QCAP, which involves the assigning of a number of quay cranes to each vessel to be handled.

The BAP can be represented in a two-dimensional form (see Figure 2), where the horizontal axis represents the time and the vertical axis represents the berth length.

We consider the following assumptions and limitations:

Assumptions: The vessel information is previously known, each vessel has a draft less or equal than the quay, and that the berthing and departure do not consume much time. Simultaneous berthing is allowed and the safety distance between vessels is not considered. The number of cranes assigned to the vessels does not change during the berthing time, once a QC starts a task on a vessel they must complete it without any pause, all QCs assigned to vessel i have the same working time ($t_{ik}=h_i, \forall k\in QC, u_{ik}=1$). The model assigns the berthing place to each vessel, that is, preference zones are not allowed for berthing a vessel.

Limitations: The quay length must be limited, the available crane number in each quay is five. There is a safety distance between cranes which must to be maintained (35 m). The maximum crane number that can be assigned to a vessel is four. The number of crane movements performed in a given time is 2.5. There must be at least one QC assigned to each vessel.

In this work we use the following notation, which was taken from [12]:

• $H$: Planning Horizon.

Let be $Q$ the set of existent quays in an MCT, where $q\in Q$ is a quay.

• $Q{C}_q$: Available cranes at quay q. All $QCs$ perform the same number of movements per unit time $\left(movsQC\right)$, data provided by the MCT.
• $Q{C}_{iq}^{ + }$: The maximum number of QCs assigned to each vessel i at quay q
• $L_q$: Length of quay $q$.

Let be $V$ a set of vessels which arrive at the port, the data for each vessel $i\in V$ is given by:

• $a_i$: Vessel arrival time at port.
• $l_i$: Vessel length.
• $c_i$: Required number of movements to load or unload the containers from the vessel.

The decision variables are:

• $B{M}_{iq}$: Takes the value of 1 when vessel i berths at quay $q$, and takes 0 otherwise.
• $m_{iq}$: Berthing time of vessel i at quay $q$. The waiting time $(w_i)$ is computed as $(w_i=m_{iq}-a_i)$
• $p_{iq}$: Position at quay $q$, where vessel i must berth.
• $n_{iq}$: Number of $QCs$ in quay $q$ assigned to vessel i.
• $u_{ik}$: Indicates whether k, (which belonging to $Q{C}_q$) works (1) or not (0) on vessel i.

The variables that are deduced from above are:

• $h_i$: Vessel handling time $i$. $h_i=c_i/(n_{iq}\ast movsQC).$
• $t_{ik}$: Working time of k (belonging to $Q{C}_q$) which is assigned to vessel $i$.
• $d_i$: Vessel departure time $\left(d_i= m_{iq} + h_i\right).$
• $s_{iq} ; e_{iq}$: Indices for the first and last $QC$, on quay $q$, used in vessel $i$, respectively.
• $M$: Is a sufficiently large number.

2.4. Fuzzy Optimization Model for the BAP + QCAP with Two Quays

The arrival, berthing, handling, waiting and departure time of each vessel are considered imprecise (fuzzy), they are denoted by: $\tilde{a}$, $\tilde{m}$, $\tilde{h}$, $\tilde{w}$ and $\tilde{d}$ respectively.

Assuming the presence of imprecision in some parameters and decision variables, the fuzzy optimization model for BAP + QCAP is introduced. This model is based on the deterministic model developed in [12]. The objective function minimizes the waiting and handling time for all vessels.

$$min{\displaystyle \sum }_{i\in V}\left({\tilde{w}}_i + {\tilde{h}}_i\right)$$

(11)

Subject to:

$${\displaystyle \sum }_{q\in Q}B{M}_{iq}=1 , \forall i\in V ,\forall q\in Q$$

(12)

$${\tilde{m}}_{iq}\ge {\tilde{a}}_i , \forall i\in V, \forall q\in Q$$

(13)

$${\tilde{w}}_i + {\tilde{a}}_i ={\tilde{m}}_{iq} , \forall i\in V, \forall q\in Q$$

(14)

$${\tilde{d}}_i={\tilde{m}}_{iq} + {\tilde{h}}_i , \forall i\in V, \forall q\in Q$$

(15)

$$p_{iq} + l_i \le L_q , \forall i\in V, \forall q\in Q$$

(16)

$$n_{iq}={\displaystyle \sum }_{k\in Q{C}_q}{u}_{ik} , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(17)

$$1\le n_{iq} \le Q{C}_{iq}^{ + } , \forall i\in V, \forall q\in Q$$

(18)

$$1\le s_{iq} , e_{iq} \le \left|Q{C}_q\right| , \forall i\in V, \forall q\in Q$$

(19)

$$s_{iq}\ge e_{iq} , \forall i\in V, \forall q\in Q$$

(20)

$$n_{iq}=e_{iq}- s_{iq} + 1 , \forall i\in V, \forall q\in Q$$

(21)

$${\displaystyle \sum }_{k\in Q{C}_q}{t}_{ink}\ast movQC\ge c_i , \forall i\in V, \forall q\in Q$$

(22)

$${\tilde{h}}_i= ma{x}_{k\in Q{C}_q}{t}_{ik} , \forall i\in V, \forall q\in Q$$

(23)

$$t_{ik}-M\ast u_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(24)

$${\tilde{h}}_i-M\left(1- u_{ik}\right)-t_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(25)

$$u_{ik} + u_{jk} + z_{ijq}^x<2 , \forall i,j\in V , i\ne j, \forall k\in Q{C}_q, \forall q\in Q$$

(26)

$$M\left(1-u_{ik}\right) + \left(e_{iq}-k\right)\ge 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(27)

$$M\left(1-u_{ik}\right) + \left(k-s_{iq}\right)\ge 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(28)

$$p_{iq} + l_i \le p_{jq} + M\left(1-z_{ij}^x\right) , \forall i,j\in V, i\ne j, \forall q\in Q$$

(29)

$$e_{iq} + 1\le s_{jq} + M\left(1- z_{ijq}^x\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(30)

$${\tilde{d}}_i\le {\tilde{m}}_{jq} + M\left(1-z_{ijq}^y\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(31)

$${\tilde{m}}_{iq} + {\tilde{h}}_i \le H , \forall i\in V, \forall q\in Q$$

(32)

$$z_{ijq}^x + z_{jiq}^x + z_{ijq}^y + z_{jiq}^y\ge 1 , \forall i,j\in V, i\ne j, \forall q\in Q$$

(33)

$$z_{ijq}^x, z_{ijq}^y\in \left\{0,1\right\} , \forall i,j\in V, i\ne j, \forall q\in Q$$

(34)

There are two auxiliary variables: $z_{ijq}^x$ is a decision variable that indicates whether the vessel $i$ is located to the left of vessel $j$ at quay $(z_{ijq}^x=1)$, while $z_{ijq}^y$ indicates that vessel $i$ berths before vessel $j$ in time, at quay $q$ $(z_{ijq}^x=1)$ (see restriction 34).

Constraint (12), assigns each vessel to a quay $q$. Constrain (13), indicates that all vessels can berth once they arrive at the port. Constraints (14) and (15) set the waiting and departure times for vessel $i$, according with the berthing time. Constraint (16) ensures the berth position of vessel $i$ does not exceed the length of the quay $q$.

Constraints (17)–(21) assign a number of working $QCs$ at quay $q$ for vessel $i$. Constraint (22) sets the minimum handling time required to load or unload containers, according with the assigned number of $QCs$ to vessel $i$. Constraint (23) assigns the handling time to vessel $i$. Constraint (24) ensures that unassigned $QCs$ have a value of $t_{ik}=0$. Constraint (25) forces the assigned $QCs$ in vessel $i$ to work the same number of hours. Constraint (26) prevents a $QC$ from being assigned to different vessels at the same time. Constraints (27) and (28) force $QCs$ to be contiguously assigned (from $s_i$ to $e_i$) at quay q. Constraint (29) takes into account the safety distance. Constraint (30) prevents a vessel from using a $QC$ which breaks through other $QCs$. Constraint (31) prevents vessel $j$ from berthing at the quay, while vessel $i$ is still berthed at the same quay. Constraint (32) indicates that the berthing and handling time for vessel $i$ must not exceed the planning horizon (H). Finally, constraint (33) establishes the relationship between each pair of vessels.

2.5. Fuzzy Optimization Model Solution

It is assumed that, in the fuzzy optimization model from Section 2.4, the imprecision times related with the vessel operations: arrival, berthing, handling, and departure, are represented by the triangular fuzzy numbers $\tilde{a}=\left(a1,a2,a3\right)$, $\tilde{m}=\left(m1, m2, m3\right)$, $\tilde{h}=\left(h1,h2,h3\right),$ and $\tilde{d}=\left(d1,d2,d3\right)$, respectively. $\tilde{h}$, is considered a fuzzy singleton; also, if the model parameters are linear, the model is a FFLP.

The method used to transform the FFLP model into a classic linear programming model requires the application of fuzzy arithmetic operations; restrictions (25), (31), and (32) show such operations between fuzzy and real numbers, in order to perform such operations, real numbers are considered to be fuzzy singletons (See Definition 6), for example, in constraint (31),

$${\tilde{d}}_i\le {\tilde{m}}_{jq} + M\left(1-z_{ijq}^y\right)$$

Is transformed into,

$$\left(d{1}_i,d{2}_i,d{3}_i\right)\le \left(m{1}_{jq}, m{2}_{jq},m{3}_{jq}\right) + M\left(\left(1,1,1\right)-\left(z_{ijq}^y, z_{ijq}^y, z_{ijq}^y\right)\right)$$

Simplifying,

$$\left(d{1}_i,d{2}_i,d{3}_i\right)\le \left(m{1}_{jq} + M\left(1-z_{ijq}^y \right), m{2}_{jq} + \left(1-z_{ijq}^y \right), m{3}_{jq} + \left(1-z_{ijq}^y \right)\right)$$

As indicated in the Nasseri method (see Section 2.2), operations between fuzzy numbers which involve sum, difference, and multiplication operations, are performed in the FFLP model; an ordering method is applied to the objective function, in this case, the First Yagger Index is used (see Definition 8), obtaining the following MILP model.

$$min{\displaystyle \sum }_{i\in V}\frac{1}{3}\left((m{1}_{iq}-a{3}_i + h{1}_i\right) + \left(m{2}_{iq}-a{2}_i + h{2}_i\right) + (m{3}_{iq}-a{1}_i + h{3}_i))$$

(35)

Subject to:

$${{\displaystyle \sum }}_{q\in Q}B{M}_{iq}=1 , \forall i\in V ,\forall q\in Q$$

(36)

$$m{1}_{iq}\ge a{1}_i , m{2}_{iq}\ge a{2}_i ,m{3}_{iq}\ge a{3}_i , \forall i\in V, \forall q\in Q$$

(37)

$$m{3}_{iq}\ge m{2}_{iq} , m{2}_{iq}\ge m{1}_{iq} , \forall i\in V, \forall q\in Q$$

(38)

$$\begin{gathered} w{1}_i + a{1}_i =m{1}_{iq}, w{2}_i + a{2}_i=m{2}_{iq} , \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}w{3}_i + a{3}_i =m{3}_{iq} , \forall i\in V, \forall q\in Q \end{gathered}$$

(39)

$$d{1}_i=m{1}_{iq} + h{1}_i , d{2}_i=m{2}_{iq} + h{2}_i , d{3}_i=m{3}_{iq} + h{3}_i , \forall i\in V, \forall q\in Q$$

(40)

$$p_{iq} + l_i \le L_q , \forall i\in V, \forall q\in Q$$

(41)

$$n_{iq}={\displaystyle \sum }_{k\in Q{C}_q}{\tilde{u}}_{ik} , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(42)

$$1\le n_{iq} \le Q{C}_{iq}^{ + } , \forall i\in V, \forall q\in Q$$

(43)

$$1\le s_{iq} , e_{iq} \le \left|Q{C}_q\right| , \forall i\in V, \forall q\in Q$$

(44)

$$s_{iq}\ge e_{iq} , \forall i\in V, \forall q\in Q$$

(45)

$$n_{iq}={\tilde{e}}_{iq}- {\tilde{s}}_{iq} + 1 , \forall i\in V, \forall q\in Q$$

(46)

$${\displaystyle \sum }_{k\in Q{C}_q}{t}_{ink}\ast movQC\ge c_i , \forall i\in V, \forall q\in Q$$

(47)

$$\begin{gathered} h{1}_i= ma{x}_{k\in Q{C}_q}{t}_{ik} , h{2}_i= ma{x}_{k\in Q{C}_q}{t}_{ik} , \\ h{3}_i= ma{x}_{k\in Q{C}_q}{t}_{ik} , \forall i\in V, \forall q\in Q \end{gathered}$$

(48)

$$t_{ik}-M\ast u_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(49)

$$h{1}_i-M\left(1- u_{ik}\right)-t_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(50)

$$h{2}_i-M\left(1- u_{ik}\right)-t_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(51)

$$h{3}_i-M\left(1- u_{ik}\right)-t_{ik}\le 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(52)

$$u_{ik} + u_{jk} + z_{ijq}^x<2 , \forall i,j\in V , i\ne j, \forall k\in Q{C}_q, \forall q\in Q$$

(53)

$$M\left(1-u_{ik}\right) + \left(e_{iq}-k\right)\ge 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(54)

$$M\left(1-u_{ik}\right) + \left(k-s_{iq}\right)\ge 0 , \forall i\in V, \forall k\in Q{C}_q, \forall q\in Q$$

(55)

$$p_{iq} + l_i \le p_{jq} + M\left(1-z_{ij}^x\right) , \forall i,j\in V, i\ne j, \forall q\in Q$$

(56)

$$e_{iq} + 1\le s_{jq} + M\left(1- z_{ijq}^x\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(57)

$$d{1}_i\le m{1}_{jq} + M\left(1-z_{ijq}^y\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(58)

$$d{2}_i\le m{2}_{jq} + M\left(1-z_{ijq}^y\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(59)

$$d{3}_i\le m{3}_{jq} + M\left(1-z_{ijq}^y\right) , \forall i\in V, i\ne j, \forall q\in Q$$

(60)

$$m{1}_{iq} + h{1}_i \le H , m{2}_{iq} + h{2}_i \le H , m{3}_{iq} + h{3}_i \le H , \forall i\in V, \forall q\in Q$$

(61)

$$z_{ijq}^x + z_{jiq}^x + z_{ijq}^y + z_{jiq}^y\ge 1 , \forall i,j\in V, i\ne j, \forall q\in Q$$

(62)

$$z_{ijq}^x, z_{ijq}^y\in \left\{0,1\right\} , \forall i,j\in V, i\ne j, \forall q\in Q$$

(63)

3. Results

For both the case study and the model evaluation efficiency a set of uniform distributed instances were used. The model was implemented with the IBM ILOG CPLEX Optimization Studio (CPLEX) solver and was programmed on a personal computer: Intel(R) Core (TM) i7-8550U CPU, 1.80 GHz processor and 8 GB of RAM. The experiments were conducted within a 60 min timeout.

3.1. Study Case

To evaluate the model, an instance of ten vessels was used as a case study (see

Table 1. Instance of ten vessels with imprecise arrival times.

Vessel	a1	a2	a3	l (m)	QC Mov.
V1	14	16	20	260	4160
V2	18	31	48	232	9680
V3	56	68	86	139	3640
V4	81	82	97	193	7610
V5	92	105	119	287	6860
V6	106	116	133	318	6300
V7	126	138	147	366	8110
V8	155	167	176	166	1560
V9	159	163	164	109	7830
V10	162	179	186	251	2220

); the input data are the vessel imprecise arrival times (early, exact, late), its length, and the number of crane movements required to handle it. Figure 3 displays the vessels’ uncertain arrivals from Table 1, represented as triangular fuzzy numbers.

For example, for vessel V7, the arrival time with all its possibilities gives us 138 time units, but it can also arrive early or late by up to 126 and 147 time units, respectively. The vessel length is 366 m and 8110 crane movements are required to handling it.

The model assumes that: two quays are used; each having five cranes, with a minimum of one crane operating on each vessel to a maximum of four cranes per vessel. Each quay has a length of 700 m.

The obtained results are shown in

Table 2. The fuzzy berthing plan and crane allocation obtained by the model.

Vessel	m1	m 2	m 3	h	d1	d2	d3	p	Cranes	Q
V1	14	16	20	694	708	710	714	440	2	0
V2	927	928	943	1076	2003	2004	2019	0	3	1
V3	56	68	86	607	663	675	693	193	2	1
V4	81	82	97	846	927	928	943	0	3	1
V5	583	600	607	763	1346	1363	1370	0	3	0
V6	708	710	714	1050	1758	1760	1764	382	2	0
V7	663	675	693	1352	2015	2027	2045	334	2	1
V8	409	426	433	174	583	600	607	0	3	0
V9	1346	1363	1370	870	2216	2233	2240	0	3	0
V10	162	179	186	247	409	426	433	189	3	0

. Within an hour of computing time, an objective value of 15,570 was obtained, this is not the optimal value; but the best value obtained within that time.

Berthing time (m1, m2, m3) and output time (d1, d2, d3) are triangular fuzzy numbers. Q = 1 refers to quay 1, and Q = 0 to quay 2.

For example, vessel V7 can berth between the time units 663 and 693, with more possibility at time unit 667. It can depart between units 2015 and 2045, with more possibility at unit 2027. Additionally, it must berth at position 334 of the quay. Finally, two cranes are assigned to this vessel, and they must berth at quay 1.

The fuzzy berthing plan is displayed in Figure 4. Vessels are represented as polygons (not as rectangles as in the deterministic problem). The red polygon line represents the allowed early time which can be tolerated for the vessel to berth; meanwhile the green line is the tolerated late arrival time. The small triangle represents the berthing with more possibility.

In Figure 4, in quay 2, the blue circle suggests a conflict area between the departure and berthing of vessels V5 and V9. However as explained in Figure 5, such conflict is not real. For example, if vessel V5 departures late from the quay at time unit 1365, vessel V9 could berthing between the 1365 and 1370 units.

In order to verify the fuzzy model robustness, incidences were simulated in the vessel arrivals (see

Table 3. Arrival incidence in vessels.

Vessel	Incidence	Time
V1	Exactly	0
V2	Exactly	0
V3	Early	6
V4	Late	10
V5	Exactly	0
V6	Late	15
V7	Exactly	0
V8	Late	8
V9	Early	3
V10	Exactly	0

). With these incidences, the final berthing plan is obtained (see

Table 4. Final berthing plan.

Vessel	m	l	h	d	p	Cranes	Q
V1	16	260	694	710	440	2	0
V2	938	232	1076	2014	0	3	1
V3	62	139	607	669	193	2	1
V4	92	193	846	938	0	3	1
V5	583	287	763	1346	0	3	0
V6	710	318	1050	1760	382	2	0
V7	669	366	1352	2021	334	2	1
V8	409	166	174	583	0	3	0
V9	1346	109	870	2216	0	3	0
V10	179	251	247	426	189	3	0

). Figure 6 shows that the final berthing plan is part of the fuzzy berthing plan.

3.2. Model Efficiency Analysis

In order to analyze the model efficiency, data from vessels 5 to 35 with 100 instances each were used. The results are shown in

Table 5. Instances evaluation results.

Vessels	Average Objective Value	Average Processing Time (Minutes)	Optimal	No Optimal
5	4922.0	0.2	100	0
6	6655.3	1.3	47	53
7	9401.3	60.0	0	1
8	10,499.7	60.0	0	1
9	13,724.7	60.0	0	1
10	15,522.3	60.0	0	1
11	19,741.3	60.0	0	1
12	23,714.3	60.0	0	1
13	28,762.3	28.9	2	98
14	36,194.7	19.0	3	97
15	39,153.0	38.8	1	99
16	42,786.3	30.8	2	98
17	49,753.0	25.0	2	98
18	57,623.7	28.5	2	98
19	68,661.0	24.7	2	98
20	71,727.3	22.8	2	98
21	80,968.3	21.2	2	98
22	92,723.3	22.1	2	98
23	88,050.0	60.0	0	1
24	96,369.0	60.0	0	1
25	110,842.3	60.0	0	1
26	108,655.3	60.0	0	1
27	128,984.3	60.0	0	1
28	116,706.3	60.0	0	1
29	172,058.3	39.7	1	99
30	142,758.3	60.0	0	1
31	158,178.0	60.0	0	1
32	177,955.3	60.0	0	1
33	156,638.7	60.0	0	1
34	200,806.0	60.0	0	1
35	-	60.0	0	0

.

Optimal values were found for five vessels, with an average objective value and processing time of 4922 and 0.2 min, respectively; this was the unique number of vessels for which an optimal solution was obtained in all its instances. For six vessels, an average objective value of 6655.3 was obtained within an average processing time of 1.3 min, where a total of 47 instances were optimally solved. No optimal solutions were found for 7 to 12 vessels. Instead, a single non-optimal solution was found in just one instance. For 13 to 22 vessels there were between one and three optimal solutions. No optimal solutions were found for 23 to 34 vessels, rather, just a best solution in the given processing time. However, for 29 vessels an optimal solution was obtained. For 35 vessels onwards, no solution was found.

Regarding the processing time (see Figure 7), is noted that for six vessels, the average time it took to find a solution was 1.3 min. Meanwhile, for seven vessels the average processing time drastically increases until it matches the allowed processing time (60 min). For 13 to 22 vessels the time varies between 19 and 38.8 min, respectively. For 23 vessels onwards, the time processing was 60 min, except for 29 vessels, which takes 39.7 min.

Figure 8 shows a trend of increase in the target value, alongside, the polynomial curve which best adjusts the data; this is a quadratic curve, with a determination coefficient R² of 0.9775.

4. Discussion

A fuzzy optimization model for the BAP + QCAP with two quays was developed, considering the vessels’ imprecise arrivals. The model was implemented in the CPLEX solver. For previously known input data from a set of vessels, the model obtained a fuzzy berthing and crane allocation plan, which supported incidences of vessels being early or late in their arrival times.

The model efficiency was evaluated with a benchmark of 100 instances for each number of vessels, with one hour of computing processing. Only for a very small number of instances (five vessels) were optimal solutions obtained in all instances. For 6 to 34 vessels, only some optimal and non-optimal solutions were obtained: these do not follow a defined behavior regarding the processing time, this is because the problem is NP-hard and solutions will only be obtained in medium instances that have a structure that allows the algorithm used by the Solver CPLEX to find a solution. For 35 vessels onwards, the CPLEX solver was unable to find solutions within the allowed processing time. The same results applied for large instances.

The model was designed for “n” quays, but in this work is applied to only n = 2. Each time the quay number increases, the complexity will increase as well.

For medium and large instances, the model must be solved with metaheuristics, because it is a highly complex combinatorial optimization problem.