Variable Neighborhood Search for Multi-Port Berth Allocation with Vessel Speed Optimization

: This paper delves into the multi-port berth allocation problem (MBAP), enriching the traditional berth allocation problem (BAP) with vessel speed optimization (VSO). In the MBAP, it is assumed that there is cooperation between the port and the shipping companies, and the operation of the vessels and the ports is planned to maximize the overall benefits. Exploring this potential collaboration between ports and shipping entities has the potential to mitigate, or even resolve, the challenges plaguing maritime transportation, e.g., port congestion and suboptimal vessel schedules, to ultimately enhance the efficiency of maritime trade. In this paper, a new mixed-integer linear programming (MILP) model for the MBAP is formulated, which attempts to minimize the total cost incurred during operations, with various constraints such as vessel sailing, the vessel space–time relationship in ports, and the planning period. Meanwhile, an innovative variable neighborhood search (VNS) algorithm is presented, in which the initial solution generation method and neighborhood structures are proposed according to the MBAP characteristics. Furthermore, two sets of MBAP instances are generated to test the proposed MILP and VNS, of which the first set is based on real-world port data and the second on existing studies. The numerical experiments verify that the VNS can efficiently and reliably solve instances of all scales, with each neighborhood structure contributing uniquely to the iterative process. In addition, by analyzing the impact of varying oil prices on the MBAP, the study offers valuable management insights. Finally, a case study based on real data from a port group in the Yangtze River Basin is presented to further demonstrate the necessity of considering vessel service time window and planning period in the MBAP as well as the important role of the VSO in scheduling.


Introduction
According to the United Nations Conference on Trade and Development (UNCTAD), shipping and trade patterns have begun to change since 2023, influenced by global trade policy and the geopolitical situation [1].Meanwhile, the global supply chain crisis has not fully recovered, the maritime market is weak, and the global shipping industry continues to face multiple challenges.As maritime trade began to gradually recover, shipping companies and ports struggle to cope with the surge in throughput.In addition, the UNCTAD has proposed higher standards for the transition to the decarbonization of maritime transportation.Usually, fleet renewal and port expansion are the most straightforward ways to solve shipping dilemmas such as unpunctual vessel schedules and port congestion; however, the high cost is prohibitive.Therefore, the technical improvement in the scheduling level has become the main means to solve these problems.
Port scheduling planning needs to complete the service of arriving vessels, among which the most basic port equipment is berths.Since Lim [2] first defined this problem as the BAP, this problem has been a research hotspot in the past two decades.Early research mainly focuses on the model formulation and algorithm innovation of the discrete and continuous BAP [3][4][5].Further, more efficient BAP models were developed, such as the multi-objective programming model [6] and bi-layer programming model [7].According to the reviews on the BAP, Bierwirth and Meisel [8] stated the heuristic algorithm is dominant in solving algorithms [9][10][11].In addition, Li et al. [12] illustrated that more problem features are considered in the BAP, such as tidal effects [13] and service order [14,15].Considering environmental sustainability, some scholars have begun to pay attention to emissions during vessel and port operations.Among them, some studies introduce the goal of minimizing carbon emissions in port resource scheduling problems [16][17][18], and some scholars study the promising emissions reduction technology of on-shore power supply [19][20][21][22][23].In addition, some studies have introduced the VSO into BAP-related problems to reduce emissions during vessel sailing by optimizing vessels' speeds.Hu [24] studied the relationship between delay times and emissions during vessel sailing and mooring in the BAP, and the experimental results showed that improving the vessel handling efficiency helps to reduce the delay time and emissions.Xia et al. [25] proposed a ship scheduling method with speed reduction to reduce emissions and tested its effectiveness through experiments.Yu et al. [26] proposed a bi-layer multi-objective model by integrating the berth allocation and quay crane assignment problem and vessel speed optimization problem and considering vessel service differentiation.They further developed a nested genetic algorithm, and the experimental results showed that this model can effectively reduce vessel emissions.However, the above studies are still limited to the planning of a single port.
With the evolution of maritime transportation, Venturini et al. [27] realized that independent port scheduling could no longer meet the actual trade demand, and for the first time extended the BAP to multiple ports and introduced the VSO and then proposed the MBAP, which assumes established port-shipping cooperation.Although they proposed a mixed-integer programming (MIP) model for the MBAP, an effective solving method was lacking.Consequently, other scholars further studied the MBAP and proposed the branchand-cut-and-price algorithm to solve the MBAP [28], a branch-and-cut method [29], and an adaptive large neighborhood search algorithm [30] to solve the MBAP with continuous berth layout.
In summary, research on the MBAP remains limited.Building upon existing MBAP studies [27][28][29][30], this paper formulates a new MILP model that incorporates the vessel spacetime relationship constraints proposed by Türko gulları et al. [31].Additionally, a VNS algorithm is proposed, tailored to the characteristics of the MBAP through modifications to the initial solution generation method and neighborhood structures.Furthermore, due to the scarcity of publicly available instances, two sets of MBAP instances based on real port data [5] and Agra and Rodrigues [32] are generated.These two sets of instances are used in VNS performance comparison experiments, and some instances in the first set are further applied to oil price sensitivity analysis, yielding management insights based on the results.Finally, an MBAP case study based on real data from a port group is performed and illustrates the importance of enforcing time dimension constraints in scheduling.
In what follows, Section 2 describes the problem characteristics of the MBAP.Section 3 presents the formulated MILP model, while Section 4 elaborates the proposed VNS algorithm.Section 5 showcases all the numerical experiments, and Section 6 provides the conclusion.

