Mathematical model and simulation of cooperative manoeuvres among a ship and tugboats

Summary Assisting harbour tugboats are essential for safe navigation in ports and for berthing/unberthing operations. Because the port infrastructure needs to be updated in the coming age of autonomous ships, autonomous tugboats or remotely controlled tugboats are expected to be part of the navigation assistance provided by future ports. Although there have been some studies on the design of cooperative control of tugboats, few studies have focused on the mathematical model of their cooperative manoeuvres. Particularly, in the case of pushing assistance, tugboats have often been treated simply as a kind of side thrusters attached to the assisted ship. Thus, we present a new framework of the mathematical model for cooperative manoeuvres that considers the coupled motions among tugboats and a ship as precisely as possible. Solving the dynamics of all tugboats as well as the assisted ship can render the model more advanced and realistic, and is the most significant contribution of this study. The simulation tool based on the proposed model can be used as a plant model in designing and verifying the tugboat’s manoeuvring control system in the future. In this study, as examples, considering a tentative control method, some unique scenarios were simulated to demonstrate the cooperative manoeuvres.


Introduction
A project in the 1980s aimed to develop a highly reliable and intelligent ship [1], and some innovative elemental technologies suggested the realisation of automated navigation ships in the future. Presently, communication and sensing technologies, which were considered challenging then, have made significant progress, and intensive research and development are being promoted worldwide. Regarding manoeuvring control, emphasis has been placed on the development of technologies that allow a ship to complete manoeuvring on its own. For example, automatic collision avoidance systems [2,3,4] and automatic berthing systems [5,6] are of interest and have received significant attention. However, the control force of a ship, such as the rudder or side-thruster force, is not necessarily sufficient for its large inertial force, which makes it challenging to establish a control system for various scenarios. In order to compensate for the lack of a control force, assistance from harbour tugboats is essential for safe ship handling, especially for large sea-going ships that usually have poor manoeuvrability. Indeed, the methodology to determine the suitable tugboat based on the type of propulsion system [7] or optimal structure of the tugboat fleet [8] are still being studied these days. This fact indicates that the importance of tugboats will not change with the coming age of autonomous ships.
Devaraju et al. [9] insisted that port infrastructure also needed to be upgraded to satisfy the requirements of developed autonomous ships and referred to autonomous tugboats as one of navigation assistance provided by future ports. The development of autonomous tugboats or remotely controlled tugboats is expected to reduce the risk involved when the tugboat operates near a ship that is underway [10] and the risk of clashing into the berth due to human error during berthing operations [11]. Autonomous tugboat assistance is also expected to save time and reduce costs for shipping companies and port operators [12]. Based on these advantages, several companies have been developing autonomous tugboats and related technologies, and their progress has been announced in web articles [13,14,15]. Sea trials of the developed tugboats have also been reported [16,17], but their purpose was to test autonomous and remote navigation systems. A few studies have discussed cooperative manoeuvres of multiple tugboats engaged in assistance operations so far. According to Devaraju et al. [9], navigation assistance by autonomous tugboats is categorised as being at the mid-term (i.e., 10-25 years) technology readiness level. Therefore, this area is still a research hotspot and motivates academic studies.
Du et al. [18] classified several recent studies based on autonomous ship berthing ways (i.e., self-berthing and assisted berthing), control methods, and disturbance considerations. Regarding the studies on assisted berthing, many of them on pushing assistance assumed that the tugboats were fixed around the ship and that the locations and orientations of each tugboat in relation to the ship were given. For example, simulation results for position tracking [19,20,21], berthing [22], and conveying a large surface vessel [23] have been presented. Feemster and Esposito [24] demonstrated their algorithm in an open-water experiment. Several studies have been conducted on towing assistance. Hajieghrary et al. [25,26] solved the cooperative control of two leading autonomous surface vehicles (ASV) towing a buoyant load. Du et al. [18] solved a similar problem for multiple autonomous tugboats for ship towing with and without environmental disturbances [27]. Lee et al. [28] also designed the control system to convey a barge by four tugboats where two tugboats towed the barge, while the other two pushed it. Meanwhile, Chen et al. [29] solved a formation-tracking problem in which ASVs connected to an object by towlines transported that object while maintaining its formation. Thus, although there have been many considerations in the design of cooperative control of ASVs, few studies have focused on the mathematical model of the cooperative manoeuvres of multiple tugboats. Particularly, in the case of pushing assistance, tugboats have often been treated simply as a kind of side thrusters attached to the assisted ship, and their dynamics have not been solved.
Thus, we present a new framework of the mathematical model for cooperative manoeuvres that considers the coupled motions among tugboats and a ship as precisely as possible, and that is able to handle their time-changing modes, i.e., pushing mode, towing mode, and mutually-unconstrained free-moving mode, switching during the assistance. Solving the dynamics of all tugboats as well as the assisted ship taking their geometrical interaction into account can render the model more advanced and realistic. It is the most significant contribution of this study. Because the simulation platform based on the proposed model can evaluate manoeuvring motions of every player in the same field, it can be used as a plant model to design and verify the tugboat's control system, assuming autonomous or remotely controlled tugboats in the future. In this study, as examples, considering a tentative control method, some scenarios were simulated to demonstrate cooperative assisting manoeuvres.

