A multi-objective fuzzy programming model for port tugboat scheduling based on the Stackelberg game

The overall cargo throughput of Chinese ports still maintains rapid growth. In 2022, the total cargo throughput of ports nationwide reached 15.68 billion tons, which is a year-by-year increase of 0.9%. Against the backdrop of slowing economic growth and accelerating industrial restructuring, China’s port throughput still maintains a steady growth trend¹. The Shanghai International Shipping Center has been built, and the maritime liner connectivity index is firmly ranked first in the world. According to the recent statistics, maritime transportation accounts for over 95% of China’s foreign trade cargo transportation volume².

In trade transactions, ports, as an important component of modern logistics chains, are the starting and ending points of maritime transportation and are also a comprehensive logistics system that should provide safe entry, exit, and berthing for vessels, including green port investments³, various buildings, facilities, machinery, tools, and personnel⁴. The operations at port waters (referring to waters surrounding berths, waterways, and anchorages) cannot be separated from tugboats. As a tool, tugboats assist ships entering and leaving the port by providing necessary berthing, shifting, and piloting services. They are an important port resource for achieving port functions. However, the availability of tugboats and the ability of each tugboat to provide services to the arriving vessels are limited. Therefore, in order to improve the efficiency of tugboat operations, tasks are expected to be dispatched through scientific methods and analytical tools. Tugboat scheduling is an important issue that urgently needs to be addressed to improve the efficiency and competitiveness of port operations.

In the process of tugboat operations, objective entities, such as vessels, shipping agents, port dispatchers, pilot companies, barge operators, pilots, tugboats, etc., are involved, among which the port dispatchers and tugboat companies are the main entities in formulating tugboat dispatch plans for port operations⁵. At present, in the actual planning operation of tugboat scheduling at the Guangdong Province ports, the port scheduling party processes large vessels into each independent operation based on customer requests and determines the tugboat power and quantity matching plan required for each operation according to the “Guangdong Province Port Ship Berthing, Departure, Pilot or Transfer Tug Ship Allocation Standard”⁶, and then the key information is entered into the scheduling management system. The tugboat operator prepares a specific tugboat operation plan based on this information and directs the tugboat deployment. The initial matching plan formulated by the port scheduling party needs to be adjusted based on the number of tugs owned by the tugboat operator and the actual operational conditions of a given day. As such, there are direct interactions between the tugboat operator and the port scheduling party.

Under the guidance of the tugboat and task matching scheme provided by the port dispatching party, the tugboat operator often uses the shortest distance tugboat principle, the first available tugboat principle, and the principle of fairness in the task quantity of each tugboat to assign tugboats. Among them, the principle of shortest distance tugboat operation (TSD) refers to first selecting the tugboat with the shortest distance (referring to the distance between the tugboat’s ___location after completing the previous task and the current task ___location or the distance between the tugboat’s berthing base and the current task ___location) as the scheduling plan according to the matching rules between the tugboat and the vessel and then dispatching that tugboat to provide services to the vessels to be served. The First Available Tug Principle (FAT) refers to first selection of tug schemes that meet the requirements of the vessel based on the matching rules between tugs and vessels. From the available tug schemes, the tug scheme that completes the previous operation as quickly as possible and is idle is selected as the adopted scheme, allowing it to provide services for the vessel. The principle of shortest distance tugboat operation is very similar to the principle of the first available tugboat in some cases, with the difference being that the focus of the TSD principle is on how to quickly operate for the next vessel, while the FAT principle considers the overall management of all tugboat operations. The principle of fairness in the task volume of each tugboat (UWAT) is to ensure that the same tugboat driver does not experience excessive fatigue in the same shift and ensures sufficient rest time. Many tugboat operators provide compensation to tugboat drivers based on the number and importance of tasks they participate in. Therefore, to ensure fairness in the salaries of tugboat drivers, tugboat dispatchers allocate tugboat tasks evenly over time, ensuring that the number of tasks each tugboat driver participates in does not differ significantly from other tugboat drivers. Most literature considers the first two principles^6,7, and this research will simultaneously consider the above three principles for the dispatch of tugboats.

After a vessel enters the port for operation, it is possible to estimate the arrival time and berthing time of the vessel, but it is difficult to predict the departure time from the berth because the loading and unloading time is difficult to determine, so the departure time may fluctuate. Hence, we can only know the notification time and determine the end time of this departure task. For the situation where this random task arrives, the strategy that will be explored in this research is to maximize the buffer time of this dynamic task and to facilitate the actual deployment of tugboats. Furthermore, the speed of vessels berthing, disembarking, and loading and unloading cargo mainly depends on the coordination of various buildings, facilities, machinery, tools, and personnel resources at the dock, and has great uncertainty. Moreover, the historical data makes it difficult to predict all operational situations. This article uses fuzzy numbers to describe the operation time of each task. In addition, busy ports need to serve dozens or even hundreds of vessels every day. When the problem scale increases, the model built in this study is difficult to solve and has certain limitations. Therefore, this paper proposes a metaheuristic algorithm to solve it.

In summary, it is necessary to weigh various factors to make scientific and reasonable tugboat scheduling decisions for busy ports. This paper investigates the optimization problem of tugboat scheduling to assist vessels in entering and exiting ports in uncertain environments and establishes a game-theoretic fuzzy model for tugboat scheduling in port operations to solve the trade-off problem between tugboat operators and port dispatchers in the tugboat scheduling process. In addition, this study incorporates three commonly used scheduling rules, task random arrival, and fuzziness into the modeling and develops a seagull optimization algorithm that sets priority encoding, fitness function, and genetic operators based on the characteristics of the tugboat scheduling problem addressed herein. The main purpose of this study is to assist tugboat operators and port dispatchers in making daily real-time decisions in vessel scheduling.

The innovation of this research lies in (1) considering the diversity of tugboat matching schemes to increase the optimality of tugboat scheduling; (2) taking into consideration of actual operational situations to make the model closer to reality; (3) using the principle of fair dispatching for tugboat scheduling; (4) exploring the situations where the task can arrive at any time; and (5) identifying the impact of game-theoretic decision-making on tugboat scheduling plans for the first time. The rest of this study is organized as follows. Section “Related work” reviews and summarizes the relevant work in the area of tugboat scheduling. Section “Problem description and modeling” provides the background of the problem and describes the composition of the model. Section Solution method outlines the main steps of the proposed solution algorithm. Section “Numerical experiments” formally presents the numerical experiments, and section “Conclusion” summarizes the results of this research.

Some scholars have studied the optimization of tugboat resources from the perspective of classic allocation problems^8,9,10,11,12. In the field of maritime transportation research, scholars have also conducted extensive research on the decision-making problems of berths or terminals, but there are few studies related to the optimization of tugboat resource allocation^12,13,14. Early research often used heuristic algorithms to solve tugboat resource optimization problems. As an example, Wang and Meng¹⁵ developed a heuristic algorithm that combined a genetic algorithm and ant colony optimization for container terminal task resource allocation problems. In the first stage, the genetic algorithm used strings to represent the chromosomes of the allocation plan and found the optimal allocation through self-learning. In the second stage, an improved ant colony optimization algorithm was introduced to optimize scheduling tasks based on the genetic algorithm allocation plans. Finally, the performance of the algorithm was tested using a tugboat allocation problem, and satisfactory results were obtained.

In 2012, Wang et al.^16,17,18 proposed a tugboat allocation optimization model from the perspective of matching the captain with the horsepower and quantity of tugboats and applied an improved trust-based ant colony optimization method to obtain the optimal scheduling plan. In the same year, they added time logic constraints to the model and established a mixed integer programming model that combined scheduling rules. They then provided a method to determine the optimal solution of the problem and analyzed the impact of the number of tugs and service capacity on ship turnover time through numerical experiments. However, the models in these three papers could not be directly solved until 2018, when Zhen et al.¹⁹ studied the tugboat allocation problem in seaports located at estuaries and proposed a mixed integer programming model to optimize the allocation scheme of barges to tugboats. They also proposed an accurate solution approach based on the branch-and-bound method to solve the proposed model. Wu et al.²⁰ studied the evacuation planning problem of barges that may arise in reclamation projects. The authors made decisions from strategic, tactical, and operational perspectives on the allocation and scheduling of tugs for a team of heterogeneous barges working in offshore land reclamation areas that need to evacuate to coastal shelters before a storm arrives and proposed a customized heuristic method.

In terms of incorporating time slot resource allocation into the allocation of tugboats, Omar et al.⁵ considered the berth allocation problem with all required resources (tugs, pilots, pilot vessels, and mooring teams), safety considerations (within and between vessels), and berth ___location availability. They proposed a mixed integer programming model and precise solution method based on the constraint separation strategy for port entry. It was found that the decision on the start and end times of departure and relocation showed significant improvements in the management of port resources, especially during congestion periods. Jia et al.²¹ determined the schedule for tugs to serve vessels and the schedule for vessels to berth or depart during berthing and turning operations. A case study was considered based on the operational data from a container port in Shanghai to demonstrate that the proposed iterative algorithm combining Lagrange relaxation and Benders decomposition outperformed the benchmark algorithm in terms of solution quality. For scheduling operations in canals, Petris et al.²² took the medium-sized Port of Venice in Italy as an example and proposed four heuristic algorithms. The developed solution algorithms were found to be efficient and could provide high-quality solutions for most of the instances in a reasonable computational time.