Problem Statement
The MBAP is an amalgamation of the BAP and the VSO, encompassing berth allocation across multiple ports while considering the speed decisions of vessels between ports.Based on existing studies [27,28], the MBAP assumes that there is cooperation between ports and shipping companies, facilitating information sharing and striving for the maximization of common interests.
Figure 1 illustrates the main components of the MBAP, which comprises vessels and ports.Specifically, adhering to the designated port sequence, vessels sequentially visit ports, completing respective loading and unloading operations.Meanwhile, the port allocates berths for arriving vessels and determines the commencement of work.This process can be regarded as a combination of multiple BAPs, typically depicted using a two-dimensional space-time diagram, as shown in Figure 2. Here, the vertical axis is the time axis, and the horizontal axis is the space axis (the berths of the port).The rectangles symbolize vessels, with their lower end indicating the start time and the upper end representing the departure time.Furthermore, the overlapping rectangle of vessels 4 and 5 in berth 3 illustrates the space-time conflict between vessels.ports and shipping companies, facilitating information sharing and striving for the maximization of common interests.Figure 1 illustrates the main components of the MBAP, which comprises vessels and ports.Specifically, adhering to the designated port sequence, vessels sequentially visit ports, completing respective loading and unloading operations.Meanwhile, the port allocates berths for arriving vessels and determines the commencement of work.This process can be regarded as a combination of multiple BAPs, typically depicted using a two-dimensional space-time diagram, as shown in Figure 2. Here, the vertical axis is the time axis, and the horizontal axis is the space axis (the berths of the port).The rectangles symbolize vessels, with their lower end indicating the start time and the upper end representing the departure time.Furthermore, the overlapping rectangle of vessels 4 and 5 in berth 3 illustrates the space-time conflict between vessels.The relationship among multiple BAPs is established through vessel sailing, and the VSO is further introduced.In the MBAP, the vessel speed is directly related to the fuel consumption and therefore becomes an important part of the total cost.In this paper, the relationship between ship speed and fuel consumption is calculated according to the formula (1) proposed by Venturini et al. [27], where ℓ is the vessel fuel consumption per unit distance at the speed ,  is the designed speed, and ℱ is the vessel fuel consumption per hour at the design speed.
Based on the above description, the MBAP follows the following assumptions.
(a) The vessel berthing time window is composed of a strict start time window and a soft departure time window; that is, the vessel must start operations after the expected start time but can leave after the expected finish time, the latter of which will incur a delay penalty.(b) The potential speed of vessels is discretized over a given range [27].(c) Berth depth and width can accommodate all vessels.ports and shipping companies, facilitating information sharing and striving for the maximization of common interests.Figure 1 illustrates the main components of the MBAP, which comprises vessels and ports.Specifically, adhering to the designated port sequence, vessels sequentially visit ports, completing respective loading and unloading operations.Meanwhile, the port allocates berths for arriving vessels and determines the commencement of work.This process can be regarded as a combination of multiple BAPs, typically depicted using a two-dimensional space-time diagram, as shown in Figure 2. Here, the vertical axis is the time axis, and the horizontal axis is the space axis (the berths of the port).The rectangles symbolize vessels, with their lower end indicating the start time and the upper end representing the departure time.Furthermore, the overlapping rectangle of vessels 4 and 5 in berth 3 illustrates the space-time conflict between vessels.The relationship among multiple BAPs is established through vessel sailing, and the VSO is further introduced.In the MBAP, the vessel speed is directly related to the fuel consumption and therefore becomes an important part of the total cost.In this paper, the relationship between ship speed and fuel consumption is calculated according to the formula (1) proposed by Venturini et al. [27], where ℓ is the vessel fuel consumption per unit distance at the speed ,  is the designed speed, and ℱ is the vessel fuel consumption per hour at the design speed.
Based on the above description, the MBAP follows the following assumptions.
(a) The vessel berthing time window is composed of a strict start time window and a soft departure time window; that is, the vessel must start operations after the expected start time but can leave after the expected finish time, the latter of which will incur a delay penalty.(b) The potential speed of vessels is discretized over a given range [27].(c) Berth depth and width can accommodate all vessels.The relationship among multiple BAPs is established through vessel sailing, and the VSO is further introduced.In the MBAP, the vessel speed is directly related to the fuel consumption and therefore becomes an important part of the total cost.In this paper, the relationship between ship speed and fuel consumption is calculated according to the formula (1) proposed by Venturini et al. [27], where l s is the vessel fuel consumption per unit distance at the speed s, s d is the designed speed, and F is the vessel fuel consumption per hour at the design speed.
Based on the above description, the MBAP follows the following assumptions.
(a) The vessel berthing time window is composed of a strict start time window and a soft departure time window; that is, the vessel must start operations after the expected start time but can leave after the expected finish time, the latter of which will incur a delay penalty.(b) The potential speed of vessels is discretized over a given range [27].(c) Berth depth and width can accommodate all vessels.(d) The berthing safety margin is included in all berths.