Coupled manoeuvring model of a ship and tugboats: Pushing assistance
A situation where number of tugboats push a large ship to assist its manoeuvring is considered. Typically, a tugboat pushes a predetermined point on the side hull, at which the Mathematical modelling and simulation of cooperative manoeuvres Masaaki Sano among a ship and tugboats structure of the ship is strengthened, to prevent the hull from denting. Therefore, coupled equations of motion of a ship and multiple tugboats are formulated assuming that each tugboat is coupled with the ship via a one-pin connection.
The right-hand coordinate system used in this study is illustrated in Fig. 1. The spacefixed coordinate system with the plane referring to the still water surface is defined as -, where the heading angle of the ship 0 and that of the i th tugboat are defined. Hereafter, the variable with subscript 0 belongs to the ship, and the variable with belongs to the i th tugboat. Notably, these subscripts are omitted when they are written in a unified notation. A hull-fixed coordinate system, that is, -is also considered, in which the origin is defined at midship, and the x-and y-axes are on the bow and starboard sides of the hull, respectively.
is the longitudinal position of the centre of gravity of each hull. The velocity components in the x and y directions are expressed as and , respectively, and is the yaw angular velocity around the z-axis. The velocity and hull drift angle are denoted as and , respectively.
Regarding the geometrical variables, and represent the distance and plane angle from the midship of the ship to the i th pushed point on 0 -0 0 0 , respectively. Meanwhile, is the distance from the midship of the i th tugboat to the pushing point on the tugboat that is assumed at the centreline, that is, is half the length of the tugboat. In addition, and are the pushing forces acting on the ship owing to the i th tugboat or the reaction force of that tugboat. They are defined in the longitudinal and lateral directions on -, respectively. is the thrust of the i th tugboat and its direction is denoted by .

Motion equations
The framework of the manoeuvring mathematical model is referred to as the MMG model [30], which is a standard model and still widely used for manoeuvring simulations. (e.g., [31]). In this study, it is developed to be able to deal with cooperative manoeuvres.
Where 0 is the mass and 0 and 0 are the added masses in the surge and sway of the ship, respectively. 0 is the moment of inertia, and 0 is the added moment of inertia. On the righthand side (RHS), 0 , 0 , and 0 represent the hydrodynamic surge force, sway force, and yaw moment due to the manoeuvring motions, respectively. The second term on the RHS refers to the total forces and moment acting on the ship as it is pushed by number of tugboats. Because the ship is assumed to be completely assisted by tugboats, the rudder and propeller forces are not included in these equations; however, they can be easily considered by adding extra terms in the RHS.

