Energy-E ﬃ cient and Integrated Allocation of Berths, Quay Cranes, and Internal Trucks in Container Terminals

: Despite a signiﬁcant number of studies that have focused on the operational e ﬃ ciency of container terminals, existing literature has paid little attention to improving energy e ﬃ ciency, e.g., energy consumption and negative externalities in container terminals. Most researchers consider energy-saving goals when allocating berths and quay cranes to vessels, assuming that internal trucks are in su ﬃ cient supply. Furthermore, recent studies have revealed that shortage of internal trucks has become an issue that greatly a ﬀ ects the operational and energy e ﬃ ciencies of container terminals. This work presents a planning model that integrates berth allocation, quay crane assignment, and internal truck assignment problems. The developed model contributes to existing literature by including energy-saving goals in the integrated planning of these problems, as well as including important realistic factors such as shortages of internal trucks and handling time estimations, thus producing a reliable handling plan that achieves energy and cost savings without additional truck investment. To solve realistic problems, a Lagrangian relaxation-based method is developed. Furthermore, the beneﬁts of the developed approach are demonstrated by comparing it to an existing approach. On average, our approach could improve the solutions of the integrated problem with di ﬀ erent numbers of internal trucks by 6% compared to the solutions obtained using the existing approach.


Introduction
Today, maritime transport has become an essential element of supply chains of services and products: more than 80% of global freight transport is by sea. It is expected that global maritime trade will increase, with an average annual growth rate of 3.5%, within the 2019-2024 period due to the increase in containerized, dry bulk, and gas cargos [1]. This in turn will lead to more pressure on seaports to accommodate the growing demand and provide adequate services with limited resource capacities. At a container terminal, once the vessel is settled at the quay, quay cranes (QCs) and internal trucks (ITs) are used to handle the containers on the vessel. The QC loads or unloads the containers to or from the IT, which transports the container from or to the storage area, where the IT is loaded or unloaded by the yard crane. The IT will repeatedly make this round trip until all containers on the vessels are handled. Accordingly, the handling time of the vessel is greatly affected by its assigned numbers of QCs and ITs. Existing studies have used two different names for the manned trucks that move containers between the quay side and yard side: yard trucks [2] and internal trucks [3]. In this paper, we use "internal trucks" to describe manned container trucks. Berths, QCs, and ITs are expensive resources at container terminals [4]. Therefore, there are considerable studies that have focused on the planning decisions related to these resources [5,6]. The berth allocation problem (BAP) determines the berthing schedule for each vessel. Typically, the berthing schedule is developed based on the estimated values of vessels' handling times. The berthing schedule of a vessel describes its berthing position, berthing time, and completion time. The quay crane assignment problem (QCAP) calculates how many QCs are assigned to each vessel, while the internal truck assignment problem (ITAP) calculates how many ITs are assigned to the QCs that handle the containers on the vessel. The QCAP and ITAP together determine the handling times of the vessels, so that a specific objective is achieved.
With increasing numbers of containers, more resources are needed for container handling operations. However, purchasing additional resources is not an easy decision due to the high investment costs of these resources. Therefore, optimizing terminal operations has become necessary for container terminals to improve vessel service time, which is as an important measure of service quality [7]. The increasing demand for maritime transport has also led to more energy consumption and its related environmental impacts, such as greenhouse gas (GHG) emissions. For seaports, significant overhead costs are caused by energy consumption, and therefore, it is essential to reduce energy consumption to achieve more cost reduction and profits. At the same time, reducing energy consumption can support the transition toward greener and more sustainable container terminals, which has been of global concern [8,9]. Handling equipment and container vessels are often operated by diesel engines, which are significant sources of emissions and energy consumption in the port environment. The operational efficiency of container terminals (e.g., short service times) is highly dependent on effective planning of the handling operations of vessels using the available berths, QCs, and ITs [10]. Some studies have shown that improving the operational efficiency (minimizing vessel service time) might lead to energy savings. This is possible if there is a positive correlation between the two objectives, for example, vessels usually consume fuel while in idle mode for lighting purposes, refrigerating goods, and other onboard service purposes. A vessel is in idle mode when it is at port waiting for its service to start. Therefore, reducing the time that a vessel waits until its service starts can directly reduce the energy consumed during idle mode [11]. However, a negative correlation may exist between the two objectives. For example, vessels may arrive at a terminal late, and to reduce their waiting times, vessels may be allocated to available berths on a first-come, first-serve basis. Although this can reduce idle-mode energy consumption, this can cause more energy consumption on the part of ITs if the assigned berths are far from the dedicated yard areas for containers from vessels [12]. Therefore, recent literature has considered the tradeoff between energy-saving objectives and time-saving objectives when making plans for vessel handling operations [13].
There is a tight interrelation between the BAP, the QCAP, and the ITAP [14]. Existing studies have mostly focused on combining the BAP and the QCAP into a unified model called the integrated berth allocation and crane assignment problem (BACAP) [12,[15][16][17]. For a detailed review on BACAP studies, we refer the readers to the recent survey in Reference [5]. Energy-saving objectives are often implemented by adapting green practices and policies to reduce carbon dioxide emissions [3]. Few studies have addressed the BACAP while considering energy savings to be a part of objective functioning [11,13,18,19]. These studies have aimed to create effective handling plans that achieve a good service level (minimization of the total delay in departure for vessels) as well as energy savings (minimization of energy consumption) in the operations of vessels and QCs. Reducing the amount of energy consumed by activated QCs can be achieved by reducing the total departure delays of all vessels [16]. This is because minimizing the total departure delays of all vessels can lead to minimizing the unproductive times of QCs, resulting in improving QCs' energy consumption. Other studies have considered the energy consumption of QCs through the minimization of the number of QCs assigned to the vessels [13,19]. In order to minimize vessel fuel consumption while at port, relevant studies on BACAP have followed an emissions mitigation strategy that requires minimizing vessel waiting times [11].
In the majority of the works mentioned above, estimations of vessel handling times were made based on the number of containers to be handled, the assigned number of QCs, and the QCs' average handling rates (container moves/hour). It is worth noting that QCs' handling rates are highly dependent on the assigned number of ITs [20]. However, most studies have assumed that ITs are enough to supply the activated QCs with their truck needs [19,21]. Thus, the ITAP is ignored by this assumption.
Taking IT capacity constraints into account when planning handling operations has become crucial in adapting to recent changes in maritime logistics [22]. In the following section, we illustrate some examples of these recent changes. Many countries have recently witnessed a significant shortage in available truck drivers [23], which in turn constrains the capacity of ITs. As another example, recently it was noted that the port of Hong Kong is in serious need of more truck drivers because 59% of employed drivers have approached their retirement age and replacements are lacking [24]. Additional examples include deactivating some ITs due to maintenance activities [25]. On the other hand, labor costs are often in the range from 50% to 75% of total operation expenses in seaports [26]. Thus, it is not a cost-effective investment to maintain a complete set of drivers and/or ITs, especially when vessels' calling frequency fluctuate greatly [27].
In the literature, the shortage of ITs is rarely considered. In general, outsourcing and sharing IT strategies have been most frequently proposed as solutions for truck shortages. Wang et al. [25] have developed a mathematical formulation for truck scheduling and storage space allocation problems with consideration of outsourcing ITs to external companies. Their proposed model seeks to find the daily number of outsourced trucks to minimize the shift between the planned and actual completion times of vessels. Sharing ITs suits ports with multiple container terminals whose locations must be adjacent to each other, making it possible to move ITs between terminals in little time. He et al. [4] created a mathematical formulation to make decisions on IT sharing between multiple container terminals in a large port. Their formulation can be used to find the number of ITs that can be shared between terminals in a given time period. More recently, Karam and Attia [28] presented a mathematical model for truck scheduling, taking into account IT sharing and outsourcing strategies to solve the issue of a lack of ITs in adjacent container terminals. However, if the operation planner cannot solve the issue of a lack of ITs using sharing and outsourcing strategies, planning handling operations has to be done taking into account the shortage of ITs. In this case, because the traditional method of estimating the vessels' handling times ignores truck capacity constraints, the estimated handling times will be unreliable and often are lower than the true values. This in turn will delay the handling operations of successive vessels and may also delay the departure dates of the vessels. Therefore, the operation costs may increase due to delay penalties and higher energy consumption.
Very few studies have considered the shortage of ITs when making plans for berth allocation and quay crane assignments in container terminals. Karam et al. [29] developed a mathematical model for integrating the QCAP as well as the ITAP simultaneously. They showed the benefits of their proposed model compared to the traditional planning method. In an extended version of the model, they considered how to determine the specific QCs assigned to each vessel and proposed a heuristic-based solution to solve real problems [30]. More recently, Vahdani et al. [22] presented a two-objective optimization approach for simultaneous integration of the QCAP as well as the ITAP that considered the sharing of ITs between adjacent container terminals. Their model aimed to overcome the limited availability of ITs by sharing them among different terminals. None of these previous studies considered the BAP when integrating the QCAP and the ITAP. To date, integrating the BAP, the QCAP, and the ITAP has only been studied by Karam and Eltawil [31]. They developed a mathematical model for an integrated solution of BAP, QCAP, and ITAP. Their model aims to minimize the average service time of arriving vessels. In addition, they validated their model by comparing it to the BACAP model proposed by Meisel and Bierwith [12], assuming that ITs are available. In addition, their model was validated through a comparison to a queuing model proposed by Kang et al. [32] under limited availability of ITs. Based on the current state of the literature, it is obvious that the lack of integration of the BAP, QCAP, and ITAP, while considering energy consumption, is a gap in the existing literature. This paper presents a new approach that simultaneously integrates the BAP, the QCAP, and the ITAP with important practical aspects, such as the berth deviation factor, the handling time estimation, energy-saving goals (for vessels, QCs, and ITs), and limited time-variant quay crane assignments. Furthermore, a Lagrangian relaxation-based method was developed for handling problems of real size. The integration of these dependent problems led to an integrated berth allocation, QC assignment, and IT assignment problem that we named as the berth allocation crane assignment truck assignment problem (BACATAP). The proposed integrated approach considers an operational bottleneck problem (the limited availability of ITs) and different objectives such as the service level for vessels (minimizing delayed departures) and emission-related externalities from vessels at port, QCs, and ITs. To date, environmental aspects have not been considered in the integrated planning of the BAP, QCAP, and ITAP. In addition, this paper provides terminal planners with a new approach that can support the achievement of energy and cost savings at the operational level with a limited availability of ITs.
The rest of this paper is structured as follows: Section 2 presents the modeling approach. Section 3 illustrates the proposed Lagrangian method. The numerical experiments are explained in Section 4. The last section presents conclusions.

The Modeling Approach
This section presents a problem description with an illustrative example and the proposed mathematical model for the BACATAP.

Problem Description
For a known set of vessels V coming during a given planning horizon T and a known set of available QCs Q and ITs I, the solution of the BACATAP describes the berthing positions, berthing times, and completion times of all vessels V. Moreover, it identifies the QC assignment profiles for all incoming vessels and the IT assignment profiles for the QCs serving each vessel. From a practical point of view, the main objective of the BACATAP is to minimize the departure delays of arriving vessels. Due to increasing concerns for green ports, the cost to the environment was considered through additional objectives focused on minimizing the energy consumption of the handling equipment (QCs, YCs, and ITs) and all vessels (while at the port) in the planning horizon.
To illustrate our approach, an instance is solved using the traditional BACAP with the condition that there are not enough ITs to supply each QC with five trucks. Figure 1 shows a berthing and QC assignment plan for five vessels with 6 QCs and 30 ITs. As can be seen, the solution can be illustrated in a space-time diagram in which the quay as well as the planning horizon are partitioned into equal segments (berthing positions and time periods). In Figure 1, a vessel is represented by a rectangle whose left lower corner represents the berthing time and berthing position of the vessel, while the length as well as the handling time of the vessel are indicated by the height and width of the rectangle, respectively. In a feasible solution, the rectangles must not overlap with each other. Moreover, the assignment of QCs and ITs to each vessel can be found in the first and second rows, starting from the bottom of the vessel rectangle.
As shown in Figure 1, periods 4 and 5 are peak periods. If there is a shortage in the IT requirement of five ITs, the handling rate of the QCs in the peak periods will be lower than what is projected in the BACAP plan. Thus, this will increase the handling times of vessels A and B, which are served at peak periods. Moreover, this will also delay the service of subsequent vessels C, D, and E. Overall, this will impair service levels and cause more energy consumption due to increasing waiting times. To demonstrate the benefits of the developed approach, the example is solved by our proposed model for the BACATAP with 6 QCs and 25 ITs. The solution is shown in Figure 2, in which it can be noted that no delay in the service of the vessels occurs. The idea is that the proposed new model reallocates the QCs and ITs to vessel A such that the available capacities of the QCs and ITs are taken into account. Compared to the traditional BACAP, the proposed approach can produce shorter vessel service times and less energy consumption when there is a shortage of ITs.   The planning horizon is partitioned into equal time periods, and each period is one hour; 2.
The berth is discretized by being divided into equally sized segments of 50 m; 3.
The service of each vessel starts directly after berthing and cannot be stopped until the vessel's workload is fully handled; 4.
It is assumed that water depth along the berth is suitable for arriving vessels; 5.
Minimum and maximum numbers of QCs are determined for each vessel, and a vessel cannot be moored until the minimum number of QCs are ready to handle its containers; 6.
The change in the number of QCs serving each vessel is restricted to one change at most, as was done in Reference [33]. It is worth noting that strategies for the QC assignments include time-invariant QC assignments, time-variant QC assignments, and limited time-variant QC assignments. Time-invariant QC assignments imply that the number of QCs serving each vessel must be the same during the service of the vessel [16,34], while time-variant QC assignments imply that the number of QCs can be varied dynamically according to the assignable range of QCs allowed for the vessel [12,15,17]. Note that varying the number of QCs serving the vessel is sometimes essential for improving the utilization of QCs as well as the service times of vessels. This is because some QCs serving vessels are reallocated to start operations on other vessels. However, many changes cause a long setup time, which is unacceptable in reality. Therefore, a well-balanced tradeoff is limited time-variant QC assignments, in which the variation in the number of QCs is limited. Referring to vessel B in Figure 2, it can be noted that the assigned numbers of QCs change only one time during service in time period 7. Therefore, the service duration of the vessel is divided into two stages. A stage is the number of successive time periods, and the number of QCs assigned to the vessel does not change during the same stage. A new stage starts whenever a change in the number of QCs serving the vessel occurs. Referring to Figure 2, it can be noted that vessel B has two stages. The first stage runs from time period 4 to time period 6, while the second stage starts in time period 7, in which the assigned number of QCs changes from 3 to 2; 7.
Minimum and maximum numbers of ITs are specified for each QC so that the waiting times of the QCs and ITs are suitable. For more details on how such numbers are determined, we refer the reader to the work of Murty et al. [20], in which the authors developed a queuing network model to capture the interactions between QCs, ITs, and yard cranes. Their model can be used to determine the numbers of ITs assigned to each QC such that that the waiting times of the QCs and ITs are minimized.

Notations
zb jvt : 1, if j and t are equal to the first berthing position (r v ) and berthing time (S v ) of vessel v, respectively, and 0 otherwise.

Estimation of the Vessel Handling Time
In this section, we illustrate how to estimate the handling time of a vessel considering its assigned number of QCs and ITs. The handling time of the vessel mainly depends on the following factors: The assigned number of QCs; II. The assigned number of ITs serving each QC; and III. The distance between the vessel's berthing position and the storage yard for the containers on the vessel.
It is worth noting that a preferred berthing position b v is determined for a vessel v that is near its storage yard. If the actual berthing position of a vessel is far from its preferred position, the distances traveled by the ITs increase, and therefore, the terminal planner should employ more ITs to reduce the time the QCs wait for the ITs. Therefore, as can be seen in the work of Meisel and Bierwirth [12], if a vessel berthed ∆r v far from its preferred berthing position, its workload is increased to (1 + β ∆r v )N v TEUs. In practice, the effect of the far berthing position of the vessel can be alleviated by increasing its assigned numbers of QCs and/or ITs. Therefore, in our model, if a vessel is allocated to a berthing position far away from its preferred berthing position, its assigned numbers of QCs and/or ITs are increased. Through a simultaneous consideration of the three factors described above, the handling time for vessel v can be estimated as follows: where q i v is the number of QCs assigned to the vessel v if i ITs are assigned to each QC. To illustrate how this formula can be used, we provide the following example, in which we assume that a vessel can be assigned three QCs simultaneously and that each QC can be assigned four ITs. We assume further that N v , ∆r v , µ, τ, and γ are 300 TEUs, 10 berth segments, 0.035, 0.055, and 0.022 h, respectively. Here, µ + 2τ + γ designates the average time needed by an IT to transport one container from the QC serving a vessel to its storage block and to return to the QC to transport another container. Therefore, the average handling rate of one QC per one IT can be calculated by 1/( µ + 2τ + γ), expressed in containers per hour. If the berth deviation factor is 0.02, the handling time of the vessel is calculated to be approximately 5 hour, according to the above formula. Meanwhile, if i, q i v , and ∆r v are changed to 2, 3, and 0, respectively, the handling time of the vessel becomes 9.2 hour. Compared to existing approaches, this proposed handling time estimation enables the consideration of the three factors simultaneously, which greatly affects the vessel handling time. This can produce a more reliable estimate of the vessel handling time, which is necessary for reliable planning. Thus, our model can improve handling plans in terms of their reliability as well as their applicability.

Formulation of Objective Functions
As stated before, the proposed approach considers improving the service level as well as the energy consumption of vessels, QCs, and ITs. The service level is considered through the minimization of the total delay in departure of the vessels. According to He [13], the minimization of the delayed departure of all vessels can also ensure that the energy consumption of activated QCs is minimized. This is because the shorter the delay in departure of all vessels is, the less the unproductive time of the QCs is, which can consequently improve the QCs' energy consumption. The energy consumption of ITs mainly depends on achieving the shortest travel distance between quayside and storage yard [18]. These travel distances are highly affected by the distance from the vessel's berthing position to the storage yards specified for its containers [10,12]. Typically, the berthing position of a vessel is optimal if it is the one closest to the storage yard areas for the containers. If the actual berthing position is far from the optimal position, the travel distances of the ITs increase. This in turn results in high congestion and more fuel consumption. Figure 3 illustrates how traffic congestion occurs when containers on vessel A are to be stored in yard 2 while containers on vessel B are to be stored in yard 1. Such traffic congestion increases the amount of fuel consumed by trucks as well as lengthening the total handling time of the vessels. Therefore, we propose a second objective function that penalizes the deviation between the desired berthing positions and the actual berthing positions of the vessels. This minimizes the energy consumption of the ITs. Finally, we use a third objective function that minimizes the waiting times of all vessels, thus minimizing the fuel consumption of vessels at port following the work of Golias et al. [11].

The Mathematical Model
∆C The objective function of the proposed model aims to minimize the weighted sum of the penalty and energy costs of the handling operations. The objective function includes three cost parts. The first part penalizes deviating the vessel's berthing position to somewhere distant from the storage yard for its containers. The yard implies the desired berthing position. The second part reflects the energy consumption costs if there is a time delay between the berthing and arrival times of the vessels. The third part is a penalty cost that arises if the vessel's completion time occurs after its expected departure time. Constraint set (3) requires that in each time period, each berthing position is taken by only one vessel. Constraint sets (4) and (5) require that in each time period, the activated number of QCs and ITs do not exceed the available capacities of QCs and ITs. Constraint set (6) requires that the containers on each vessel must be entirely serviced. Constraint set (7) defines the variable x vt as the berthing duration and is calculated as the length of time between the berthing and completion times. This constraint also guarantees that the handling operation cannot be stopped once it begins. Constraint sets (8) and (9) describe the completion and berthing times for all vessels, respectively. Constraint set (10) ensures that within the planning horizon, each vessel must be assigned only one berthing position and one berthing time. Constraint sets (11)- (14) require that for each vessel, xb jvt = 1 only in the grids of the time-space diagram covered by the rectangle of vessel v. The grids, which are covered by the vessel, are [S v , . . . ,C v ] and [r v , . . . ,r v +l v -1]. Constraint sets (15)- (16) require that the vessel's assigned number of QCs must be within the specified boundaries. Constraint sets (17)-(18) define the variable zb jvt . Constraint sets (19)- (20) guarantee that the number of QCs serving a vessel is adjusted at most one time during vessel service. It should be noted that the change in the assigned number of QCs can be simply adjusted using constraint set (21). For example, if the port planner prefers to set the number of crane changes to one, the right hand side of constraint set (21) is set to three. This is because there is also a change in the number of QCs at the start and end of the handling operations for the vessel. Constraint sets (22)-(25) define the variables ∆C v , ∆S v , and ∆C v respectively. Constraint set (26) guarantees that the completion time of the vessel is larger than the planning horizon boundaries. Constraint set (27) guarantees that the berthing time of each vessel is not before its arrival time. Constraint set (28) requires that the berthing position of each vessel does not exceed the berth boundaries. Constraint sets (29)-(30) are the non-negativity and binary constraints, respectively. Obviously, the proposed model is not in the standard form for mixed integer programing due to the absolute value in constraint sets (19) and (20). Therefore, the linearization instructions presented by Bertsimas and Tsitsiklis [35] were followed to remove the absolute value. By introducing the variable lq vt , the mathematical model could be reduced to a mixed integer programming problem, as follows: subject to constraints (3)-(18), (21)- (30), and the following constraints: lq vt ≥ 0 and integer ∀v ∈ V; ∀t ∈ T.

The Lagrangian Solution Method
As a result of combining three already complex problems, the complexity of the developed model grew dramatically. An exact solver for mixed integer programming problems is able to handle up to three vessels within one hour, while an instance of five vessels cannot be solved by the same solver. Therefore, a Lagrangian relaxation-based method was developed. The design of the Lagrangian method was motivated by the following: constraints (3), (4), and (5) make the proposed model difficult to be solved due to the limited resources for berthing positions, QCs, and ITs. Relaxing these constraints leads to a much easier problem. Therefore, a relaxed version of the problem is obtained by relaxing these constraints and adding them to the objective function with multipliers. Note that each multiplier is a penalty that is incurred in the solution when its corresponding constraint is violated. It should be noted that the solution for the relaxed problem includes a set of vessels, and each vessel has an individual cost contributing to the total solution cost. Therefore, it is possible to decompose the relaxed problem into subproblems, and each subproblem is related to only one vessel. Then, the overall solution of the problem is found by combining the solutions to the subproblems together. The optimal solution obtained by solving the relaxed problem may be infeasible due to the relaxation of the three constraints, but the optimal value can be used as a lower bound on the optimal solution of the original problem. Therefore, a feasible solution to the original problem can be generated from the solution to the relaxed problem. An optimal lower bound can be calculated by using the subgradient optimization procedures illustrated in Reference [36], which also have been considered in recent literature on the BAP, the BACAP, and the integrated QCAP and ITAP [12,30,33,37]. In these works, resource capacity constraints were dealt with as hard constraints, which if relaxed, could lead to a much easier problem that can be solved efficiently in comparison to the original problem.

The Lagrangian Relaxation for the Developed Model
In this work, Lagrangian relaxation is achieved by relaxing constraints (3), (4), and (5). The relaxed problem (P) can be written as follows: subject to constraints (6)- (30) and π jt , π t , π t ≥ 0. Here, π jt , π t , andπ t are the Lagrangian multipliers of the relaxed constraints. By eliminating the constant parameters from (36), the relaxed problem can be written as follows: subject to constraints (6)-(30) and π jt , π t , π t ≥ 0. The Lagrangian dual (LD) of P is as follows: LD : Max π jt ,π t ,π t P P : Min subject to constraints (6)-(30) and π jt , π t , π t ≥ 0. To obtain the Lagrangian dual (LD v ) for the subproblem P v , P is decomposed into a set of subproblems, each of which represents a vessel. LD v is as follows: LD v : Max π jt ,π t ,π t P v P v : Min subject to constraints (6)-(30) and π jt , π t , π t ≥ 0.

The Overall Procedures for the Subgradient Optimization Method
The overall procedures described in the work of Geoffrion [36] are as follows: Step 1: Initialize the upper bound (Z UB ) as a large value and the multipliers π jt , π t , and π t as 0.
Step 2: The lower and upper bounds of the optimal solution objective are calculated as follows: Step 2.2: For each vessel v, solve LD v using initial values of the multipliers. The procedures for solving LD v will be presented in Section 3.3. The lower bound of the solution objective (Z LB ) is formed by (34), and the largest lower bound is updated, Z max LB = max Z max LB , Z LB . Step 2.3: If the solution obtained using (36) violates the relaxed constraints, a feasibility restoration method is used to solve the infeasibilities and calculate the upper bound (Z UB ). Afterwards, update the lowest upper bound, Z min UB = max Z min UB , Z UB .
Step 3: Evaluate the subgradients as well as the step sizes to calculate the multipliers π jt , π t , and π t .
The subgradients G jt , G t , and G t can be calculated as follows: The step sizes can be calculated as follows: (43) where λ is a variable for adjusting the step sizes; and λ has a value in the range from 0 to 2. In the first iteration of the procedures, λ is initialized to 2 and is halved as long as there is no improvement in the lowest upper bound (Z min UB ) for a specified number of successive iterations. The procedure is aborted if λ reaches a value below 0.005 or the duality gap between Z min UB and Z max LB is below one.

The Overall Procedures for the Subgradient Optimization Method
The following section presents the proposed systematic procedures for obtaining solutions to the decomposed primal problem (LD v ) for vessel v. The solution for LD v describes the berthing position, berthing time, QC assignments, IT assignments, and completion time for vessel v. Before introducing the detailed steps of the systematic procedures, an example is provided to illustrate some terms to be used in the explanation. Consider vessel A, which has N A = 400, q max A = 4, l A = 5, a A = 3, and b A = 4. Furthermore, the terminal data are as follows: i max = 5, H = 8, and L = 9. As shown before, the vessel is represented in the space-time diagram by a rectangle. Now, it is required to represent the feasible region within which the left lower vertex (or reference corner) of the vessel's rectangle can be located. Figure 4 shows the feasible region of vessel A, which is shaded in gray. The feasible region is determined through its respective spatial and temporal constraints. The symbol w A in Figure 4 represents the longest time period during which vessel A can be berthed and is calculated as the rounded (up) value of H − N A pr * i max * q max A + 1, where pr = 1 µ+2τ+γ . Moreover, the time-space diagram can be divided into grids, and each grid can be described by ( j, t). It can also be noted that (b A , a A ) represents the perfect location of the left lower vertex of the vessel rectangle, as the berthing cost (the sum of the second and third terms of the objective function) is zero. However, when the container terminal is congested, the reference corners of some vessels' rectangles deviate from their ideal locations due to the limited availability of berthing positions, QCs, and ITs. The proposed procedures first determine the lowest cost in terms of berthing time and berthing position. Then, assignments of QCs and ITs are generated for the identified berthing position and time. After that, the objective function in (39) is calculated for the current schedule (berthing time, berthing position, completion time, and assignments of QCs and ITs).
The procedures, from step 1 to step 3, are performed for all feasible combinations of ( j, t) at which the reference corner of vessel v can be located, where j = 1, . . . , L − l v + 1 and t = a v , . . . , H − N v pr * i max * q max v + 1.
Step 1: Assign an initial value to the minimum objective value, e.g., 15,000, and set j = 1 and t = 1. Step 2.1: Select the next grid, starting from b v through L − l v + 1. If the berthing cost in the grid is lower than the minimum objective value, go to step 3. Otherwise, go back to the beginning of step 2.1 until all berthing positions from b v through L − l v + 1 are tried once, and then go to step 2.2. It is important to point out that the minimum objective value is iteratively updated in step 3 based on the vessel schedule (berthing time, berthing position, completion time, and assignments of QCs and ITs). Therefore, if the berthing cost is higher than the minimum objective value, this implies that berthing the vessel in the grid corresponding to this berthing cost will not yield an objective value better than the current minimum objective value. Therefore, this grid is excluded, and the next grid is tried. In this way, step 3 is not performed for all the grids in the feasible region, and thus computational time can be saved.
Step 2.2: Select the next grid, starting from (b v − 1) to 1. If the berthing cost in the grid is lower than the minimum objective value, go to step 3. Otherwise, go back to the beginning of step 2.2 until all berthing positions from b v − 1 to 1 are tried once, and then go to step 2.
Step 3: The reference corner of the rectangle corresponding to vessel v is located at (j,t), which is determined by step 2. Hereafter, in this step, the QCs and ITs are allocated to vessel v. To stimulate limited time-variant quay crane assignments, the service of the vessel is divided into two stages. The first and second stages are designated by t 1 and t 2 , respectively. To try all feasible combinations of the QCs and ITs that can be allocated to vessel v, the following steps are performed.
Step 3.1: If all combinations of (q 1 , i 1 ), q 1 = q min v , . . . , q max v , and i 1 = i min , . . . , i max are tried once, go back to step 2. Otherwise, select the next combination of (q 1 , i 1 ) and set t 1 = 1, where t 1 is the first stage at which q 1 QCs are assigned to vessel v and i 1 ITs are assigned to each QC.
Step 3.2: For all combinations of (q 2 , i 2 ), q 2 = q min v , . . . , q max v , and i 2 = i min , . . . , i max , select the next combination of (q 2 , i 2 ) and set t 2 to the rounded value of max 0, , where t 2 is the second stage at which vessel v is assigned q 2 QCs and i 2 ITs are assigned to each QC so as to completely handle N v . If t + t 1 + t 2 -1 ≤ H, then obtain the objective value of (39) corresponding to the schedule just identified for the combination of {j, t, t 1 , q 1 , i 1 , t 2 , q 2 , i 2 }, and update the minimum objective value. If two or more schedules have the minimum objective value, random selection is used to select the vessel's schedule from these schedules. Otherwise, go to step 3.2 until all combinations (q 2 , i 2 ) have been tried once.
Step 3.3: If t 1 ≤ N v pr * i 1 * q 1 , increase t 1 by one and go to step 3.2. Otherwise, go to step 3.1.
The lowest objective value determined from step 3.2 represents the best solution for LD v . Equation (36) is used to obtain the lower bound of the optimal solution objective. It should be noted that the obtained solution may violate the relaxed constraints; therefore, the solution may be infeasible. In this case, another method is utilized to make the solution feasible.

Feasibility Restoration and Solution Refinement
As stated before, it may happen that some iterations of the subgradient method find an infeasible solution. In this case, the following procedures are used to change the infeasible solution into a feasible solution.
Step 1: Sort vessels sequentially in increasing orders of a v and set v = 0 and Fe = In = ∅, where Fe and In are feasible and infeasible sets, respectively.
Step 2: Increase v by 1 and perform step 2 if v ≤ n. Otherwise, go to step 3. In this step, the obtained solution of the vth vessel is checked for any violation of the relaxed constraints, given that the solutions for the preceding vessels in Fe are known. If there is no violation of the relaxed constraints, the vth vessel is inserted into Fe. Otherwise, the vth vessel is inserted into In.
Return the beginning of this step.
Step 3: Sort vessels in In sequentially in increasing orders of a v , and then the next vessel in In is selected. The procedures in Section 3.3 are used to find the least-cost feasible schedule for the selected vessel that results in no violation of the relaxed constraints, provided that the solutions of vessels in Fe are given. Afterwards, the selected vessel is deleted from In and is inserted into Fe. Redo step 3 until In contains no vessels.
Step 4: After a feasible solution is constructed by performing the previous steps, the feasible solution is refined. As noted, the solution is composed of a set of vessels. Each vessel has its own cost that contributes to the overall quality of the obtained solution. Therefore, if the schedules of vessels with relatively higher costs are refined, the overall solution cost may be improved. A straightforward procedure is applied as follows: Sort vessels sequentially in decreasing order of their individual cost. For the first 20% of vessels at the top of the list, perform the same procedure of solving the decomposed primal problems. It is important to point out that a percentage of 20% is selected because if this percentage were higher, more vessels will have to carry out the procedure of solving, which would lead to an increased computational time.

An Illustrative Example
To illustrate the proposed heuristic algorithm, the results of some iterations during solving a numerical example of the three vessels A, B, and C are provided in Figure 5. In this numerical example, the number of available QCs and ITs is set to 8 and 35, respectively. Moreover, the minimum and maximum limits of the ITs are 3 and 5, respectively. This implies that there is a shortage of five ITs, such that it is impossible to provide each QC with its maximum number of ITs when all of the available QCs are activated. Figure 5a illustrates the first iteration of the developed subgradient optimization method when all multipliers have values of zero. Note that in the first iteration, each vessel is berthed at its ideal location, (b A , a A ). As can be seen in Figure 5a, there is a violation of the crane capacity constraint (4) as well as the truck capacity constraint (5). The violation is indicated in red in Figure 5a. In addition, the rectangles overlap with each other. This means that constraint (3) is also violated. During the iterative procedure of the subgradient optimization, the multipliers corresponding to the grids, in which constraints (3), (4), and (5) are violated, increase. The objective function value of LD v in (39) is a function of the multipliers. Therefore, the objective function value belonging to the vessels, which are berthed in these grids, increases, leading to a reallocation of berthing positions, berthing times, QCs, and ITs to these vessels. The results, with updated multipliers, are shown in Figure 5b,c, in which it can be noted that violations of constraints (3), (4), and (5) are reduced compared to the first iteration in Figure 5a. Figure 5d shows the final solution for the numerical example.

Numerical Experiments
This section includes four parts. The first part shows how test instances are generated. In the second part, a number of small test instances are solved using the developed Lagrangian method, and the results are validated against the results produced by the ILOG CPLEX TM solver. Numerical experiments with large problem sizes are carried out in the third part. In the fourth part, we illustrate the added value of the proposed approach using a comparison to an existing approach. All of the experiments were carried out using a computer with a 2.5-GHz processor and 1 GB of RAM. The Lagrangian method was coded on MATLABTM.

Generating Test Instances
To generate realistic instances for a performance assessment of the developed Lagrangian method, three vessel classes, i.e., feeder, medium, and jumbo, were created, as in Reference [12]. The vessel classes have different technical characteristics, as illustrated in Table 1. In each generated instance, feeders represent 60% of the vessels, medium vessels represent 30% of the vessels, and jumbo vessels represent 10% of the vessels. The planning horizon's length is set to one week and is discretized into one-hour periods. The quay length is assumed to be 1000 m and is partitioned into 50-meter segments. Thus, the available number of berthing times (H) and berthing positions (L) is 168 and 20, respectively. The expected arrival times, the preferred berthing positions, and the expected departure times are generated from U(1, 168), U(1, L − l v + 1), and the rounded value of a v + N v /(30 * q max v ) + U(5, 10), respectively. The remaining model parameters are as follows: Q = 10, I = 40, i max = 5, i min = 3, and β = 0.01. Moreover, µ, τ, and γ are 0.035, 0.055, and 0.022 hour, respectively. Finally, ca 1 v , ca 2 v , and ca 3 v were estimated to be $1000/unit distance deviation, $1000/hour, and $1000/hour, respectively. However, the developed model can be solved using different weight values. It should be noted that different values of ca 1 v , ca 2 v , and ca 3 v will clearly influence the results. For example, decreasing ca 1 v will result in allowing vessels to be moored at any available berthing position along the quay. This in turn increases the deviation between the expected and actual berthing positions, and thus, more energy is consumed by ITs traveling between the quay side and storage areas. It is worth noting that larger value of ca 1 v is preferred for large container terminals. This is because the travel distance between the quay side and the storage areas is significant in large container terminals compared to small ones. Setting the value of ca 2 v to be higher than ca 1 v and ca 3 v will force the model to decrease emissions from vessels by minimizing their waiting times at the terminal. This means that large value of ca 2 v allows for serving vessels almost on a first-come, first-served basis. Large values of ca 2 v are preferred for container terminals serving vessels whose idle emissions are relatively high, such as vessels whose workload is mostly refrigerated containers. Setting the value of ca 3 v higher than ca 1 v and ca 2 v will ensure that most vessels depart from the terminal according to their promised departure times. This suits container terminals whose handling equipment is energy-efficient, and thus, such terminals pay more attention to improving service levels by using large values of ca 3 v . Future research should focus on more detailed analyses of ca 1 v , ca 2 v , and ca 3 v . Here, equal values are used to illustrate the proposed methodology.

Comparison to Cplex Results
To judge the quality of the solutions found using the developed Lagrangian method, the solutions determined by the developed method were compared to the optimal solutions obtained from solving the proposed model using ILOG CPLEX TM (IBM, New York, NY, USA). For instances of small sizes, the optimal solutions were found by CPLEX in short computational times. However, the computational time became relatively long with increasing instance sizes, and in some cases, it was not possible to find the optimal solution due the limited physical memory of the PC. Consequently, this section presents the solutions of small size instances. Note that each instance was solved 10 times using the Lagrangian method, and the average upper and lower bounds for each instance were calculated. The computational results for the comparison are given in Table 2. In these experiments, the planning horizon was 24 hour (one day). It is worth noting that the average upper bound represents the feasible objective value (final solution) that could be found using the Lagrangian method. As shown in Table 2, the Lagrangian method could find the optimal solutions for the instances that could be solved to optimality by solving the proposed mixed integer programming model using ILOG CPLEX TM . To judge the quality of the lower bounds found using the Lagrangian method, the lower bounds found using the proposed procedures to solve the decomposed primal problems (Section 3.3) were compared to the optimal solutions found by ILOG CPLEX TM . The results illustrate that the lower bounds found using the Lagrangian method are lower than the optimal solutions found by ILOG CPLEX TM in all instances. This implies that the procedures for solving the decomposed primal problems can produce optimal lower bounds, and thus the relative gap "Gap %" between the average upper bound and average lower bound can be used as a measure of the performance of the developed algorithm.

Solving Instances with Large Sizes
To check the performance of the developed Lagrangian method in solving generated instances of a realistic size, three sets of instances, including 20, 25, and 30 vessels, were randomly generated. Each set contained 10 instances. The numerical experiments for the differently sized instances can be seen in Table 3. The worst, best, and average values of the upper bounds, lower bounds, and gaps, as well as the average computational time (obtained after running the method 10 times), are shown in Table 3. As can be noted from Table 3, the gap and the computational time increased as the problem size increased. On average, the gap and computational time were 15.66% and 16.66 min, respectively. However, we can observe that the average gap and average computational time were more than 30% and more than 40 min, respectively, for some instances of 30 vessels. This is because the higher the number of vessels is, the more the degree of overlap between the vessels becomes significant. This makes such instances difficult to be solved with limited numbers of berthing positions, QCs, and ITs. Thus, during the iterative procedure for the subgradient optimization, the multipliers related to the grids, where the vessels overlapped with each other, were increased. This led to moving the overlapping vessels to unallocated grids away from their ideal locations, resulting in an increase in the objective value for some vessels: to try more grids in the feasible region and in each grid, steps 3.1 to 3.3 are performed, which increases the computational time.
In Table 3, it can also be observed that the average statistics (in the last rows of the table) for the upper and lower bounds for the worst, best, and average cases were slightly different. This implies that in terms of consistency, the efficiency of the proposed method is fairly high. Figure 6a shows the lower bound and the corresponding upper bound generated during the iterative procedure for instance 14. The convergence of the heuristic for instance 14 is also presented in Figure 6b. The upper and lower bounds converged after 56 iterations, resulting in a solution of 24,000 with an optimality gap of 16.67%.
The performance of the developed Lagrangian method was considered to be satisfactory in solving the realistically sized instances, considering the increasing difficulty of the research problem as the problem size increase. Moreover, it should be noted that we solved these instances with a shortage in trucks, which also increases the problem difficulty.

Significance of the Proposed New Integration Framework
In this section, we demonstrate the significance of the proposed integration framework for the BAP, QCAP, and ITAP by comparing our approach to an existing approach. We solved an instance by using our approach and the approach of Ursavas [17] with different numbers of ITs. The main reason for specifically choosing the approach of Ursavas [17] is that she considered QCs with different handling rates, which facilitated a comparison to our approach, as will be illustrated later. In the following section, we refer to the approach of Ursavas [17] as the existing approach. To compare the two approaches to each other, we modified the objective function of the existing approach to be the same as our objective function. In order to solve the instance with different numbers of ITs using the existing approach, the handling rates of the QCs were reduced according to the truck shortage, as follows: 30 TEUs/h for a QC with five ITs, 24 TEUs/h for a QC with four ITs, and 18 TEUs/h for a QC with three ITs. For example, if the truck capacity was reduced from 50 to 45, then five QCs would have a handling rate of 24 TEUs/h because these five QCs would be assigned four ITs instead of five ITs, while the remaining QCs would be assigned five ITs and would work with a handling rate of 30 TEUs/h. Figure 7 shows the total handling costs obtained using the two approaches for this instance if the number of ITs varied from 40 to 50 and all other parameters were held constant. The two approaches could obtain a total handling cost of $30,000 with no shortage in trucks (I = 50). It can be observed that compared to our approach, the handling costs found by the existing approach obviously increase with decreasing numbers of available ITs. Interestingly, in some cases, e.g., I = 45, 44, and 43, the handling costs obtained using our approach are the same, while the handling costs obtained using the existing approach increase with decreasing available ITs. The main reason is that our proposed approach can handle the shortage in ITs such that this shortage does not affect the handling cost of the vessels, as illustrated in detail in Section 2.1. On average, our approach could improve the solutions for the instances with different numbers of ITs by 6% compared to the solutions obtained using the existing approach. This ensures that the integration of the BAP, QCAP, and ITAP has a significant added value for container terminals. In addition, integrating the BAP, QCAP, and ITAP leads to a more reliable estimation of vessel handling times, and therefore, such integration is an essential factor in managing handling operations in container terminals.

Conclusions
The operational efficiency of container terminals (e.g., short service times and low operating costs) is highly dependent on effective planning of the handling operations of vessels using the available berths, QCs, and ITs, which are scarce and expensive resources. Because significant overhead costs for container terminals are also caused by energy consumption, it is essential to reduce energy consumption to achieve cost reduction and higher profits. Therefore, recent studies have considered the tradeoff between energy-saving objectives and time-saving objectives in planning handling operations for container terminals.
This paper presented a mixed integer programming model for integrating the BAP, the QCAP, and the ITAP, considering energy-saving goals. As a result of integrating three already complex problems, the complexity of the resulting model increased dramatically, and therefore, there was a need for an efficient solution method. Thus, a Lagrangian relaxation-based method was developed to solve the model.
The main contribution of the presented model is the integration of the BAP, the QCAP, and the ITAP, considering the energy-saving goals of vessels, QCs, and ITs. In addition, we estimate the vessel handling time considering the number of QCs assigned to the vessel, the number of ITs assigned to each QC, and the distance between the berthing position of the vessel and its storage yard. Therefore, compared to the traditional BACAP, the integrated model yields better estimates of vessel handling times and enables reduction in energy costs. Furthermore, rather than renting external trucks to compensate for the lack of ITs as suggested by existing studies, which would increase the operating and energy costs of the handling plan, the integrated model can handle the limited availability of ITs and provides a handling plan with a lower cost than what would be obtained by solving these problems separately.
Numerical experiments were carried out to assess the performance of the developed Lagrangian method. The solutions of the developed method were validated against those found by an exact solver. The exact solver failed to reach an optimal solution or a lower bound for instances of a practical size, while the developed method could obtain an upper bound and a lower bound for all test instances. This reflects the benefits of developing a heuristic solution method for instances with a large size to provide a good solution within a practical amount of time. The developed method could obtain near optimal solutions with an average gap of 15.66% and an average computational time of 16.60 min.
In addition, the results showed that compared to the existing approach, the solutions obtained by our approach in the case of a shortage of ITs caused no delay in the service of the vessels. The idea is that the proposed new model reallocates the QCs and ITs to vessels in a way that the available capacities of the QCs and ITs are not violated. This enables the production of better solutions in the case of a limited availability of ITs, and this therefore confirms the need to integrate the BAP, QCAP, and ITAP, as proposed in this work.
In future research, the developed model can be extended to determine the assignment of QCs to each specific vessel. In the current work, the solution time was a bit consuming. Therefore, other heuristic and metaheuristic approaches can be investigated to improve the solution time of the proposed problem. Finally, future work may address the possibility of alternative fuels for IT and QC engines/motors, or autonomous vehicles that would reduce labor costs and limitations in terms of drivers. All of these proposals could save on costs and emissions.

Conflicts of Interest:
The authors declare no conflicts of interest.