Model Formulation
An MILP model for the MBAP is presented in this section, which is developed from the MIP model proposed by Venturini et al. [27].In particular, the space-time relationship constraints of vessels are modified according to Türko gulları et al. [31] and the planning period constraint is introduced.

Notations
The sets, parameters, and decision variables involved in the proposed MILP model for the MBAP are described in Table 1.Binary variable, which equals 1 if vessel j starts handling after vessel i (i, j ∈ V, i ̸ = j) departs in port p ∈ P i P j and 0 otherwise.
Binary variable, which equals 1 if vessel j berths at the right berth of vessel i (i, j ∈ V, i ̸ = j) in port p ∈ P i P j and 0 otherwise

Proposed Optimization Model
The MILP model for the MBAP is represented as follows, where Z is the objective value of the problem.
∑ b∈B p φ vpb = 1, ∀v ∈ V, p ∈ P v (7) The objective function (2) attempts to minimize the sum of vessel idle cost, delay cost, sailing cost, and handling cost for all vessels, where "p, p ′ ∈ P v {p ≺ p ′ }" means vessel v sails to port p ′ after visiting port p [28].Constraint (3) ensures that the vessel starts working after arrival at the port, and constraint (4) further limits it not earlier than the expected start time.Constraints ( 5) and (6) determine the delay time of the vessel, which is non-negative.Constraint (7) stipulates that the vessel is berthed in a certain berth in the visiting port.Constraints ( 8) to (10) prevent conflicts of any two vessels in the same visiting port.Constraints (11) and (12) mandate that a vessel sails to the next visiting port immediately after completing the work.Constraint (13) is the planning time limit, and constraint (14) guarantees that the vessel sails at a certain speed between any two visiting ports.Constraints (15) to (18) set the range of the decision variables.

Solution Method
The VNS is a local search meta-heuristic, which was proposed by Mladenović and Hansen [33] and further explained by Hansen and Mladenović [34].The main ideas of the VNS include shaking and local search (also known as variable neighborhood descent, or VND), and the process is shown as follows.
Step 1: Initialization.Select the set of neighborhood structures N k , k = 1, . . ., k max , to be used in the search; find an initial solution x; and set a stopping condition.
Step 3: Shaking.Generate a point x ′ at random from the k th neighborhood of x Step 4: Local search.Apply some local search method with x ′ as the initial solution; denote x ′′ as the obtained local optimum.
Step 5: Move or not.If this local optimum x ′′ is better than the incumbent (currently best known) solution, move to the new optimum ( x ← x ′′ ) and set k ← 1 ; otherwise, set k ← k + 1 .
Step 6: If k ≤ k max , skip to Step 3; otherwise, the stopping condition is reached.
In this section, we propose a VNS algorithm tailored specifically for the MBAP.We introduce a set of neighborhood structures designed to suit the unique characteristics of the MBAP, along with an integrated initial solution generation strategy.

VNS Framework
In the MBAP, the main decision variables of any vessel v ∈ V in each visiting port p ∈ P v include the sailing speed (s vp ), the work start time (ST vp ), and the choice of berth (b vp ).The remaining variables can be calculated according to these.Therefore, the main decision variable set K vp = s vp , ST vp , b vp of vessel v in port p is defined, and the decision information of the solution x of MBAP includes all K vp .
Further, Algorithm 1 provides the VNS framework for solving the MBAP.Specifically, based on the given shaking neighborhood structure N k i (x), i = {1, 2, . . . ,k max }, local search neighborhood structure N l j (x), j = {1, 2, . . . ,l max }, and the generated initial solution x o , the VNS starts iterating until it reaches the termination condition (cumulative I max iterations).In each iteration, the incumbent solution x is performed shaking and local search until no better solution is produced.Meanwhile, the local search is repeated whenever a better local solution is found, and the local search for each N l j (x) is performed t l times.
Algorithm 1 The VNS framework for the MBAP Input: a set of neighborhood structures N k i (x), i = {1, 2, . . . ,k max } for shaking, a set of neighborhood structure N l j (x), j = {1, 2, . . . ,l max } for local search, an initial solution x o .Set the best solution x = x o .For t = 1 to I max do Find a local best solution x ′′ in t l iterations from N l j (x).Move or not:

End while Move or not:
End for Output: Best solution x.
In the following subsections, the initial solution generation strategy and all neighborhood structures applied in the VNS are further described.

Initial Solution Generation
The MBAP is complex, and the random initial solution generation method is not only inefficient but also may obtain an infeasible solution.Therefore, a greedy initial solution generation method based on the characteristics of the MBAP is proposed, which is described as follows.
Step 2: If n = |P| + 1, skip to Step 11; otherwise, define set A pn , which contains vessels whose n th visiting port is p.
Step 3: Randomly select a vessel v in A pn .
Step 4: Determine the sailing speed s vp by equation (19), and the arrival time a vp is further calculated by equation (20).
Step 5: Determine the start time ST vp = max a vp , EST vp .
Step 6: Define set C according to ST vp , which contains vessels that conflict with vessel v. Also, the berths of non-conflicting vessels are recorded in set B.
Step 7: If B = ϕ, delay ST vp to the earliest departure time of the vessels in C, then skip to Step 6; otherwise, determine b vp as the berth with the smallest h vpb in B.
Step 8: Based on the final ST vp , update s vp and a vp by Equations ( 21) and (20).At this point, the initialization of vessel v is complete, and then A pn − {v}.
Step 11: Obtain an initial solution for the MBAP.