Liu et al.²³ proposed a two-stage stochastic mixed integer linear programming model for harbor berth and channel planning, considering the limited availability of berths and channels. The aim was to minimize the expected total weighted completion time of vessels in the presence of uncertainty in vessel arrival time and vessel processing duration. The first stage determined the allocation of berths for vessels under multi-scenario conditions. In the second stage, channel planning, including channel selection, tugboat allocation, and vessel sequencing, was determined in the presence of uncertainty. The authors proposed a stage decomposition method and decomposition-based heuristic algorithm to solve the problem. Kasm et al.²⁴ established a mixed integer programming to simulate vessel scheduling under channel constraints and different tugboat allocation strategies and achieve better resource utilization. Jiang et al.²⁵ established a mathematical optimization model from the perspective of time allocation, considering constraints, such as restricted channels and complex navigation rules for scheduling resources (including tugboats), and proposed an elite selection genetic algorithm encoded with vessel service sequences to solve the optimization model. It was pointed out that the modeling from the perspective of allocation issues required adding variables and increasing the complexity of the model.

Another group of scholars study tugboat scheduling optimization problems from the perspective of classic vehicle routing problems. For instance, Wei et al.⁷ formed a point line network of berths, tugboat bases, and artificial waypoints in defined port waters, establishing a mixed integer linear programming model for tugboat assisted container vessel berthing, moving, and unberthing processes, and developing six effective inequality groups that were incorporated into the standard branch-and-bound algorithm. However, this study simply considered that the tugboat base is a fixed starting and ending point, without considering that tugboats can dock at any tugboat base. In terms of rivers and seaports, Zhu et al.²⁶ designed a variable neighborhood search algorithm to address the social responsibility of vesselping companies in reducing greenhouse gas emissions by arranging tugs and barges reasonably and determining the transportation route for barge transportation. The calculation results show that compared with the scheduling rules, the model reduced carbon emissions by approximately 46.93%. In addition, consideration of barge transportation in the model could reduce carbon emissions by approximately 10.46%. In addition, the variable neighborhood search algorithm demonstrated an optimality difference of approximately 0.29% in a short period. Hao et al.²⁷ considered the joint scheduling of barges and tugs at river sea transfer ports with tidal influence channels, established an integer programming model based on a state spatiotemporal network to minimize the time required for bulk cargo transfer, and designed a customized variable neighborhood search algorithm. Wang et al.²⁸ studied the tugboat scheduling problem considering multiple waypoints and multiple service modes, established a mixed integer linear programming model that can better utilize limited tugboats, and developed an efficient adaptive large neighborhood search algorithm. Li et al.²⁹ proposed a new fuzzy programming optimization model for tugboat scheduling, considering multiple berthing bases, time windows, and operational uncertainties.

The remaining studies on tugboat scheduling optimization addressed the problem from the perspective of a flow shop. Kang et al.³⁰ considered the uncertainty of container vessel arrival and tugboat process times in large container ports. The uncertain vessel arrival and tugboat operation process times were represented as a limited set of discrete scenarios, and uncertainty was addressed by integrating active and reactive scheduling strategies. The authors established a mixed integer linear programming model for the tugboat scheduling problem and designed a self-organizing algorithm to generate tugboat chains and solve large-scale problems. Although this study considered the uncertainty of vessel arrival and tugboat operation processes, the proposed approach would increase the difficulty of solving the model, and it is challenging to estimate the start time of an unpredictable operation for the vessel. Therefore, the present study adopted fuzzy numbers to replace the uncertainty of tugboat operation time and set a goal to reserve as much buffer time as possible for the operation to start. Taking into account natural factors, Zhong et al.⁶ studied the green tugboat scheduling problem for tidal ports with multiple port areas in long waterways. Considering the factor of the tidal time window, the schedule of tugboat service vessels was determined based on task sequence numbering and captured the trade-off between time and fuel consumption goals. The tidal time window was only a simple upper threshold value. However, modeling from the perspective of the flow workshop makes it difficult to address the issue of vessels not entering the port in the predetermined order.

Tug scheduling is a complex multi-party decision-making process, which has been studied by some scholars from the perspective of multi-stage decision-making. Musus et al.³¹ studied the strategic issue of locating and allocating tugs along the Norwegian coast to optimize sea readiness in emergency rescue and minimize the sum of the costs of public-operated tugs in emergency towing services and the expected penalty costs incurred due to inadequate preparation. A two-stage stochastic programming model was established, with the first stage positioning tugs to meet nominal coverage requirements, and the second stage deploying positioned tugs to assist vessels in distress. However, tugboats were allocated by port dispatchers and dispatched by barge operators, which can be described with a two-party game-theoretic process, and such a methodological approach is lacking in the existing literature.

In the process of tugboat dispatching, decision makers need to pay attention to fairness to balance the workload of tugboat drivers. Orgut et al.³² proposed by definition that each region within the service area should receive fairness constraints that are fully proportional to the needs of that region. The study also proposed mathematical models to promote fair and effective distribution of donated food by food banks among populations at risk of hunger. Se and Oyz³³ studied the optimal distribution of influenza vaccines in heterogeneous populations composed of multiple subpopulations and proposed a mathematical model based on fairness constraints to help public health authorities incorporating fairness in vaccine allocation decisions. However, to the best of the authors’ knowledge, there are no significant research efforts on fair dispatch of tugboats.

In summary, the allocation and scheduling of tugboat resources is a relatively new topic. Although many researchers have devoted their studies to research certain aspects of the tugboat scheduling problem, there is little research on tugboat allocation and scheduling from the perspective of bipartite games, considering three commonly used scheduling rules, task random arrival, and fuzziness in relevant operational features. To this end, this study proposes a new mathematical model for optimizing the allocation and scheduling of tugboats from the perspective of a Stackelberg-based game-theoretic model, taking into account the well-known scheduling rules, task random arrival, allocation fairness, and fuzziness in operation times of tasks.

Problem definition

Firstly, we clarify the main four definitions used in this study:

1. (1)

   Task: Before the tugboat starts operations, the port scheduling party arranges the operation plan for the day through customer requests to receive the vessel’s operation requirements. As shown in Fig. 1, when a vessel enters the port, there are four processes on the port waters, including entering the channel, berthing, unberthing, and exiting the channel. Thus, vessels that require tugboat assistance in operation are generally divided into four tasks before entering the port. The service content, start time, and start and end locations of individual vessel operations can be known in advance, and the minimum number and minimum power of tugboats required can be determined based on the type of vessel. This study defines these four aspects of information as tasks. Specifically, when the port scheduling party makes berth plans, it is possible to predict the start time of vessel berthing, and the start time of entering the waterway can also be predicted through the length of the waterway. However, due to the difficulty in determining the loading and unloading time, the start time of departure is difficult to predict, and only the time and deadline for the task can be predicted. To this end, there are two types of tasks defined in this article: static tasks (i.e., tasks 1 and 2) and dynamic tasks (i.e., tasks 3 and 4). In addition, the task content also includes top flow, downstream turning anchoring, and anchor turning operations.

   

   Entry and exit processes of a single vessel.

   Full size image

2. (2)

   Tugboat operation: When working as a tugboat in a port and completing a given task, if there is no assignment for the next task available, the tugboat needs to proceed to a certain ___location within the port serving as a tugboat mooring base. If there is an assignment to handle the next task available, the tugboat needs to immediately sail to the ___location of the next task.

3. (3)

   Initial state: It refers to the base where the tugboat docked after the previous day’s task was completed.

4. (4)

   Scheduling plan: This refers to arranging tugboat services for each task based on the initial status of each tugboat and the set of tasks for a given day.

Schematic diagram of the Stackelberg game-theoretic model structure.

Full size image

The tugboat scheduling problem studied in this paper is the allocation and scheduling problem between multiple berthing bases, multiple tugboats, and multiple tasks. The tugboat can dock at any berthing base but must start from a certain base (initial state), complete a task, return to a certain base, or continue to complete subsequent tasks. As shown in Fig. 2, considering the scarcity of tugboat resources in the port, the port scheduling party allocates a tugboat matching plan with the minimum number of tugboats for each task according to the “Standard for the Allocation of Tug Ships for Berthing, Leaving, Piloting, or Moving”, and then feeds back the plan to the tugboat operator’s scheduling plan. The tugboat operator can make decisions based on the status of the tugboat and the principle of shortest distance tugboat. The first available tugboat principle and the principle of fairness in the task volume of each tugboat are limited and implemented for dispatching from the ___location where the previous task was completed to the mooring base or dispatching a tugboat to the current task that needs to be executed. Hence, the decisions of both parties interact and form a game-theoretic relationship.

Tugboat services are charged based on the number of tugboat operations, and the total amount of fees is also determined after the task is determined. The quality of the scheduling plan depends on the level of operation costs, which are mainly based on fuel costs. Therefore, the goal of this study is to minimize the total fuel consumption, which consists of the driving fuel consumption (calculated by time) from the base to the start of the task, from the end of the task to the base, and from the end of the preceding task to the start of the succeeding task, as well as the operating fuel consumption (calculated by time) to complete the task. For dynamic tasks, the tugboat operator must also reserve enough time to dispatch suitable tugboats to serve them, while also taking into account the total tugboat operation time of the day, which could be affected by uncertainties. In the actual operational environment, due to the interference of uncertain factors, dispatchers cannot accurately grasp the port environment and operational changes. Therefore, the operation time of the task set in this study is uncertain, and the changes in parameters are represented by symmetric trapezoidal fuzzy numbers (as discussed more in detail in the following sections of the manuscript).