i th tugboat engaged in pushing assistance
Similarly, the 3-DOF motion of each tugboat was considered. The equations of motion for the i th tugboat are as follows.
Where the first term on the RHS of each equation is the hydrodynamic force owing to the manoeuvring motions. The reaction forces received by pushing the ship are considered in the second term of the RHS. The third term refers to the thrust, assuming an azimuth-type thruster positioned at ( , 0) on -. Although a tugboat may be equipped with multiple thrusters, a single thrust vector term is considered in these formulations. This signifies that they are assumed to be simply controlled in parallel and to generate the same thrust; thus, they are assumed to function as one thruster.

Hydrodynamic forces acting on the hull
The ship and tugboats may experience large manoeuvring motions during their assisting operation. In general, when an inflow angle to the hull becomes large, the formulation of its hydrodynamic forces becomes complicated [32]. Here, , , and are expressed in reference to the so-called crossflow drag model developed by Yoshimura [33], which can cover a wide range of manoeuvring motions with a small number of coefficients. Thus, it is practically useful.
Mathematical modelling and simulation of cooperative manoeuvres Masaaki Sano among a ship and tugboats Where is the water density and and are the ship length and draft, respectively. The constants , , , , and , are called the hydrodynamic force derivatives of manoeuvring. They represent the magnitudes of the manoeuvring force and moment. Regarding , the resistance is considered as the first term, and the constant 0 switches according to the forward or backward direction. and have a characteristic nonlinear force term in the third term. They are calculated by integrating the crossflow drag along the longitudinal direction, where 0 is the drag coefficient and and are the modification coefficients of the crossflow due to yawing.