Neighborhood Structures
The neighborhood structures of the VNS are mainly applied in shaking and local search steps [34].The former is mainly used to enrich the diversity of the solutions, while the latter attempts to find local optimal solutions.According to the characteristics of the MBAP and the relationship between decision variables, the following neighborhood structures are proposed in this section.

Shaking
For a given MBAP solution, the shaking step disturbs the speed decision information of some vessels, and the specific neighborhood structure is as follows, where N SK is the number of affected vessels.
Speed Cut: Select N SK vessels with the highest average speed and reduce their maximum speed to the lowest.Higher speed will bring high sailing cost and fuel emissions, which are undesirable.However, blindly reducing the speed may destroy the excellent characteristics of the solution.Therefore, this neighborhood structure considers the average speed but only changes the highest speed among them.
Speed Up: Select N SK vessels with low average speed and increase their minimum speed by one level.Although lower speed is beneficial for cost savings and environmental protection, it can lead to a loss of time.Thus, this neighborhood structure appropriately accelerates vessels with lower average speeds.
Random Speed: Select N SK ships and randomly reset a segment of their speed.This neighborhood structure fully implements randomness and greatly enriches the diversity of solutions.
For vessels whose sailing speed is disturbed, the arrival time is updated according to Equation (20), and the work start time is re-decided according to constraints (3) and (4).However, the above operations may result in conflicts between the affected vessels and other vessels.Taking the disturbed vessel v SK in port p ∈ P v SK as an example, Figure 3a shows two potential conflict situations, and Figure 3b illustrates the corresponding adjustment methods.Specifically, if vessel v SK conflicts with its predecessor vessel v C , the commencement of vessel v SK is delayed until the departure time of vessel v C .Otherwise, the vessel v C is adjusted accordingly.After that, the speed and arrival time of the adjusted vessel are updated according to Equations ( 21) and (20), respectively.other vessels.Taking the disturbed vessel  in port  ∈  as an example, Figure 3a shows two potential conflict situations, and Figure 3b illustrates the corresponding adjustment methods.Specifically, if vessel  conflicts with its predecessor vessel  , the commencement of vessel  is delayed until the departure time of vessel  .Otherwise, the vessel  is adjusted accordingly.After that, the speed and arrival time of the adjusted vessel are updated according to Equations ( 21) and ( 20), respectively.

Local Search
For the solutions generated after each shaking, the local search steps are mainly to find better decision information for vessels with different characteristics by applying multiple neighborhood structures.In particular, each neighborhood structure follows the same destroy-repair process as follows, which improves the destroy-repair neighborhood structure proposed by Wagner and Mönch et al. [35].
Step 1: According to the criteria set by the neighborhood structure, select  vessels and their decision information will be reset.
Step 2: Destroy.Removes information for  vessels from the given solution.
Step 3: Repair.According to Steps 3 to 8 in Section 4.2, re-determine the decision information for these  vessels in visiting ports.Taking the decision-making process of a vessel  ∈  in visiting port  ∈  with three berths as an example, Figure 4 more intuitively shows the repair process.Specifically, as shown in Figure 4a, under the original  , vessel  will have conflicts at each berth in port .Since vessel 1 is the first of the potentially conflicting vessels to depart, the  is delayed until the departure time of vessel 1.Further, as shown in Figure 4b, the vessel  finds a feasible berth 1.As can be seen from the above, the difference between the different neighborhood structures lies in the selection criteria for the destroy-repaired vessels, which are described in detail below.

Local Search
For the solutions generated after each shaking, the local search steps are mainly to find better decision information for vessels with different characteristics by applying multiple neighborhood structures.In particular, each neighborhood structure follows the same destroy-repair process as follows, which improves the destroy-repair neighborhood structure proposed by Wagner and Mönch et al. [35].
Step 1: According to the criteria set by the neighborhood structure, select N DR vessels and their decision information will be reset.
Step 2: Destroy.Removes information for N DR vessels from the given solution.
Step 3: Repair.According to Steps 3 to 8 in Section 4.2, re-determine the decision information for these N DR vessels in visiting ports.Taking the decision-making process of a vessel v ∈ N DR in visiting port p ∈ P v with three berths as an example, Figure 4 more intuitively shows the repair process.Specifically, as shown in Figure 4a, under the original ST vp , vessel v will have conflicts at each berth in port p.Since vessel 1 is the first of the potentially conflicting vessels to depart, the ST vp is delayed until the departure time of vessel 1.Further, as shown in Figure 4b, the vessel v finds a feasible berth 1.
shows two potential conflict situations, and Figure 3b illustrates the correspondin justment methods.Specifically, if vessel  conflicts with its predecessor vessel  commencement of vessel  is delayed until the departure time of vessel  .Other the vessel  is adjusted accordingly.After that, the speed and arrival time of th justed vessel are updated according to Equations ( 21) and ( 20), respectively.