Based on the above description, the main decision problem studied in this research is the game-theoretic optimization process between the allocation of tugs by the port dispatcher and the dispatch of tugs by the tugboat operator. In the environment where dynamic tasks and fuzzy numbers coexist, the goal is to utilize the existing tug resources of tugboat operators to more economically and timely complete the operations involved in vessel entry and exit from the port.

Model establishment

Assumptions

Based on the problem description, this study proposes the following assumptions: (1) the scheduling plan period is 1 day; (2) Each tugboat docks at the base before all tasks begin on that day; (3) The difference between external and internal berths of the base is not considered; (4) The standby time of the tugboat is short, and the time it takes to dock and leave the base berth is negligible; (5) The tugboat has a fixed speed (which is equal to the average speed); and (6) There is no limit to the number of tugboats that can be docked at each base.

Symbol description

For the convenience of establishing the model, this article defines the following symbols:

1. (1)

   Set:

   I is the set of tasks,\(i,j \in I=\{ 1,2, \ldots NI\}\),\(NI\)is total task amount;\(I^{\prime}\)is the set of dynamic tasks,\(i,j \in I^{\prime},I^{\prime} \subseteq I\); J is the union set of 0 tasks and other tasks,\(i,j \in J=\{ 0,1,2, \ldots ,NI\}\), 0 represents the task at 0 o’clock before the start of the day, which means that the task is not executed and a tugboat is only docked at the base; K is the set of tugboats,\(k \in K=\{ 1,2, \ldots NK\}\),\(NK\)is the total number of tugboats; Lis the set of bases,\(l \in L=\{ 1,2, \ldots ,NL\}\),\(NL\)is the total amount of the bases; R is the union set of 0 base and other berthing bases,\(l \in R=\{ 0,1,2, \ldots ,NL\}\), 0 means no docking at the base。

2. (2)

   Parameters:

   \({\tilde {w}_i}\)is the fuzzy service time for task i, \(i \in J\), measured in min;\(c_{{li}}^{1}\)is the distance from base l to the starting point of task i,\(l \in L,i \in J\), measured in km; \(c_{{li}}^{2}\)is the distance from base l to the end of task i,\(l \in L,i \in J\), measured in km;\({B_{kl}}\)is the position of tugboat kat base lat the end of the previous day’s mission,\(k \in K,l \in R\), which is a 0–1 matrix, if the tugboat stops at the base position, it is 1, otherwise it is 0; Mis a sufficiently large value;\({D_i}\)is the number of tugboat requirement for task i,\(i \in I\);\(vc\)is the average speed of the tugboat, measured in km/h;\({P_k}\)is the power of tugboat k, \(k \in K\), measured in HP;\({m_i}\)is the required tugboat horsepower for task i,\(i \in I\), measured in HP;\({\gamma _k}\)is the fuel consumption per unit time of operation,\(k \in K\), measured in kg/min;\(C{o_k}\)is the unit fuel consumption during navigation,\(k \in K\), measured in kg/min;\({Q_i}\)is the number of tugboats with specified power requirements for task i,\(i \in I\);\(c_{{ij}}^{3}\)is the distance from the end point of task ito the starting point of task j,\(i,j \in I\), measured in km;\({A_i}\)is the estimated start time of task i,\(i \in I\), specifically, when\(i \in I^{\prime}\), this symbol represents the time when the notification for task iexists, measured in min;\({\tau _i}\)is the maximum waiting time for task i,\(i \in I\), specifically, when\(i \in I^{\prime}\), this symbol represents the deadline period for task i, measured in min;\({\tilde {\pi }_i}\)is the fuzzy maximum delay time of tugboat service for task i,\(i \in I\), measured in min; N is the fairness coefficient, with the value of\(\left[ {0,1} \right]\).

3. (3)

   Upper level decision variables:

   \({x_{klij}}\) is a 0–1 variable, and if the tugboat kcompletes task iand returns to base lto complete task j, it is 1; otherwise, it is 0, \(k \in K,l \in R,i,j \in J\);\({s_{ki}}\)is the start time of task iserved by tugboat k,\(k \in K,i \in J\);\({y_{ki}}\)is the service time of tugboat kfor task i,\(k \in K,i \in J\).

4. (4)

   Lower level decision variables:

   \({x^{\prime}_{ki}}\)is a 0–1 variable, if tugboat kis assigned to task i, it is 1; otherwise, it is 0,\(k \in K,i \in I\).

Upper level programming model

The upper-level decision maker is the tugboat operator, whose goal is to reasonably dispatch tugboats to carry out the tasks, to minimize the daily total fuel consumption of their company, and to refer to the requirements of the port’s tugboat matching plan. The upper-level planning model M1 has three objective functions, namely minimizing the total daily fuel consumption of tugboat operations, maximizing the total buffer time of dynamic tasks, and minimizing the total completion time:

1. (1)

   Objective function 1 (total fuel consumption for tugboat operations) includes minimizing the sum of fuel consumption when sailing from the base to the start point of the task, from the end point of the task to the base, from the end point of the previous task to the start point of the next task, and the waiting time from the start time of the task to the end time of the task:

   $$\begin{gathered} \hbox{min} \;{f_1}=60\sum\limits_{{k \in K}} {\sum\limits_{{l \in L}} {\sum\limits_{{i \in J}} {\sum\limits_{{j \in I}} {{{C{o_k}c_{{lj}}^{1}{x_{klij}}} \mathord{\left/ {\vphantom {{C{o_k}c_{{lj}}^{1}{x_{klij}}} {vc}}} \right. \kern-0pt} {vc}}} } } } +60\sum\limits_{{k \in K}} {\sum\limits_{{l \in L}} {\sum\limits_{{i \in I}} {\sum\limits_{{j \in J}} {{{C{o_k}c_{{li}}^{2}{x_{klij}}} \mathord{\left/ {\vphantom {{C{o_k}c_{{li}}^{2}{x_{klij}}} {vc}}} \right. \kern-0pt} {vc}}} } } } \hfill \\ \quad \quad \quad \quad +60\sum\limits_{{k \in K}} {\sum\limits_{{l \in \left\{ 0 \right\}}} {\sum\limits_{{i \in I}} {\sum\limits_{{j \in I}} {{{C{o_k}c_{{ij}}^{3}{x_{klij}}} \mathord{\left/ {\vphantom {{C{o_k}c_{{ij}}^{3}{x_{klij}}} {vc}}} \right. \kern-0pt} {vc}}} } } } +\sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {{\gamma _k}{y_{ki}}} } \hfill \\ \end{gathered}$$

   (1)
2. (2)

   Objective function 2 (total buffer time of dynamic tasks) is expressed as maximizing the sum of the difference between the start time of the dynamic task and the end time of the immediately preceding task:

   $$\hbox{max} \;{f_2}=\sum\limits_{{k \in K}} {\sum\limits_{{l \in L}} {\sum\limits_{{i \in I}} {\sum\limits_{{j \in I^{\prime}}} {\left[ {\left( {{A_j}+{\tau _j} - {y_{kj}}} \right){x_{klij}} - \left( {{s_{ki}}+{y_{ki}}} \right){x_{klij}}} \right]} } } }$$

   (2)
3. (3)

   Objective function 3 (total completion time) represents minimizing the maximum task end time:

   $$\hbox{min} \;{f_3}=\mathop {\hbox{max} }\limits_{{k \in K,i \in I}} \left\{ {{s_{ki}}+{y_{ki}}} \right\}$$

   (3)

According to the model assumptions, the decision variables need to meet the following constraints:

Among them, Eq. (4) ensures that each tugboat is positioned at the initial mooring base after completing the tasks from the previous day. Equation (5) assures that a given tugboat is dispatched to the next task after completing the previous task. Equation (6) indicates that the start and execution times of the previous task, the total time to return to a certain base after completing this task, and then depart from this base to reach the starting position of the next task should be earlier than the start operation time of the next task. Equation (7) ensures that the start and execution times of the previous task, and the total time required to directly reach the starting position of the next task should be earlier than the start operation time of the next task. Equation (8) limits the task start time to a certain range. Equation (9) limits the service time to a certain allowable range. Equation (10) imposes constraints on the actual dynamic task end time based on the maximum task waiting time. Equation (11) restricts a given tugboat from stopping at non-zero base at the end of the day. Equation (12) represents a fair constraint on the task quantity for each tugboat, which adjusts the fairness of the task quantity assigned to the tugboat by changing the value of N. Equation (13) represents the integrality constraints for the 0–1 variables; Eq. (14) represents the integrality constraints for the variables that are not less than zero.

Moreover, when N is 0, to maintain the feasibility of Eq. (12), Eq. (12) can be replaced with Eq. (15) that can be formulated as follows:

Lower level programming model