Geometrical position relationship in pushing assistance
In this study, the ship and tugboats were coupled using a pin connection, and their positions were dependent on each other. The coordinates of the i th tugboat ( , ) in the space-fixed coordinate system have the following relationship with the coordinates of the ship ( 0 , 0 ).
Taking the second order derivative with respect to time, the continuity of the plane accelerations of the ship and i th tugboat are guaranteed by the following equations. Where and are the surge and sway accelerations, respectively, in the space-fixed coordinate system. As ̇ and ̇ in Eqs. (12) and (13) (17) for the n number of tugboats, and Eqs. (12) and (13) for their geometrical position relationship are considered as linear simultaneous equations and are summarised in the matrix equation as follows. [ The solution vector consists of the component vector of the ship 0 , which contains five variables related to the motion of the ship, and that of each tugboat, which contains seven variables related to the motion of the tugboat and pushing forces.
The coefficient matrices are composed of block matrices. is a 5 × 5 matrix that consists of the ship variables, and is a 7 × 7 matrix that consists of the variables of the i th tugboat. They are positioned diagonally in the coefficient matrix. Meanwhile, (5 × 7) and (7 × 5) are matrices representing the geometrical dependence between the ship and i th tugboat. The block matrices are represented as follows: It is assumed that there is no dependence between the tugboats. Therefore, the block matrix ( ≠ ≠ ) (7 × 7), which represents the direct interference from the j th pushing tugboat to the i th pushing tugboat, is set as a zero matrix.

Coupled manoeuvring model of a ship and tugboats: towing assistance
In this section, a situation in which multiple tugboats provide towing assistance to a ship is considered. Because the cable is always taut during towing, it is modelled as a pin joint element, and the tension force is supposed to act on the ship and tugboat in the direction of attraction. The right-hand coordinate systems assuming this situation are shown in Fig. 2, where the tension force acting on the cable of the i th tugboat is denoted as and its direction from the heading angle of the tugboat is denoted as . These are unknown variables that must be addressed. The spatial angle of this cable is denoted as , which is assumed to be constant owing to the taut cable. Regarding the geometrical variables, and represent the distance and plane angle from the midship of the ship to the i th towing point on 0 -0 0 0 , respectively.
is the longitudinal distance from the midship of the i th tugboat to the towing point on that tugboat, which is assumed to be at the centreline, that is, ( , 0) on -. The descriptions of the other variables are presented in Fig. 1.     The relationship between the acceleration terms individually defined in the space-fixed and hull-fixed coordinate systems is given by Eqs. (14) and (15) for the ship and Eqs. (16) and (17) for the i th tugboat, respectively.

= [̇̇̇̇̇̈]
The expression for is presented in Eq. (19). Thus, when n number tugboats tow the ship, the 5 + 7 matrix equation is solved. Regarding the constant vector, the component vector of the i th tugboat is expressed as follows.
The expression for is given by Eq. (21). Meanwhile, the coefficient matrix is composed of block matrices, where , , ⋯ are positioned diagonally in the matrix. is a 7 × 7 matrix composed of the variables of the i th tugboat, and 0 (5× 7) and 0 (7 × 5) are matrices representing the geometrical dependence between the ship and i th towing tugboat. These block matrices, excluding those of (Eq. (23)), can be expressed as follows. ( ≠ ≠ ) (7 × 7) is set as a zero matrix assuming no direct interference from the j th towing tugboat to the i th towing tugboat.

Manoeuvring model of a ship and tugboats while tugboats move freely
A situation in which the n number of tugboats manoeuvre freely near the ship is considered. Each vessel followed 3-DOF motions, and it was assumed that their motions did not interfere with each other. Consequently, the coefficient matrix consists of only a block matrix of diagonal elements.

Manoeuvring model of a ship and tugboats under various assistance modes
We consider a situation in which an arbitrary number of pushing tugboats and towing tugboats assist a ship cooperatively and free-manoeuvring tugboats are near for their backup. The block matrix corresponding to the state of each tugboat was selected as a module and combined to form a coefficient matrix representing the entire situation. The solution and constant vectors were constructed in a similar manner. This is a modular-type manoeuvring model that can deal with a variety of time-changing states of tugboats. In this chapter, the timing to switch modes is explained first, and examples of the modular-type model are shown next.

Between free-moving and pushing-assistance modes
The coordinate values of the waterline at the designed draft of each tugboat and ship are stocked during the simulation. The contact timing before the pushing-assistance mode is detected using the winding number algorithm [34], known as the algorithm for judging the inclusion of a point in a planar polygon. The winding number of a point with respect to a polygon with k vertices is expressed as follows.
Where is the angle formed by two lines drawn from a certain point to two adjacent vertices. Monitoring this number of the bow point of each tugboat with respect to the waterline of the ship, the moment when this number changes from 0 to 1 is considered as the contact timing and the coordinate values at the intersection is detected as the contact point. Meanwhile, the end of pushing-assistance is detected by monitoring the direction (sign) of the pushing force which is solved and updated time by time. When the tugboat begins to have a propulsion not repulsion, the mode switches to the free-moving one.
5.1.2 Between free-moving and towing-assistance modes Since in reality the cable length is adjusted by the tugboat's winch and the cable can be winded up at any timing, the start-timing of towing-assistance is arbitrarily determined by the author in the simulation. Meanwhile, during towing, the towing force on the tension cable always acts on both tugboat and ship as an attraction force. Therefore, the end of the towing-assistance or the timing of the cable looseness is detected by monitoring the direction of the towing force. The towing-assistance mode switches to the free-moving mode after the detection. Fig. 3 shows examples of the conversion of matrix equations obtained by changing the combination of modules. The block matrices and component vectors are coloured depending on the state of the tugboats; those of the pushing mode, towing mode, and independent freemanoeuvring mode are highlighted in light blue, pink, and yellow, respectively. The matrices and vectors attributed to the ship manoeuvring are highlighted in green, and the zero matrix is in grey. In addition, the matrices related to the coupling terms between a tugboat and ship is highlighted in mixed colours. Note the size of image of each block matrix corresponds to its number of rows and columns written in parentheses. When the first tugboat tows, the second tugboat pushes the ship, and the third tugboat moves freely, the matrix equation is presented in (a). It is converted into (b) if the cable of the first tugboat becomes loose. Furthermore, when the first tugboat joins to push the ship and the third tugboat also starts towing the ship, the

Trial simulation of assistance operations
Some cases were simulated just as examples to demonstrate multiple tugboats cooperated to assist the ship manoeuvring based on the proposed model. Under the quasi-steady-state assumption, the matrix equation, which was updated from time to time, was solved by the Gauss elimination method, and time integration was carried out using the Runge-Kutta method.

Outline of the subject ship and tugboat
A full-load bulk carrier with a real-scale length of 178 m was used for the subject ship. The experimental coefficients used in the manoeuvring model, such as the hydrodynamic force derivatives and crossflow drag coefficients, were obtained from [33]. As there are few published data on the hydrodynamics of a tugboat, we have used the data of a pusher boat [35], which has relatively similar dimensions to domestic standard tugboats. However, the crossflow drag coefficients, which were not published in above reference, were obtained by analogy. The principal dimensions of the subject ship and tugboat are listed in Table 1. The engine power of the tugboat was assumed to be 3000 ps, and the maximum thrust was set to 30 tf. In the simulation, the azimuth-type thruster was assumed to swing at a target angle of 2.32 deg/s.

Masaaki Sano
Mathematical modelling and simulation of cooperative manoeuvres among a ship and tugboats 6.2 Controlled thrust based on a pseudo-inverse matrix Compared to course-keeping control (e.g., [36]), there seem few references explicitly explaining a way of control for complex assistance operation by tugboats. In this study, the thrust is controlled in the following steps. It is tentative and sort of simple but gives a solution to this control problem.
First, the force ( , ) and moment required to move the ship as intended are defined based on the ship-fixed coordinate system and estimated by PID control with respect to the command speed, yaw rate, and heading angle of the ship. They should be equivalent to the total force provided by the n number of tugboats. In Eq. (54), ( _ , _ ) is the force exerted by the i th tugboat on the ship based on the ship-fixed coordinate system. In this study, this was determined using the Moore-Penrose pseudo-inverse matrix. Here, the matrix in Eq. (54) is defined as . Because the inverse matrix of cannot be defined, the pseudo-inverse matrix is defined instead.

= ( ) −
The solution vector that satisfies the minimum norm in a system with multiple solutions is calculated as follows. The magnitude of the command thrust vector i and its direction are calculated as follows.
Where takes 1 in pushing and -1 in towing because towing while moving backward. Note the thruster is assumed to reach to this command direction by a specified swing speed.

Simulation of a ship moving transversely while being pushed by multiple tugboats
We considered some tugboats with automatic control collaborating to push the ship in the lateral direction on a berthing stage. Approximately one ship length was assumed to be running, and the final target speed was set to 15 cm/s in the direction. Two scenarios were simulated, where 0 was the initial heading angle and 0 was the final one of the ship.  This scenario assumes that the ship is not anchored parallel to the berth ( 0 = −20 o ), and the three tugboats try to push the ship such that it is parallel to the berth while correcting its heading angle ( 0 = 0 o ). Fig. 5 shows the time series and trajectories, and we can see that the smooth movement of the ship is achieved by the coordinated-assistance movement of the three tugboats. Eventually, the ship reached a 0 o heading angle and moved in the lateral direction at the target speed.

Simulation of a ship moving transversely while being towed by multiple tugboats
Assuming an unberthing stage, we considered two tugboats collaborating to pull out the ship at the anchor and tow it at 15 cm/s in the lateral direction while moving backward. The two scenarios described below were also simulated. The initial angle of the thrust vector of each tugboat was set to 180 o .
6.4.1 The number of tugboats: = 2, 0 = 0 , and 0 = 0 o : In this scenario, two tugboats tow a ship while maintaining a constant heading angle of the ship. The bow-side tugboat in red was initially positioned at a slight deviation and the cable was not in parallel to the moving direction before towing. Fig. 6 shows the bow-side tugboat tried to adjust its position to tow the ship straight backward, and a small resultant oscillation was observed at the beginning of towing. It could be improved by updating the control method. But finally, the oscillation was converged and this scenario was realized by these tugboats.  Fig. 7 that they move dynamically to each other and undergo complex coordinated manoeuvres to achieve their goal. Eventually, the ship was towed in the lateral direction while changing toward the target heading angle. 6.5 Example of berthing simulation of the ship assisted by the two tugboats A simulation based on a simple scenario was performed to obtain an image of how multiple tugboats could assist in the berthing operation of a large ship. Total 15 scenes through the simulation were shot and shown in Fig.8. The following stages were considered for this scenario and the number in parentheses corresponds to the number of scene shots.


Stage A (1) / Both tugboats in free-moving mode: The two tugboats followed a ship that approached the berth area slowly. Because each tugboat and ship were assumed to be connected by a slack cable, they could manoeuvre freely.  Stage B (2, 3, and 4) / Both tugboats in towing-assistance mode: The cable became taut owing to the tugboats going astern. The timing of this astern was set by the author. The tension force of the tense cable acted on the ship as the breaking force until the ship stopped. The same control method as 6.2 was adopted.  Stage C (5, 6, 7, 8 and 9) / Both tugboats in free-moving mode: The tugboats moved around the ship to the port side. Then, they turned toward the ship and moved ahead. The waypoints and control of each tugboat were adjusted by trial and error in order to guide the tugboats to the ship.  Stage D (10) / One tugboat in pushing-assistance mode and another tugboat in freemoving mode: One tugboat contacted with the ship. Since the contact point was not predetermined, the touching place by the tugboat to the ship was considered as the contact point. The contact timing was detected by the way as explained in 5.1.1. After detection, this tugboat switched pushing-assistance mode and started pushing the ship. Meanwhile, another tugboat still run heading for the ship, which was considered in free-moving mode.  Stage E (11, 12, 13, 14 and 15) / Both tugboats in pushing-assistance mode: After the late tugboat contacted the ship, the two tugboats started pushing the ship in cooperation so as to push the ship slowly while adjusting it parallel to the berth. They were controlled by the method of 6.2. The simulation continued until the ship touched to the berth.
During the berthing assistance, we can verify that the tugboat manoeuvres continuously switch to different modes, that is, free-moving, towing and pushing. In particular, we see that the contact judgement was effective and that these tugboats started pushing the ship in cooperation after contact smoothly. The transition of the mathematical model, which was updated every stage, is shown in Fig. 9.

Fig. 9
Transition of the manoeuvring mathematical model corresponding to each stage

Conclusions
In the future, intelligent tugboats are expected to automatically provide manoeuvring assistance to large ships. This motivated us to develop a new mathematical model as a plant model of the entire system for berthing situation where multiple tugboats and a ship work together in a variety of modes in the same field. For this purpose, the dynamics of all tugboats and the ship were solved considering the coupling effect on their motions. It would be useful for building and verifying the control system of tugboats. The main contributions of this study are summarised as follows.  The manoeuvring mode of a tugboat was classified into three modes: pushing assistance, towing assistance, and free moving. The coupled motion of a ship with multiple tugboats under each mode was individually formulated, in which the pushing and towing forces acting between them are determined by their relative manoeuvring motions. These forces also influence their manoeuvring motions as well. The equations of motion were organized in a matrix. Some examples were simulated based on the proposed mathematical model. The control method adopted here was tentative just to demonstrate the cooperative manoeuvres and check the behaviours of each player. The relevant conclusions are summarized as follows.  Pushing and towing assistances were simulated. We verified that the multiple tugboats flexibly adjusted their heading and positions, which resulted in the adjustment of the pushing or towing force to move the ship as requested.  In the berthing simulation, it was shown that multiple tugboats could cooperate to assist the ship manoeuvres. We verified that the proposed model could handle various situations according to their assistance modes.
Because a tugboat has a high degree of freedom of movement, its control must be complex. In reality, as skilled experts manoeuvre tugboats carefully based on their experiences, the control method should be sophisticated in order to be able to cope with more realistic situations. It will be continuously studied using this simulation platform.