Local Search
For the solutions generated after each shaking, the local search steps are mai find better decision information for vessels with different characteristics by applying tiple neighborhood structures.In particular, each neighborhood structure follow same destroy-repair process as follows, which improves the destroy-repair neig hood structure proposed by Wagner and Mönch et al. [35].
Step 1: According to the criteria set by the neighborhood structure, select  sels and their decision information will be reset.
Step 2: Destroy.Removes information for  vessels from the given solution.
Step 3: Repair.According to Steps 3 to 8 in Section 4.2, re-determine the decisi formation for these  vessels in visiting ports.Taking the decision-making proc a vessel  ∈  in visiting port  ∈  with three berths as an example, Figure 4 intuitively shows the repair process.Specifically, as shown in Figure 4a, under the or  , vessel  will have conflicts at each berth in port .Since vessel 1 is the first potentially conflicting vessels to depart, the  is delayed until the departure ti vessel 1.Further, as shown in Figure 4b, the vessel  finds a feasible berth 1.As can be seen from the above, the difference between the different neighbor structures lies in the selection criteria for the destroy-repaired vessels, which ar scribed in detail below.As can be seen from the above, the difference between the different neighborhood structures lies in the selection criteria for the destroy-repaired vessels, which are described in detail below.
Delay Destroy-Repair: This neighborhood structure removes vessels that have high delay cost.Delays are one of the least desired situations in maritime transportation, and there is a possibility of delays spreading in MBAP [28].
Berth Destroy-Repair: This neighborhood structure considers vessels that have been handling in ports for too long.In addition to the handling time of a vessel in the port being directly related to its own loading and unloading volume, another important factor is the choice of berth.In this paper, the influence of different berths is considered by setting the handling time of different berths, among which the best berth has the shortest handling time.The criterion of this neighborhood structure is based on the relative deviation between the actual handling time and the minimum handling time.
Idle Destroy-Repair: This neighborhood structure gives priority to the selection of vessels with waiting times that are too long.The existence of waiting time for vessels not only generates additional cost but also reflects unreasonable vessel speed optimization.
Random Destroy-Repair: This neighborhood structure is performed at the end of the local search to further enrich the diversity of solutions.Compared with the Random Speed neighborhood structure in the shaking step, this neighborhood structure may apply a better solution.

Numerical Results
The VNS algorithm for the MBAP presented above is coded in C++, while the proposed MILP is solved by the commercial solver CPLEX 12.8.0.Since the existing relevant research instances [27,28] cannot be obtained.In this paper, based on the real port data [5] and the instances proposed by Agra and Rodrigues [32], two sets of MBAP instances are generated respectively for numerical experiments.Specifically, each BAP instance proposed by Cordeau et al. [5] and Agra and Rodrigues [32] is taken as the data for a port, and then a complete MBAP instance is generated by combining the data for multiple ports.
The first set of MBAP instances (A1) based on Cordeau et al. [5] is divided into three scales based on the number of vessels and further refined into six categories based on the number of ports.Moreover, the number of port berths and the generation range of distance between ports are also different for each scale instance, as shown in Table 2.In addition, the handling time of vessels at each berth is derived from the Cordeau et al. [5].Since Cordeau et al. [5] considered the unavailable berths (handling time is 0) that are not considered in this paper, the handling time of these berths in the instances of A1 are the minimum value of the handling time of other berths.Finally, three instances of each type are generated, and each instance in A1 is named according to the format sCale|P|_k, where sCale = {S, M, L}, |P| is the number of vessels, and k = {1, 2, 3} is the serial number.The second set of MBAP instances (A2) is based on the instances proposed by Agra and Rodrigues [32], which are derived from Rodrigues and Agra [36] with the number of vessels covered {6, 7, . .., 15}.In this section, the instances in A2 consider twice the number of vessels, i.e., the number of vessels covering {12, 14, 16, 18, 20, 22, 24, 26, 28, 30}.The number of ports is set to three in all instances, and the distance between ports is varying in [400, 600].Other parameters are generated according to the methods proposed by Agra and Rodrigues [32] and Rodrigues and Agra [36].Specifically, Rodrigues and Agra [36] set the number of berth sections as 34, and the vessel occupancy sections vary in {5, 6, 7}.Based on this, the number of port berths in A2 takes a random value in {5, 6, 7}.Meanwhile, according to the deterministic handling time of vessels proposed by Agra and Rodrigues [32], the disturbance obeying U(1.0, 1.5) is further applied to generate the vessel handling time at each berth in port required in this paper.Finally, the instances in A2 are named R|V|, where |V| is the number of vessels.More details are available in Agra and Rodrigues [32] and Rodrigues and Agra [36].
For the instances in A1 and A2, the port visiting sequence of each vessel is randomly generated according to the number of ports.Meanwhile, the time window of vessels in each visiting port is generated according to (22) and (23), where s is a random speed and ζ is the disturbance coefficient following U(−0.3,0.3).Further, according to Venturini et al. [27], the vessel speed information is discretized into 11 levels with a spacing of 0.5 knots from 14 to 19 knots, and the vessel design speed is set to 19 knots.In addition, the oil price O = USD 250/ton, the idle cost I = USD 200/h, the delay cost E = USD 300/h, the handling cost H = USD 200/h, and the planning period T is set to 1 week.
Based on the generated MBAP instance sets, the performance comparison experiment for the VNS and CPLEX, the effect test for the VNS neighborhood structures, and the sensitivity analysis for the oil price are performed.All numerical experiments below are performed on a Windows PC with Intel (R) Core (TM) i5-8500 CPU and 8.00 GB RAM.For the VNS, the maximum number of iterations (I max ) is set to 1000, and the number of iterations for each neighborhood structure in local search (t l ) is set to 100.Meanwhile, the number of disturbed vessels in the shaking step (|N SK |) is set to |V|/4, and the number of destroy-repaired vessels in the local search step (|N DR |) is set to |V|/2 .Each run records the objective value of the output solution (Z V ) and the solution time.The commercial solver CPLEX limits the maximum solution time to 1 h and records the objective value of the returned solution (Z C ), the gap from the optimal solution, and the return time.