As a scarce resource, tugboats need to be matched according to the tugboat matching principle. The port side of the lower-level decision maker aims to minimize the total number of different types of tugboats matched for each task and obtain the optimal tugboat matching scheme for each task, which will be fed back to the tugboat scheduling plan of the upper-level decision maker. The planning model M2 can be formulated as follows:

Wherein, Eq. (16) aims to minimize the total number of tugboats for the available tasks. Equation (17) ensures the number of tugboats required is allocated to serve a given task, and \({\lambda _i}\) is the dual variable of this constraint condition. For a task, according to the requirements of the port, a total of D tugboats may be necessary, among which Qtugboats with Phorsepower can work together to complete the tugboat auxiliary phase. However, sometimes, as long as the total power of the tugboat serving the task is greater than the\(Q \times P\)horsepower, the power supply demand can also be met. Therefore, Eq. (18) restricts the number of tugboats with rated horsepower required for a given task to not be lower than the minimum required for the task or to require horsepower greater than the total rated horsepower.\(\mu _{i}^{1}\)and\(\mu _{i}^{2}\)represent the dual variables of each constraint condition. Equation (19) represents the relationship between upper and lower decision variables, with\({\eta _{ki}}\)being the dual variable of the constraint condition. Equation (20) indicates that the main variables are 0–1 variables.

Bi-level model framework

Stackelberg game

Generally speaking, the Stackelberg game planning problem involving two players can be characterized using a bilevel programming model. Because in the game-theoretic process, the dominant player makes decisions first, and the upper level of the model takes the dominant player as the decision-making subject along with its decision variables. The follower makes decisions after the leader, so it corresponds to the lower-level optimization problem and its decision variable. Due to the sequential nature of decision-making, for the lower-level model, the dominant decision variable is considered as the known parameter that is passed from the upper level. In addition, in the proposed optimization framework, the upper-level model imposes constraints, reflecting the impact of follower decision-making on the leader. Based on this, it can be seen that the bilevel programming model can well reflect the mutual constraint relationship between the leader and follower in the Stackelberg game planning process. The general expression for the above Stackelberg game-theoretic problem can be described using the bilevel programming model (21).

Where\(F\left( {x,y,\lambda ,\mu } \right)\)and\(f\left( {x,y} \right)\)are the objective functions for the upper- and lower-level optimization problems, respectively.\(h\left( {x,y} \right)=0\)is the equality constraint of the lower-level optimization problem, and its corresponding dual variable is\(\lambda\).\(g\left( {x,y} \right) \geqslant 0\)is the inequality constraint of the lower-level optimization problem, whose corresponding dual variable is\(\mu\). \(H\left( {x,y,\lambda ,\mu } \right)=0\)and\(G\left( {x,y,\lambda ,\mu } \right) \geqslant 0\) are the equality and inequality constraints for the upper-level optimization problem, respectively. The optimization variable of the lower level model is the decision variable of the follower. The optimization variables of the upper model include not only the decision variable y of the leader, but also the optimization variable x of the lower model, as well as the dual variables of the lower model.

Stackelberg game-theoretic models are very different from traditional optimization models and require specific solution approaches. At present, the commonly used algorithms for solving Stackelberg game-theoretic models include the following³⁴: (1) the Karush–Kuhn–Tucker (KKT) optimality condition method, which uses the KKT optimality conditions to replace the lower-level decision problem, thereby transforming the hierarchical game-theoretic model into traditional linear programming; (2) the Penalty function method, which converts the original problem into a single-layer optimization problem through the principle of penalty function; (3) Descending direction method, which utilizes the commonly used steepest descent method and rotation scale-based methods in nonlinear programming to find local or global optimal solutions; (4) Global optimization method, which converts the original problem into a concave minimization problem and uses global optimization methods to solve it. In addition, some intelligent optimization algorithms and metaheuristics, such as particle swarm optimization, are also used to solve Stackelberg game-theoretic models. This article uses the KKT optimality condition method to solve the Stackelberg game-theoretic model (21).

If the lower level model is a linear programming problem, it can be equivalently replaced with the corresponding KKT condition, and then the original bilevel model can be transformed into a single-layer optimization problem, such as model (22).

Wherein,\(0 \leqslant \mu \bot g\left( {x,y} \right) \geqslant 0\; \Leftrightarrow \;\mu \geqslant 0,g\left( {x,y} \right) \geqslant 0,\mu \leqslant Mz,g\left( {x,y} \right) \leqslant M\left( {1 - z} \right)\),\(z \in \left\{ {0,1} \right\}\).

At this point, the original bilevel programming model can be transformed into a mixed integer linear programming problem as follows:

The optimization model proposed in this study can be solved using the aforementioned method, that is, the lower-level linear programming with the KKT optimality conditions (reserved Eq. (19)):

Fuzzy parameters

Many scholars use fuzzy theory to deal with uncertainty problems, such as Huang et al.³⁵ using a compromise fuzzy programming method to construct regional sewage system planning. However, according to the study conducted by Dehshiri et al.³⁶ and their proposed definition of the Me measure, the following equivalent formula can be derived:

If\(\lambda ^{\prime}\) = 1, then Me = Pos (where Pos is the possibility measure), indicating that the decision maker considers the possibility measure and uses the maximum probability. If\(\lambda ^{\prime}\) = 0, then Me = Nec (where Pos is the necessity measure), indicating that the decision maker considers the necessity measure and uses the minimum probability; If\(\lambda ^{\prime}\) = 0.5, then Me = Cr, the decision maker considers the feasibility measure and uses the average probability.

Based on the above analysis, Eq. (9) can be replaced by Eqs. (30) and (31) as follows:

Wherein, \(\alpha\) and \(\beta\) are the confidence levels.

Transformation of the objective functions

The upper-level programming model is a multi-objective programming model that could be challenging to solve directly. Therefore, the appropriate steps need to be taken to convert the multi-objective model into a single-objective model, which include the following.

Step 1: Linearize objective function 3.

Objective function 3 can be updated as follows:

And a new inequality should be added in the optimization model:

Step 2: Linearize objective function 2.

Objective function 2 contains nonlinear components (multiplication of two variables). Let \(Y_{{klij}}^{1}={y_{kj}}{x_{klij}}\),\(Y_{{klij}}^{2}={s_{ki}}{x_{klij}}\)and\(Y_{{klij}}^{3}={y_{ki}}{x_{klij}}\), then objective function 2 can be repressed as follows:

New constraints should be added in the optimization model:

\(Y_{{klij}}^{1},Y_{{klij}}^{2},Y_{{klij}}^{3} \geqslant 0,\quad \forall k \in K,l \in L,i \in I,j \in I^{\prime}\)

Step 3: Convert the three objectives into a single objective function.

Firstly, define each objective function as a linear membership function:

Where \({\mu _j}(v)\) represents the satisfaction of the j-th objective function.

In addition, according to the method proposed by Bilgen³⁷, a payoff table can be used (see Table 1), where\({f_t}(X)\)is the t-th objective function and\({X^{(t)}}\)is its optimal solution, to ensure the feasibility of each objective range.

The lower and upper bounds of the minimization and maximization problems can be determined from the matrix in Table 1 using Eqs. (47) and (48), respectively:

Then there are \(f_{1}^{{PIS}}=\hbox{min} \left\{ {{f_{11}},{f_{21}},{f_{31}}} \right\}\),\(f_{1}^{{NIS}}=\hbox{max} \left\{ {{f_{11}},{f_{21}},{f_{31}}} \right\}\),\(f_{2}^{{PIS}}={\hbox{max} _2}\left\{ {{f_{12}},{f_{22}},{f_{32}}} \right\}\),\(f_{2}^{{NIS}}=\hbox{min} \left\{ {{f_{12}},{f_{22}},{f_{32}}} \right\}\),\(f_{3}^{{PIS}}=\hbox{min} \left\{ {{f_{13}},{f_{23}},{f_{33}}} \right\}\) and \(f_{3}^{{NIS}}=\hbox{max} \left\{ {{f_{13}},{f_{23}},{f_{33}}} \right\}\).

Then, using the aggregation function of Torabi and Hassini³⁸, the three objective functions can be transformed into the following formula:

The following constraints should be added in the optimization model:

Where\({\lambda _0}=\mathop {\hbox{min} }\limits_{j} \left\{ {{\mu _j}(v)} \right\}\),\(\psi\)is the minimum satisfaction level for the objectives, \(\psi \in \left[ {0,1} \right]\), \({\theta _j}\)are the coefficients determined based on the decision maker’s preference, and\(\sum\nolimits_{{j=1,2,3}} {{\theta _j}} =1,{\theta _j} \geqslant 0\).

In summary, the Stackelberg game-theoretic models M1 and M2 can be transformed into a single-objective mixed integer linear programming model M3 with Eq. (49) as the objective, equations from (4) to (8), equations from (10) to (14), Eq. (19), equations from (24) to (27), equations from (30) to (31), Eq. (33), equations from (35) to (46), and equations from (50) to (51) as constraints.

The model M3 constructed in the previous section is a mixed integer linear programming model, which can be directly solved using popular commercial solvers, such as CPLEX, LINGO, and GAMS. However, as the scale of the problem increases, the efficiency of existing exact optimization solvers becomes quite low. In this case, heuristic and metaheuristic algorithms, such as genetic algorithms, nested segmentation algorithm, and simulated annealing algorithm, can be used to improve the solving efficiency and obtain good-quality solutions in a short period of time.

The seagull optimization algorithm (SOA) is a new swarm intelligence optimization algorithm proposed by Dhiman and Kumar³⁹, which mainly simulates the migration of seagulls in nature and their attack behavior (foraging behavior) during migration. This algorithm simulates the behavior of seagulls during migration and predation, constructs a mathematical model, and achieves the goal of finding good-quality solutions to the combinatorial optimization problem. After the model construction is completed, the stability and superiority of SOA are verified through a large number of test functions and comparative experiments. Finally, SOA has been introduced to solve large-scale industrial engineering problems, and its strong practicality has been verified using different industrial optimization problems. In addition, SOA has also been applied to a variety of decision problems, such as multi-reservoir systems⁴⁰, cogeneration systems⁴¹, predictive models⁴², overcurrent relays², and shield tunneling prediction⁴³.

However, like other swarm intelligence optimization algorithms, SOA also has certain drawbacks, such as premature convergence, relatively weak global optimization ability, poor population diversity, and susceptibility to local optima. The main reason for this is the lack of communication between individuals in the seagull population in the SOA algorithm. Therefore, to overcome the existing limitations of the seagull optimization algorithm, this study utilizes genetic operators, such as crossover and mutation operations widely deployed in traditional genetic algorithms^44,45, to increase the exchange of information between individuals, thus expanding the search scope of seagull groups and improving the search capability of SOA.

The priority encoding method has been successfully applied in different optimization algorithms for solving large-scale strategic logistics network design problems^46,47,48,49. This encoding method has the characteristics of simple operation and easy generation. Based on the characteristics of tugboat scheduling problems and the structure of priority encoding methods, this paper proposes a priority encoding scheme with mooring bases. In summary, this study presents a Seagull Optimization Algorithm based on Priority Encoding and Genetic Operators (SOAPG) to solve model M3. The algorithm flow is shown in Fig. 3.

The main steps of the SOAPG algorithm.

Full size image

Individual encoding and decoding operations

Individual encoding operation

Considering that this study focuses on the allocation problem between tugboats, tasks, and docking bases, after studying the solution structure, it is necessary to determine which tugboats are serving each task, and whether a given tugboat will directly go to the next task or dock at a certain base after serving. Therefore, we have designed a customized coding structure inspired by the investigation of priority encoding methods in logistics networks. Based on the symbols defined in section “Symbol description”,\(\left| K \right|\)represents the total number of tugboats and\(\left| I \right|\)represents the total number of tasks. Randomly generate a row matrix\({v_{\left| K \right| \times \left| I \right|}}\)of 1 row\(\left| K \right| \times \left| I \right|\)containing integers from 1 to\(\left| K \right| \times \left| I \right|\), known as the priority matrix, and add a number after each priority, which is selected from the set \(R=\{ 0,1,2, \ldots ,NL\}\). In addition, due to the distinction between tugboats and tasks, in individual decoding operations, it is necessary to form a row-column integer matrix\({v_{\left| K \right|,\left| I \right|}}\). The specific operation method is to arrange the numbers of\({v_{\left| K \right|,\left| I \right|}}\)in order from left to right and then from top to bottom. Randomly generate a row-column integer matrix \({v_{ki}}\) containing the values from 1 to \(\left| K \right| \times \left| I \right|\), which is the priority matrix, and add a number after each priority, which is selected from set\(R=\{ 0,1,2, \ldots ,NL\}\).

If\(\left| K \right|\)is 3,\(\left| I \right|\)is 5, and\(NL\)is 2, then the feasible priority encoding schemes\({v_{\left| K \right| \times \left| I \right|}}\)and\({v_{\left| K \right|,\left| I \right|}}\) can be represented as the ones shown in Tables 2 and 3, respectively. More information regarding the decoding operations (i.e., how the information presented in Tables 2 and 3 can be translated into a solution for the decision problem studied herein) is provided in section “Individual decoding operation” of the manuscript.

Individual decoding operation

This operation continues to use the symbolic meaning described in section “Symbol description” and proceeds as follows:

Step 1: Determine the input parameters, including\({D_i}\),\({P_k}\),\({m_i}\),\({Q_i}\), N, and priority\({v_{ki}}\),\(k \in K\),\(i \in I\). Determine the output variables, including\({x^{\prime}_{ki}}\).

Step 2: Determine the selection of 0 or 1 for the matrix\({x^{\prime}_{ki}}\)based on the inequality relationship in formulas (12), (17), and (18) and the priority relationship in the direction of the priority column\({v_{ki}}\). In addition, because formula (18) contains two formulas, one of them is selected as the matching constraint with a 50% probability during this step of operation.

For example, continuing to use the data from Table 3, add the values for\({m_i}\),\({Q_i}\),\({D_i}\),\({P_k}\)as shown in the first to third rows and second columns of Table 4. In addition, assuming N is 0.5, according to formula (12), the value range of\(\sum\limits_{{i \in I}} {{{x^{\prime}}_{ki}}}\)is between 1.333 and 6.666, indicating that the number of tugboats serving the task ranges from 2 to 6. Finally, obtain a 0–1 matrix from the fourth to sixth rows and the third to seventh columns in Table 4, which is the matrix\({x^{\prime}_{ki}}\). For the fourth to sixth rows of the third column, if it is logical to choose the numbers 15, 7, and 3 in the second column of Table 3, then the position where the maximum number 15 is located is 1, and the rest are 0. The numbers in the fourth to sixth rows of the third column can be represented as [1,0,0]^T. The same process can be applied to the other tasks.

Step 3: Calculate the start time by inputting parameters, including\({B_{kl}}\),\(c_{{li}}^{1}\),\(c_{{li}}^{2}\),\(c_{{ij}}^{3}\), formula \(\frac{{(\alpha - \lambda )w_{i}^{d}+(1 - \alpha )w_{i}^{c}}}{{1 - \lambda }}\), and priority\({v_{ki}}\)(horizontal direction):

Firstly, use Eq. (52) to calculate the time from a certain tugboat’s base to the first task:

Where i represents the task with the highest horizontal priority under the selected tugboat matching scheme, and l represents the initial docking base of tugboat k.

For subsequent tasks, Eq. (53) can be used as follows:

Where nrepresents the first scheduled arrangement of preceding task n,\(n+1\)represents the current task arrangement,jrepresents the task corresponding to the second highest priority after task i, and lrepresents the docking base corresponding to the priority of task i.

For example, based on the data in Tables 3 and 4,\(vc\)is 5,\(\alpha\)is 0.9,\(\lambda\)is 0.1, and

For tugboat 1, due to the selection of tasks 1 and 4, and based on their priorities 15 and 14, the order in which tugboat 1 executes tasks is 1, 4, and the following relationships will apply:

,

.

For tugboat 2, due to the selection of tasks 2, 3, and 5, based on their priorities 13, 11, and 10, the order in which tugboat 2 executes tasks is 2, 3, and 5

,

,

.

Features of the seagull optimization algorithm

The standard seagull optimization algorithm obtains the optimal flight position by simulating the migration and attack behavior of seagulls within a specific search area. Seagulls are a globally distributed group of seabirds that migrate between different regions in search of food as the seasons change. The predation process of seagulls mainly includes the migration stage and the attack stage. During the migration process, when all seagulls fly towards the position with the most abundant food, individuals maintain independence according to a certain flight pattern to avoid positional collisions. Meanwhile, to obtain more food, seagulls usually engage in spiral-flight attacks on the prey.

The standard seagull optimization algorithm is only suitable for optimizing in a continuous solution space. However, the solution space of the problem studied in this article is discrete, and the standard seagull optimization algorithm needs to be discretized to make it suitable for solving tugboat scheduling problems. For metaheuristic algorithms that solve continuous problems, it is necessary to redefine the position update strategy in the algorithm to adapt to discrete encoding problems. In the cuckoo bird optimization algorithm, Zheng et al.⁵⁰ used a path reconnection strategy instead of the Lévy flight to generate new paths through swap and insert methods. Teymourian et al.⁵¹ used the 2-opt method to replace the Lévy flight. Alssager and Othman⁵² used local search algorithms to search in large neighborhoods instead of the Lévy flight but did not specify which local search algorithms were used. In addition, Gezici and Livatyali⁵³ proposed a corresponding metaheuristic optimization algorithm for discrete optimization problems based on a population-based metaheuristic optimization algorithm for continuous optimization problems. The algorithm modified the way new solutions were produced, and its application was demonstrated using the CEC2019 test function and a three-dimensional packing problem (3D-BPP) dataset with 320 samples. After studying the discretization process of various continuous algorithms, this paper adopted a discrete update strategy suitable for the priority encoding of seagull optimization algorithms.

Migration operations

Seagulls may not directly search for the neighborhood of the optimal solution during their migration process, but they tend to explore good-quality solutions around them. In order to explore good-quality solutions, facilitate the global evolution direction of the population, and ensure the flight independence of seagulls, this study deployed the strategy of inserting and reversing the best obtained path scheme with a 50% probability of developing the neighborhood of the best solution identified, so that seagulls can migrate towards the neighborhood of the best position. Taking the data in Table 2 as an example, the insert operation (see Table 5) randomly selects two nodes i and j in the best path scheme identified, inserts node j next to the position of node i, and forms a new path. The reverse operation (see Table 6) aims to randomly select two nodes i and j in the best path scheme identified, invert all nodes between node i and node j, and form a new path.