Experimental Comparisons
All the generated MBAP instances are solved by VNS and CPLEX, where the VNS is run 20 times and the objective value of the best solution (Z V B ), the average objective value (Z V A ), and the average solution time are recorded.Detailed comparison results are shown in Tables 3 and 4, where the relative deviations (R B and R A ) between the VNS best solution and the average solution and the CPLEX returned solution are, respectively, according to Equation (24).As shown in the results in the A1 instance set in Table 3, for the results of the CPLEX, the optimal solution can be returned within 1 s when solving the small-scale instances, but it cannot be solved within 1 h when the instance scale is enlarged.For the VNS, when solving small-scale instances, the same optimal solution as CPLEX can be obtained within 10 s, but for some small-scale instances (S2_01 and S3_02), the optimal solution cannot be obtained stably.However, the performance advantages of the VNS become obvious in medium-and large-scale instances, where they can not only obtain a better solution than the CPLEX in 1 min but also further improve the stability of the solution.Therefore, the proposed VNS algorithm for solving the MBAP can meet the solving requirements for small-scale instances and has better solving performance than the CPLEX in medium-and large-scale instances.
The results for the A2 instance set shown in Table 4 show the change in solution performance of VNS compared with CPLEX as the number of vessels increases.Specifically, for the R12, R14, and R16 instances, the VNS is able to return the optimal solution (gap is 0) obtained by CPLEX in about 5 s, and the average solution of the 20 runs deviates from the optimal solution by a maximum of 0.06% and a minimum of 0. Further, for the R18, R20, and R22 instances, CPLEX takes a long time to obtain the optimal solution, with a maximum time of 2895.12 s.However, the VNS can still return the optimal solution within 30 s, although the average solution has a maximum relative deviation of 1.41%.As the number of vessels expands further, the advantages of the VNS become significant.For the R26, R28, and R30 instances, the CPLEX can no longer return the optimal solution within 3600 s, while the VNS can obtain a better solution in 1 min, and the average solution is also better than the return solution of the CPLEX.Therefore, the instance set A2 again shows that the proposed VNS can efficiently and stably solve the large-scale MBAP instances.
Further, the effect of the proposed neighborhood structures in the VNS is tested by comparing the change in solution after removing a neighborhood structure.Specifically, first, a neighborhood structure is removed, and then the VNS without it is repeated 20 times to record the best and average solutions and further calculate their relative deviations from the results in Table 3 by Equation (24).The first instance of each medium and large scale in A1 is selected in the experiment.The results are shown in Figures 5 and 6, where each neighborhood structure is represented by its initial letter combination, and "-" means to remove the neighborhood structure.
As shown in Figures 5 and 6, the neighborhood structures applied in shaking and local search play different roles in the VNS optimization process.For the shaking step, the Random Speed (RS) neighborhood structure plays the most significant role, followed by the Speed Cut (SC) and the Speed Up (SU).This is because while most vessels tend to sail at the desired low speed, individual vessels sailing at high speed will bring greater overall benefits.According to the results of local search neighborhood structures, the vessel delay has the greatest influence on the MBAP, followed by vessel waiting and berth selection.
first, a neighborhood structure is removed, and then the VNS without it is repeated 20 times to record the best and average solutions and further calculate their relative deviations from the results in Table 3 by Equation (24).The first instance of each medium and large scale in A1 is selected in the experiment.The results are shown in Figures 5 and 6, where each neighborhood structure is represented by its initial letter combination, and "-" means to remove the neighborhood structure.As shown in Figures 5 and 6, the neighborhood structures applied in shaking and local search play different roles in the VNS optimization process.For the shaking step, the Random Speed (RS) neighborhood structure plays the most significant role, followed by the Speed Cut (SC) and the Speed Up (SU).This is because while most vessels tend to sail at the desired low speed, individual vessels sailing at high speed will bring greater overall benefits.According to the results of local search neighborhood structures, the vessel delay has the greatest influence on the MBAP, followed by vessel waiting and berth selection.
Therefore, in the actual scheduling process, special attention should be paid to those vessels that have delays and exceedingly long wait times.Meanwhile, it is not desirable to blindly reduce the vessel speed only to reduce the fuel consumption, which may bring more serious time costs.Hence, reasonable speed optimization should allow some vessels to sail at higher speeds.