Attack operations

When seagulls attack their prey, the method of using a spiral flight is relatively common. To expand the local search range and improve the ability to discover new paths, based on the final migration ___location of seagulls, swap and 3-opt operations are performed on the current path scheme to replicate seagull attack methods. The swap operation (see Table 7) is conducted by randomly selecting two nodes i and jin the selected path scheme and exchanging the positions of node i and node jto form a new path. The 3-opt operation (see Table 8) is performed by randomly selecting three nodes i, j, and kfrom the selected path scheme and exchanging the positions of the three nodes to form a new path.

Fitness function and selection operations

Due to the need for multiple tugboats to be used simultaneously in the task, the encoding method mentioned above may cause start time conflicts after decoding. Generally speaking, there are two types of conflicts:

1. (1)

   Conflict situation 1: For task i, check whether\(s_{{ki}}^{n}\)is within the range of the task time window\(\left[ {{A_i},{A_i}+{\tau _i}} \right]\). Moreover, for task\(i \in I^{\prime}\), check whether\(s_{{ki}}^{n}+\frac{{(\alpha - \lambda )w_{i}^{d}+(1 - \alpha )w_{i}^{c}}}{{1 - \lambda }}\)is within the range of task time window\(\left[ {{A_i},{A_i}+{\tau _i}} \right]\). If not, add the corresponding penalty function to the objective function Eq. (49) as follows:

   $$\hbox{max} \;\psi {\lambda _0}+(1 - \psi )\left[ {{\theta _1}{\mu _1}(v)+{\theta _2}{\mu _2}(v)+{\theta _3}{\mu _3}(v)} \right] - {F_1} \times \sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {\left( {s_{{ki}}^{n} - {A_i} - {\tau _i}} \right)} }$$

   (54)

   Where\({F_1}\)is the penalty value.

2. (2)

   Conflict situation 2: When each tugboat is performing the same task, it must be carried out simultaneously at a certain period. It is not allowed for one tugboat to complete the task and start the next task before another tugboat starts the task that was scheduled for simultaneous execution.

   The minimum and maximum start times for a task are\(s_{i}^{{\hbox{min} }}=\hbox{min} \left\{ {s_{{ki}}^{n}} \right\}\)and\(s_{i}^{{\hbox{max} }}=\hbox{max} \left\{ {s_{{ki}}^{n}} \right\}\), respectively. The tugboat k that performs task i, for all\(j \in I\backslash \left\{ i \right\}\), if both\(s_{{kj}}^{n}<s_{i}^{{\hbox{max} }}\)and\(s_{{kj}}^{n}+{y_{kj}}>s_{i}^{{\hbox{min} }}\)are met, then task j, intersects with task i. Firstly, set an empty matrix\({\Delta _{kij}}\)for tugboat k that performs task i. For each j, if\(s_{{kj}}^{n}<s_{i}^{{\hbox{max} }}\)and\(s_{{kj}}^{n}+{y_{kj}}>s_{i}^{{\hbox{min} }}\)are not met, then\({\Delta _{kij}}\) = 0. For tugboat kthat does not perform task i, then\({\Delta _{kij}}\) = 0 for each iand j. Otherwise, proceed to the next step:

   If\(s_{{kj}}^{n}<s_{i}^{{\hbox{min} }}\)and\(s_{{kj}}^{n}+{y_{kj}} \leqslant s_{i}^{{\hbox{max} }}\), then\({\Delta _{kij}}=s_{{kj}}^{n}+{y_{kj}} - s_{i}^{{\hbox{min} }}\);

   If\(s_{{kj}}^{n} \geqslant s_{i}^{{\hbox{min} }}\)and\(s_{{kj}}^{n}+{y_{kj}}>s_{i}^{{\hbox{max} }}\), then\({\Delta _{kij}}=s_{i}^{{\hbox{max} }} - s_{{kj}}^{n}\);

   If\(s_{{kj}}^{n}<s_{i}^{{\hbox{min} }}\)and\(s_{{kj}}^{n}+{y_{kj}}>s_{i}^{{\hbox{max} }}\), then\({\Delta _{kij}}=s_{i}^{{\hbox{max} }} - s_{i}^{{\hbox{min} }}\).

   The corresponding penalty function needs to be added to Eq. (54) as follows:

   $$\begin{gathered} \hbox{max} \;\psi {\lambda _0}+(1 - \psi )\left[ {{\theta _1}{\mu _1}(v)+{\theta _2}{\mu _2}(v)+{\theta _3}{\mu _3}(v)} \right] - {F_1} \times \sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {\left( {s_{{ki}}^{n} - {A_i} - {\tau _i}} \right)} } \hfill \\ \quad \quad - {F_2} \times \sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {\sum\limits_{{j \in I}} {{\Delta _{kij}}} } } \hfill \\ \end{gathered}$$

   (55)

   Where\({F_2}\)is the penalty value.

   Therefore, the fitness value to be used for quantifying the quality of solutions in the seagull optimization algorithm can be expressed as follows:

   $$\begin{gathered} objective=\psi {\lambda _0}+(1 - \psi )\left[ {{\theta _1}{\mu _1}(v)+{\theta _2}{\mu _2}(v)+{\theta _3}{\mu _3}(v)} \right] - {F_1} \times \sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {\left( {s_{{ki}}^{n} - {A_i} - {\tau _i}} \right)} } \hfill \\ \quad \quad \quad \quad - {F_2} \times \sum\limits_{{k \in K}} {\sum\limits_{{i \in I}} {\sum\limits_{{j \in I}} {{\Delta _{kij}}} } } \hfill \\ \end{gathered}$$

   (56)

   Due to the existence of penalty values, Eq. (56) may produce negative values. Therefore, it is necessary to use Eq. (57) to turn all negative values into positive values. Let the current population set be\(t \in T\), and\(objectiv{e_t}\)be the objective value of the individual t. Then, the fitness function can be represented as follows:

   $$fitnes{s_t}=\mathop {\hbox{max} }\limits_{{t \in T}} \left\{ {\left| {objcetiv{e_t}} \right|} \right\}+objcetiv{e_t}$$

   (57)

   The roulette wheel selection method was adopted in this study as a selection mechanism. Assuming the population size is\(pop\_szie\), due to the crossover operation applied to the population in this algorithm, the population size becomes twice the original population size, i.e., \(2 \times pop\_szie\). During the selection operation, the roulette wheel needs to be rotated \(pop\_szie\) times. The individuals identified by the roulette wheel selection mechanism will be designated as seagulls that will be moved to the next iteration and continue the search process.

Crossover and mutation operations

By taking two different nodes i and jfrom two different individuals (parents), the middle part is exchanged by the developed crossover operator to replace the mapping relationship of the exchange part for scenarios where there are no exchanges but duplicate numbers with the exchange part. As shown in Table 9, after exchanging the parts between two random nodes, the changes are made through the mapping relationship of 1\(\leftrightarrow\)4\(\leftrightarrow\)13, 9\(\leftrightarrow\)12\(\leftrightarrow\)3, and 2\(\leftrightarrow\)10, as shown by the bold and underlined numbers in Table 9. After the crossover operation, the mutation operation is performed as well. The mutation operation is mainly aimed at changing the docking base, as shown in Table 10. The developed mutation operator randomly selects node i and then changes the docking base (as shown in bold underlined numbers in Table 10).

Input data generation

To evaluate the performance of the proposed Seagull Optimization Algorithm based on Priority Encoding and Genetic Operators (SOAPG), CPLEX 12.8, traditional genetic algorithm (GA), discrete seagull optimization algorithm (SOA), and simulated annealing algorithm (SA) were considered during the computational experiments as well. For fair comparison, the encoding of GA, SOA, and SA was performed following the encoding method proposed in this research. Furthermore, the crossover and mutation operations used within GA were similar to the ones described in section “Fitness function and selection operations”. The migration and attack operations used within SOA were similar to the ones described in section “Features of the seagull optimization algorithm”. The SA operations for generating new solutions were similar to the ones described in section “Migration operations”. The codes for SOAPG, GA, SOA, and SA were compiled in MATLAB 7.0. The running environment for software MATLAB 7.0 and CPLEX was composed of an Intel (R) Core (TM) i7 2.70 GHz processor, 4GB memory, and a Windows 10 (64bit) laptop operating system. This study randomly generated 45 problem instances of different scales based on Tables 11, 12, 13, 14, 15 and 16. Among them, Tables 11 and 12, and 14 were obtained based on data from the literature⁵. Table 13 was provided by internal resources of the Guangzhou Port, and Table 15 was obtained based on the information presented in the “Guangdong Province Port Ship Berthing and Departure and Pilot or Transfer Tug Ship Allocation Standards”.

Generally speaking, the speed limit for tugboats in Guangzhou Port is 10 knots, while the maximum speed in other regions is 13 knots. From the nine average speeds in Table 13, it can be inferred that the average speed of tugboats in the Guangzhou Port on each sea route is 10.62 km/h, which meets the requirements. According to Table 13, the distance from the task to the base can be set between 7 km and 28 km, and the distance between tasks can be set between 2 km and 28 km.