Sensitivity Analysis of Oil Price
In the MBAP, the oil price determines the vessel sailing cost and is an important factor.This section further explores the impact of oil price on the MBAP through sensitivity analysis, where oil prices are set to USD 250, 400, and 600/ton [27].The experiment selects the same instances as in Figure 5. Based on the best solution of the VNS runs for 20 times, the average oil consumption and average sailing speed of the vessels are illustrated in Figure 7. first, a neighborhood structure is removed, and then the VNS without it is repeated 20 times to record the best and average solutions and further calculate their relative deviations from the results in Table 3 by Equation (24).The first instance of each medium and large scale in A1 is selected in the experiment.The results are shown in Figures 5 and 6, where each neighborhood structure is represented by its initial letter combination, and "-" means to remove the neighborhood structure.As shown in Figures 5 and 6, the neighborhood structures applied in shaking and local search play different roles in the VNS optimization process.For the shaking step, the Random Speed (RS) neighborhood structure plays the most significant role, followed by the Speed Cut (SC) and the Speed Up (SU).This is because while most vessels tend to sail at the desired low speed, individual vessels sailing at high speed will bring greater overall benefits.According to the results of local search neighborhood structures, the vessel delay has the greatest influence on the MBAP, followed by vessel waiting and berth selection.
Therefore, in the actual scheduling process, special attention should be paid to those vessels that have delays and exceedingly long wait times.Meanwhile, it is not desirable to blindly reduce the vessel speed only to reduce the fuel consumption, which may bring more serious time costs.Hence, reasonable speed optimization should allow some vessels to sail at higher speeds.

Sensitivity Analysis of Oil Price
In the MBAP, the oil price determines the vessel sailing cost and is an important factor.This section further explores the impact of oil price on the MBAP through sensitivity analysis, where oil prices are set to USD 250, 400, and 600/ton [27].The experiment selects the same instances as in Figure 5. Based on the best solution of the VNS runs for 20 times, the average oil consumption and average sailing speed of the vessels are illustrated in Figure 7. Therefore, in the actual scheduling process, special attention should be paid to those vessels that have delays and exceedingly long wait times.Meanwhile, it is not desirable to blindly reduce the vessel speed only to reduce the fuel consumption, which may bring more serious time costs.Hence, reasonable speed optimization should allow some vessels to sail at higher speeds.

Sensitivity Analysis of Oil Price
In the MBAP, the oil price determines the vessel sailing cost and is an important factor.This section further explores the impact of oil price on the MBAP through sensitivity analysis, where oil prices are set to USD 250, 400, and 600/ton [27].The experiment selects the same instances as in Figure 5. Based on the best solution of the VNS runs for 20 times, the average oil consumption and average sailing speed of the vessels are illustrated in Figure 7.As shown in Figure 7, for all instances, when the oil price increases, the average fuel consumption and average speed of the vessels will decrease.This result suggests that high oil prices will force the vessels to reduce sailing costs by slowing down, even with the potential loss of time.In reality, the oil price is time varying, so it is extremely important to optimize the vessel speed reasonably according to different oil prices.Specifically, when the oil price is high, shipping companies can reduce fuel consumption by appropriate speed reduction but need to consider the minimum time limit in operation.When the oil price is low, it is possible to increase the timeliness of operations by slightly accelerating the vessel speeds to avoid delays.
In general, in the MBAP, the VSO is not only affected by transportation demand but As shown in Figure 7, for all instances, when the oil price increases, the average fuel consumption and average speed of the vessels will decrease.This result suggests that high oil prices will force the vessels to reduce sailing costs by slowing down, even with the potential loss of time.In reality, the oil price is time varying, so it is extremely important to optimize the vessel speed reasonably according to different oil prices.Specifically, when the oil price is high, shipping companies can reduce fuel consumption by appropriate speed reduction but need to consider the minimum time limit in operation.When the oil price is low, it is possible to increase the timeliness of operations by slightly accelerating the vessel speeds to avoid delays.
In general, in the MBAP, the VSO is not only affected by transportation demand but also by time-varying oil prices.Meanwhile, the solution of the VSO is directly related to the berth allocation plan of each port.Therefore, a suitable VSO scheme is flexible and based on demand, port operations, and real-time oil prices.