In addition, in practical operations, the total power method is also used, which means that the total power provided by the available number of tugs with a certain power multiplied by their respective horsepower should be greater than the required power. As shown in Table 14, when 270 ≤ L < 390, the total power of tugs serving the vessel needs to be greater than 2\(\times\)4000HP = 8000HP.

Overall, the maximum operating time for one tugboat was set to 15–45 min, for two tugboats it was set to 25–60 min, and for three tugboats it was assumed to comprise 30–75 min. When working with the available data for planning purposes, the port side generally knows in advance the vessel type and related operational information of the incoming vessels, so they can determine in advance the required configuration of tugboats for these vessels.

Evaluation of the candidate solution approaches

In order to facilitate the verification of the performance of the proposed SOAPG, objective 1 was adopted for all the four considered solution approaches, namely:

In addition, the fitness formula follows Eq. (57), and the objective value function of CPLEX follows Eq. (1). The ambiguity coefficients were assumed to be \(\alpha\)=\(\beta\)=\(\lambda ^{\prime}\)=0.5, \({F_1}\)=10⁵,\({F_2}\)=10⁶. Furthermore, a set of parameter tuning experiments were conducted to set the appropriate values of the considered metaheuristic algorithms. In particular, based on the parameter tuning analysis, the crossover and mutation rates of SOAPG and GA were set to 60%~70% and 10%~15%, respectively. Each algorithm was run 20 times for each problem instance, and the results are reported in Tables 17 and 18 (bold numbers represent the optimal solution, and those marked with “*” represent feasible solutions). The sentence “Did not converge” indicates that CPLEX identified a feasible solution but could not find the optimal solution within the time limit imposed (which is 9,000 s).

As shown in Table 17, for problem instances 1–16, the objective values obtained by CPLEX and SOAPG are the same, but the solving time used by SOAPG is fairly longer than that used by CPLEX, indicating that for small-scale examples, CPLEX has high efficiency. For problem instances 17–21, although the objective values obtained by CPLEX and SOAPG are also the same, the solving time used by CPLEX is gradually much longer than that of SOAPG, indicating that for small to medium-scale examples, the efficiency of SOAPG begins to improve. For problem instances 22–32, although the objective values obtained by CPLEX are still optimal, CPLEX required too much time to discover these optimal solutions. Although SOAPG could not find the optimal objective values for problem instances 22–32, the obtained objective values are not much larger than the optimal ones and can be viewed as acceptable. Moreover, the solving time required by SOAPG is much shorter than that required by CPLEX. This indicates that for medium to large-scale examples, the accuracy of SOAPG may slightly decrease, but its efficiency in terms of computational time improves when comparing to exact optimization methods. For problem instances 33–45, the CPLEX solution time exceeded 9,000 s, and CPLEX was able to identify a feasible solution but could not converge within the CPU time limit imposed. However, SOAPG obtained superior objective values with the solution time much less than 9,000 s. For subsequent large-scale problem instances, CPLEX could not provide a feasible solution within 10,000 s, while SOAPG could still provide good-quality solutions within reasonable computational time. Hence, according to the findings from the conducted experiments, the SOAPG algorithm has a high level of efficiency, which can be confirmed by the fact that the objective values of SOAPG are the optimal ones for small to medium-scale and small-scale problem instances. Furthermore, SOAPG can provide good solutions within an acceptable time range when solving large-scale problem instances.

Based on the results reported in Tables 17 and 18, the largest problem instances that could be solved by GA, SOA, and SA to optimality are problem instances 9, 10, and 9, respectively. SOAPG demonstrated a much more competitive performance and could obtain optimal solutions for 21 problem instances. Moreover, for all the considered problem instances, the objective values obtained by SOAPG are lower than those of GA, SOA, and SA, indicating the accuracy of SOAPG in solving the decision problem studied herein. For all the generated problem instances, the solving time of SOAPG was found to be slightly higher than that of GA, SOA, and SA. Such a pattern can be explained by the fact that during the SOAPG operations, in addition to global search, individual information exchange operations are also added. Nevertheless, the solving time of SOAPG is still within an acceptable time range. Therefore, according to the findings from the conducted experiments, the SOAPG algorithm is more prone to jumping out of local optimal solutions than traditional GA, SOA, and SA algorithms and demonstrates a high level of efficiency in terms of solution quality and computational time.

Convergence patterns of the four algorithms using problem instance 21 as an example.

Full size image

To analyze the convergence speed and search capabilities of SOAPG, GA, SOA, and SA, the convergence patterns were evaluated for problem instance 21 for the four algorithms and are presented in Fig. 4. Similar tendencies in algorithmic performance were noticed for other problem instances as well. After 400 iterations, the objective values obtained are arranged in descending order as GA, SA, SOA, and SOAPG, respectively. From Fig. 4, it can be seen that under the same number of iterations, population size, and other parameters, SOA is more likely to jump out of the local optimal solutions during the search process when comparing to GA and SA, indicating the superiority of the SOA algorithm design. After adding crossover and mutation genetic operators, SOAPG shows advantages over the standard SOA in global and local search capabilities. Such a tendency can be explained by the fact that the new genetic operators can enhance communication between individuals in the population, making it easier to jump out of the local optimal solutions during the search process.

Sensitivity analysis and managerial insights

The Guangzhou Port currently has four major port areas: Inner Port Area, Huangpu Port Area, Xinsha Port Area, and Nansha Port Area. The Inner Port of Guangzhou is composed of three port areas: Panyu, Wuhe, and Xintang. As such, six berthing bases are set up, including Huangpu Port Area, Xinsha Port Area, Nansha Port Area, Panyu Port Area, Wuhe Port Area, and Xintang Port Area. A total of 15 tugboats were considered during the sensitivity analysis, and the information for 25 tasks was regenerated using Tables 12, 13 and 14 as shown in Table 19. In addition, a maximum fuzzy delay time of (10,20,30,40) was set for all tasks.

In addition, the power of tugboats 1–15 (in HP) was set to 1600, 3000, 4000, 5000, 6000, 6800, 1600, 3000, 4000, 5000, 6000, 6900, 5000, 6000, and 6900, respectively. Before the start of the day’s operations, tugboats 2, 8, and 11 were assumed to be docked at the Huangpu Port Area berthing base; tugboats 5 and 14 were assumed to be docked at the Xinsha Port Area berthing base; tugboats 4, 9, and 13 were assumed to be docked at the Nansha Port Area berthing base; tugboats 6 and 15 were assumed to be docked at the Panyu Port Area berthing base; tugboats 1, 7, and 12 were assumed to be docked at the Wuhe Port Area berthing base; and tugboats 3 and 10 were assumed to be docked at the Xintang Port Area berthing base.

Although the optimal solution can be obtained for the 32nd problem instance, the efficiency of CPLEX significantly decreases when considering objective function 2 only, and the optimal solution can only be obtained for problem instances 1–6. Therefore, to ensure the rationality of the solution, this article applies both CPLEX and SOAPG algorithms to solve this problem instance. First, CPLEX was used to solve the optimization model with objective functions 1 and 3, then the SOAPG algorithm was used to solve the optimization model with objective function 2, and finally the SOAPG algorithm was used to solve the optimization model with all the objective functions considered and associated weights. When N=\(\alpha\)=\(\beta\)=\(\lambda ^{\prime}\)=0.5, solving the model with three objectives yielded the results presented in Table 20. Thus, the maximum value of objective 1 is 561,613.96, and the minimum value is 92,264.30; the maximum value of objective 2 is 5,323.22, and the minimum value is 0; and the maximum value of objective 3 is 1,125 and the minimum value is 1,095.

Next, the SOAPG algorithm was executed for various scenarios that were developed by changing the minimum satisfaction\(\psi\)level of the objectives, the decision maker’s preference estimates\({\theta _1}\),\({\theta _2}\)and\({\theta _3}\), the fairness coefficient N, and the values of fuzzy parameters\(\alpha\),\(\beta\)and\(\lambda ^{\prime}\). The recorded objective function values are presented in Table 21 for all the considered scenarios. Furthermore, Figs. 5, 6 and 7 were developed based on the data obtained for scenarios 12–16, 6, and 17–21; Figs. 8, 9, 10 and 11 were developed based on the data obtained for scenarios 26, 6, and 27; and Figs. 12, 13 and 14 were developed based on the data obtained for scenarios 28, 6, and 29.

Behavior of the weighted total objective under the minimum satisfaction level change.

Full size image

Behavior of objective 1 under the minimum satisfaction level change.

Full size image

Behavior of objective 2 under the minimum satisfaction level change.

Full size image

Behavior of the weighted total objective under the change of measure probability.

Full size image

Behavior of objective 1 under the change of measure probability.

Full size image

Behavior of objective 2 under the change of measure probability.

Full size image

Behavior of objective 3 under the change of measure probability.

Full size image

Behavior of the weighted total objective under the confidence level changes.

Full size image

Behavior of objective 1 under the confidence level changes.

Full size image

Behavior of objective 2 under the confidence level changes.

Full size image

By observing the results obtained for scenarios 1–11 in Table 21, while keeping the coefficients other than the fairness coefficient constant, the weighted total objective changed from 0.8395 at the fairness coefficient of 0 to 0.8408 at the fairness coefficient of 0.1 and remained unchanged at 0.8408 thereafter. Objective 1 changed from 236,185.90 at a fairness coefficient of 0 to 234,939.70 at 0.1 and remained unchanged thereafter. Objective 2 and objective 3 remained unchanged at 3,698.22 and 1,095, respectively. This indicates that the change in the fairness coefficient has an insignificant impact on the fair allocation of tugboat task volume in the tugboat scheduling problem studied herein. Such a pattern can be justified by the fact that, in this example, there are a certain number of tugboats with the same horsepower that can meet the needs of 25 tasks. When the fairness coefficient changes, the tasks assigned to the original tugboat can be transferred to another tugboat with the same horsepower for operation, so the objective function value can be kept unchanged when the model solution changes. To illustrate the impact of coefficients other than the fairness coefficient, as shown in Table 21 for the data obtained for scenarios 22–25, when changing coefficients other than the fairness coefficient, the objective function trends appeared to be the same as those of recorded for scenarios 1–11.

By observing the data obtained for scenarios 12–16, 6, and 17–21 in Table 21; Figs. 5, 6 and 7, while keeping the coefficients other than the minimum satisfaction of the objective unchanged, the weighted total objective exhibited a linear decreasing trend during the change of the minimum satisfaction level of the objective from 0 to 1; objective 1 and objective 2 exhibited stable changes after increments from 0 to 0.3, respectively; and objective 3 remained stable during the change in the minimum satisfaction level from 0 to 0.9 and suddenly increased to 1,104.16 when it reached 1. This indicates that the change in the minimum satisfaction of the objectives had a significant impact on the tugboat scheduling problem studied herein. Decision makers of the tugboat company should choose a satisfactory tugboat plan suitable for port operations based on the actual conditions in the port.

By observing the data obtained for scenarios 26, 6, and 27 in Table 21; Figs. 8, 9, 10 and 11, while keeping the coefficients other than the measurement probability constant, the weighted total objective and objective 2 show an increase followed by a decrease after changing the measurement probability from 0 to 1; objective 1 shows a decreasing change; and objective 3 shows a decrease followed by an increase. This indicates that the impact of necessity measurement, possibility measurement, and feasibility measurement on the tugboat scheduling problem studied herein is obvious. Decision makers involved in the tugboat operations should fully consider the tugboat operational conditions after the vessel arrives at the port and choose a suitable measurement to develop a tugboat plan. According to the results in the provided example, the feasibility measure is the most suitable because it can maintain the maximum weighted total target satisfaction of the entire tugboat operations, lower fuel consumption, maximum total buffer time for dynamic tasks, and minimum latest end time.

By observing the data obtained for scenarios 28, 6, and 29 in Table 21; Figs. 12, 13 and 14, while keeping all the coefficients except the confidence level constant, the weighted total objective and objective 2 exhibited a quasi-linear decreasing change after the confidence level change from 0 to 1; objective 1 showed a stable and unchanged state after increasing; and objective 3 showed a steady and increasing state. This indicates that the change in the confidence level has a significant impact on the tugboat scheduling problem studied herein. Decision makers of the tugboat company should choose an appropriate confidence level based on the historical operational data of the tugboats used by the incoming vessels, as well as the hydrological and weather conditions at the port at that time, to ensure that the tugboat plan can be executed with low risk.

By comparing the solutions for objectives 1, 2, and 3, it was found that for the solutions for objectives 1 and 3, when tugboats serve the task, they can be allocated and scheduled according to the number of tugboat requirements\({D_i}\)for task i. However, for the solutions for objective 2, to maximize the total buffer time, the optimal solution for objective 2 assigns all the tugboats that meet the horsepower requirements to serve dynamic task 3. Although there is some irrationality in such a solution (referring to the allocation of extra tugboats to serve dynamic tasks), it indicates that if one wants to respond to the arrival of dynamic tasks and immediately dispatch tugboats for service, the decision makers of the tugboat company should keep all the tugboats that meet the horsepower requirements idle at the expected time of arrival of dynamic tasks to maintain the requirements of objective 2. From this perspective, the solutions satisfying objective 2 can be also viewed as reasonable.

The essence of the tugboat scheduling problem discussed in this study lies in the game-theoretic interaction between the tugboat scheduling of the tugboat company and the tugboat allocation of the port dispatcher. The Stackelberg game-theoretic method and multi-objective nature of the proposed model precisely solve the irrationality of objective 2, that is, adjusting the weights of the three objectives and the minimum satisfaction of the objectives to achieve a certain balance between the three goals of the tugboat company and the goal of port dispatcher in minimizing of the number of tugboats equipped for each task. As shown in the results for scenario 6 in Table 21, the weighted total objective satisfaction value is 0.8408, objective 1 is 234,939.70, objective 2 is 3,698.22, and objective 3 is 1,095.00. The analysis of the solution obtained for scenario 6 shows that 10 tugboat services were assigned for dynamic task 3, other tasks maintained the minimum tugboat demand, and all the coefficients had their intermediate values.

Due to the setting of the third task as a dynamic task in this section, it cannot represent generality. At the same time, when solving the 29 problem instances in this section, it took more than one hour to solve the model using CPLEX. To test the impact of dynamic tasks on the model and reduce the solving time, without changing the other input data (i.e., the coefficient values were the same as the values used in scenario 6 in Table 21), the following changes were made. Three objective functions were used as the objective functions of the M3 model, and then the 3rd, 8th, 13th, 18th, and 23rd tasks were set as dynamic tasks. The maximum waiting times for these five tasks were set to 300 min, 300 min, 180 min, 300 min, and 240 min, respectively. Finally, Eq. (58) was added as a supplementary constraint before executing CPLEX.

The results show that the optimal solution for objective 1 could be obtained in 20 min (less time than the 21 problem instances described in Table 21), while 10 min were required for objectives 2 and 3 (see Table 22). The faster solution time can be explained by a tighter upper bound for objective function 2, reducing the search range and thus reducing the overall solution time required.

Because CPLEX takes more than two hours to solve the optimization model with the total objective combining objectives 1 through 3, the SOAPG algorithm can be viewed as the most practical solution approach. However, the introduction of Eq. (59) allowed achieving competitive computational time even for the cases with three objective functions considered at the same time.

In recent years, an increase in China’s foreign trade cargo transportation has led to an increasing number of vessels docked at coastal ports, and ports are facing enormous operational pressure. As the main elements assisting vessels in entering and exiting ports, tugboats need to be scientifically managed. The current research proposed an optimization model for allocation and scheduling of tugboats. The developed model considered the impact of decisions of the tugboat operator and the port scheduling party on tugboat scheduling under the scenario of dynamic arrival of tasks and fuzzy tugboat operation time. The tugboat company aimed to minimize the total daily fuel consumption of tugboat operations, the total buffer time of dynamic tasks, and the total completion time as the objective functions, while the port dispatcher minimized the total number of tugboats allocated for the given tasks. The interactions between the tugboat operator and port dispatcher were modeled in the Stackelberg game-theoretic settings.

Based on the characteristics of the decision problem studied herein, a new priority encoding-based solution method was proposed and applied within the framework of the seagull optimization algorithm, which also incorporated genetic operators (i.e., the seagull optimization algorithm based on priority encoding and genetic operators – SOAPG). Finally, a comprehensive comparative analysis with the exact optimization method (CPLEX), genetic algorithm (GA), standard seagull optimization algorithm (SOA), and simulated annealing (SA) algorithm was conducted for a wide range of different problem instances to assess the effectiveness of the proposed algorithm. Numerous experimental results have shown that even for large-scale numerical examples, the proposed algorithm could provide optimal and high-quality solutions within acceptable computational time.

The results indicate that the proposed priority-encoding method was convenient and suitable for expressing the features of the problem studied in this research. In addition, the difference in CPU time required by SOAPG to solve the generated problem instances was not much higher than the time required by GA, SOA, and SA. In addition, compared to GA, SOA, and SA, the SOAPG algorithm achieved a smaller objective value of total fuel consumption in different situations. This indicates that the SOAPG algorithm was able to achieve better solution quality and robustness. Moreover, SOAPG demonstrated competitive performance in terms of solution quality and computational time when comparing to the exact optimization method (CPLEX). The applicability of the presented optimization model was also demonstrated through the sensitivity analysis, which was conducted by changing the fairness coefficient, uncertainty parameter correlation coefficients, and objective function correlation coefficients. The current research is of great significance for optimizing tugboat allocation and scheduling problems with uncertain situations and game-theoretic relationships. The results of this study are expected to provide scientific support for relevant port enterprises and stakeholders.

The current research can be further expanded in several aspects. In the future, we can fully consider how the tugboat companies and port dispatchers can collaborate in planning tugboat plans and develop more effective schedule alternatives. In addition, the encoding method proposed in this study can be considered universal and needs to be further deepened for its application in different maritime scenarios. Furthermore, better strategies need to be proposed to address the uncertainty of dynamic task arrivals, such as adopting a multi-scenario strategy, discretizing the period of dynamic arrival, studying the results brought by each arrival situation, and using certain robust methods to solve the resulting optimization model. Last but not least, the developed SOAPG algorithm can be compared to more advanced optimization algorithms, including customized heuristics, adaptive algorithms, island algorithms, diffused algorithms, and hyperheuristics.