An MBAP Case Study
In this section, a case study is conducted based on practical data from a port group in the Yangtze River Basin in central and western China, which is obtained from [37].This case includes four ports, with the parameters depicted in Table 5.In order to adapt to the MBAP studied in this paper, the direct linear distance between every two ports is calculated according to the relative coordinates of each port.Meanwhile, the case includes 20 vessels, and the relevant parameters are shown in Table 6.In order to apply the method to the MBAP, the port visiting sequence P v of each vessel v is randomly generated.Note that due to the limitation of berth length in the original data, ports that could not accommodate vessels are not considered when randomly generating the port visiting sequence for each vessel.The handling time h vpb each berth in each port is calculated by (25), where C v is the total number of containers (from [37]) and k C is the unit container handling time (set to 20 min according to [37]).In addition, due to the small number of port berths in this case, the planning period has been extended to 10 days, and other parameters are generated in the same way as A1 and A2 instance sets.
Constraints ( 3), ( 7)~( 12), ( 14)~ (18).Based on models M1 and M2, the VNS is called for 20 times to solve, and the comparison results between their best solutions and M0 are recorded in Table 7, including average vessel speed, total fuel consumption, and latest vessel completion time.As shown in Table 7, the average vessel speed and total fuel consumption decrease when the vessel service time window is removed.Vessels reduce sailing costs by slowing down, since there is no penalty for late departure, and this also leads to an increase in latest vessel completion time.Meanwhile, due to the planning period constraint, the vessels' deceleration is limited.When the constraint of the planning period is further relaxed, the average vessel speed and total fuel consumption continue to decline, in which the average vessel speed reaches the minimum.At this point, the vessel no longer considers the total completion time, and the latest vessel completion time rises to 273 h.As the constraints of time window and planning period are gradually ignored, the effect of vessel speed optimization gradually disappears, which further leads to a substantial delay in scheduling time.Therefore, it is necessary to integrate vessel speed optimization and consider the vessel service time window and planning period in the MBAP.
In practical scheduling, vessels can be forced to carry out reasonable speed optimization by strictly limiting the planning period and giving delay penalties, which will help alleviate or even solve maritime difficulties such as port congestion and unreliable shipping schedules.

Conclusions
This paper investigates the MBAP, offering insights into potential collaborations between ports and shipping companies.Such partnerships aim to enhance maritime transportation efficiency and address existing challenges.Building upon prior research [27,28], we develop a new mixed-integer linear programming (MILP) model for the MBAP, incorporating constraints from Türko gulları et al. [31] regarding vessel space-time relationships.Additionally, we propose a VNS algorithm to effectively solve the MBAP.
Furthermore, we generated an instance set A1 containing MBAP instances of different sizes using real port data and an MBAP instance set A2 based on BAP instances from existing studies.Numerical experiments performed on the two instance sets demonstrate the efficiency and stability of our proposed VNS algorithm, even for large-scale instances.Moreover, performance tests based on instance set A1 reveal that key factors influencing the MBAP include high-speed sailing, delayed departures, and idle waiting scenarios.
By analyzing the impact of different oil prices on the MBAP, several management insights emerge: (1) Adjusting vessel acceleration can improve the scheduling timeliness during periods of low oil prices.(2) Slowing down during high oil prices can help control overall costs.(3) Shipping companies should consider transportation demand and port operation conditions comprehensively when optimizing vessel speeds.Finally, an MBAP case study demonstrates the importance of integrating VSO in the MBAP as well as the need to consider the vessel service time window and planning period.In conclusion, portshipping cooperation proves to be an effective strategy for addressing current maritime challenges, with optimized scheduling offering significant benefits to all involved parties.Future research could explore integrating additional port equipment, such as quay cranes, trucks, and yards, into the MBAP for collaborative scheduling.Additionally, considering vessel emissions in ports aligns with the focus on green shipping.Furthermore, given the

Figure 3 .
Figure 3. Diagram of the vessel conflict adjustment.

Figure 4 .
Figure 4. Diagram of the repair process.

Figure 3 .
Figure 3. Diagram of the vessel conflict adjustment.

Figure 3 .
Figure 3. Diagram of the vessel conflict adjustment.

Figure 4 .
Figure 4. Diagram of the repair process.

Figure 4 .
Figure 4. Diagram of the repair process.

Figure 5 .
Figure 5. Shaking neighborhood structure effect experiments results.

Figure 6 .
Figure 6.Local search neighborhood structure effect experiment results.

Figure 5 .
Figure 5. Shaking neighborhood structure effect experiments results.

Figure 5 .
Figure 5. Shaking neighborhood structure effect experiments results.

Figure 6 .
Figure 6.Local search neighborhood structure effect experiment results.

Figure 6 .
Figure 6.Local search neighborhood structure effect experiment results.

Figure 7 .
Figure 7. Results of solving MBAP by VNS under different oil prices.

Figure 7 .
Figure 7. Results of solving MBAP by VNS under different oil prices.

Table 1 .
Notations and variables.Binary variable, which equals 1 if vessel v ∈ V sails from port p to the next port p ′ (p, p ′ ∈ P v := p ≺ p ′ ) at speed s ∈ S and 0 otherwise a vp Arriving time of vessel v ∈ V in port p ∈ P v φ vpb Binary variable, which equals 1 if vessel v ∈ V in port p ∈ P v is served at berth b ∈ B v , and 0 otherwise ST vp Start time of vessel v ∈ V in port p ∈ P v d vp Departure time of vessel v ∈ V in port p ∈ P v e vp Delay time of vessel v ∈ V in port p ∈ P v

Table 2 .
Parameters of the MBAP instances in A1 with different scales.

Table 3 .
Comparison results of VNS and CPLEX for instances in A1.

Table 4 .
Comparison results of VNS and CPLEX for instances in A2.

Table 7 .
Comparison results of different MBAP models.