Control of Dynamic Positioning System with Disturbance Observer for Autonomous Marine Surface Vessels

The main goal of the research is to design an efficient controller for a dynamic positioning system for autonomous surface ships using the backstepping technique for the case of full-state feedback in the presence of unknown external disturbances. The obtained control commands are distributed to each actuator of the overactuated vessel via unconstrained control allocation. The numerical hydrodynamic model of CyberShip I and the model of environmental disturbances are applied to simulate the operation of the ship control system using the time domain analysis. Simulation studies are presented to illustrate the effectiveness of the proposed controller and its robustness to external disturbances.

In view of the maneuvering difficulties caused by the weight of a ship, it is not an easy task to improve the quality of navigation, especially for ships moving at low speed (called dynamic positioning). Designing an efficient dynamic positioning (DP) system for a marine vessel is a challenging practical problem. The performance and robustness of the DP system is essential for the success of the mission. In dynamic positioning systems, the main goal is to keep the marine vessel in a steady position and at a constant heading (direction) in the horizontal plane or to follow the target trajectory using only hull-mounted thrusters. The first generation of dynamic positioning systems comes from the early 1960s, when drilling began to be performed at very great depths. The first vessel equipped with a dynamic positioning system was Eureka, owned by Shell Oil Company, which entered into operation in 1961 [46]. Currently, dynamic positioning systems are used on various types of ships to perform many marine tasks, such as hydrographic surveys, marine construction, wreck research and geodesy. In the offshore oil and gas industry, many tasks can only be performed with the assistance of DP systems. This refers to the operation of service vessels, rigs and drilling vessels, shuttle tankers, cable and pipe laying units and floating production storage and offloading (FPSO) units.
The first implemented dynamic positioning systems made use of PID (Proportional-Integral-Derivative) controllers. To counteract the excessive activity of thrusts associated with wave-caused hull motion, the controllers used cut-off filters in a cascade arrangement with low-pass filters [47]. The improvement in the quality of control system operation took place after the application of more advanced control techniques based on the optimal control theory and the Kalman filter theory [48][49][50][51]. The major disadvantage of this approach was that the kinetic equations of motion had to be linearized for certain conditions. For each linearized equation, new gains were computed for the Kalman filter and for coupling, and then these gains were modified online by gain scheduling. In the 1990s, introducing nonlinear observers and feedback control theory to the designs of dynamic positioning systems resulted in the removal of assumptions related to linearization [52][53][54].
In dynamic positioning systems, multivariable controllers usually determine the commanded forces and torque, which must then be generated by (obtained from) thrusters installed on the ship. For overactuated marine vessels, control allocation is a vital part of the DP system. Improper allocation may lead to degraded control performance, lower energy efficiency and the increased wear and tear of the actuators. There is a rich literature regarding control allocation for marine surface vessels, commonly referred to as thrust allocation [106][107][108][109][110][111][112][113][114][115][116][117][118][119][120]. In-depth reviews of the literature are given in [121,122].
A vessel operating in the ocean is subjected to disturbances caused by waves, wind and sea currents that cause the vessel to deviate from its desired position and direction. Hence, disturbance damping is one of the key problems in the design of DP control. The ocean disturbance can be divided into low-frequency (LF) disturbance caused by second-order waves, sea currents and wind and wave frequency (WF) disturbance caused by first-order waves. The low-frequency disturbance causes the ship to drift, while the wave frequency disturbance causes it to oscillate. The compensation of wave frequency disturbances would wear out the marine actuators and increase fuel consumption. On the other hand, there is no need to compensate for WF disturbances because they cause only oscillatory movements of the ship [69]. These motions should be filtered off by wave filtering algorithms from the vessel position and heading measurements before passing them to the DP control system. Several wave filtering techniques have been proposed [123][124][125][126]. Therefore, in this article, only the low-frequency components of environmental disturbances are considered in the DP control design.
The control objective in this paper is to design a ship motion control algorithm in dynamic positioning using the backstepping method with a disturbance observer of unmeasured disturbances affecting the ship's hull.

Formulation of the Problem of Steering the Ship at Low Speed
The motion of a ship sailing on the water surface in the horizontal plane is analyzed in three degrees of freedom. This motion is described using two coordinate systems, limited to two dimensions: the inertial frame (X N , Y N ) associated with the sea map and the body-fixed reference frame (X B ,Y B ) associated with the moving ship ( Figure 1). Physical quantities assumed as state variables for the ship moving in the horizontal plane can be grouped into two vectors: η = [x, y, ψ] T and ν = [u, v, r] T , where (x,y) are the coordinates of ship's position, ψ is the ship's heading, (u, v) are the linear velocity components of ship's motion in the surge and sway directions, and r is the yaw rate [127]. The velocity vector determined in the inertial frame (X N , Y N ) is related to that determined in the body-fixed reference frame (X B , Y B ) by the following kinematic relationship: where R(ψ) ∈ R 3×3 is the rotation matrix by angle ψ, determined from the relationship The properties of the rotation matrix given by Formula (2) are as follows where S ∈ R 3×3 is the skew-symmetric matrix The mathematical model of the dynamics of a ship sailing on the surface of the sea and ocean in the presence of environmental disturbances is described as follows [127]: where M ∈ R 3×3 is the inertia matrix, C ∈ R 3×3 is the matrix of Coriolis and centripetal terms, D ∈ R 3×3 is the damping matrix, and τ = [τ x , τ y , τ n ] T is the input control vector consisting of the surge force τ x , sway force τ y and yaw moment τ n produced by propellers and thrusters installed in the ship's hull. Unmodeled external low frequency (LF) forces and moments due to wind, currents and waves are connected together into an Earth-fixed constant (or slowly varying) bias term b Here, it is assumed that the changing rate of the bias is bounded, where C d is a nonnegative constant. The above assumption is reasonable because the environmental energy applied to the vessel is limited. In Equation (5), τ d represents external disturbances acting on the vessel in the bodyfixed reference frame, given as The inertia matrix M ∈ R 3×3 , which includes the hydrodynamic added inertia, can be written as [127] where m is the vessel mass, I z is the moment of inertia about the fixed z-axis of the vessel, and Xu, Yv, Yṙ, Nv and Nṙ are hydrodynamic derivatives. Zero-frequency masses are added to the surge, sway and yaw due to accelerations along the relevant axes. For control applications, which are restricted to LF motions, the wave frequency independence of the added inertia (zero wave frequency) can be assumed. This implies thatṀ = 0. The matrix of Coriolis and centripetal terms has the form For a straight-line stable vessel, D ∈ R 3×3 is a positive damping matrix due to linear wave drift and laminar skin friction. The linear damping matrix is defined as [127]

Control Algorithms in the Dynamic Positioning System
The control objective in this paper is to design a DP control system for a ship with unknown time-varying disturbances, so that the vessel's actual position (x,y) and heading converge to the desired values Assumption 1. The desired smooth reference signal η d is bounded and has bounded firstη d and secondη d derivatives. This means that the functions describing the ship's position and direction, as well as their derivatives of the first and second order, are limited.
The considered movement of the vessel is executed at low speed in the system shown in Figure 2. The input to this control system is the vector η r with reference coordinates of the ship's position and direction. The smooth and bounded desired trajectories η d with their firstη d and second-orderη d derivatives needed for the controller were generated by a second-order filter: where the reference η r is the operator input, ζ i is the relative damping ratio, and ω ni is the natural frequency. Notice that lim andη d andη d are smooth and bounded derivatives for steps in η r . The vessel is assumed to be overactuated. The trajectory of motion of such a vessel depends on the operation of all actuators installed in its hull. The controller's task is to determine the forces and moments to be applied to the ship's hull. This also requires the use of an appropriate system for the distribution of forces and torque determined by the controller.

Thrusters
Vessel dynamics

Backstepping Method with Disturbance Observer
The control algorithm used in the dynamic positioning system was derived using the backstepping method and assuming that the entire plant state vector is known. The vector of control forces τ c (t) was designed in such a way as to ensure that the state variables in vectors η(t) and ν(t) remain constrained and that the position and course are asymptotically convergent to their set constant values η(t) → η d (t) with ν(t) ≈ 0 for t ≥ 0. The classic method of backstepping is described in [72]. The desired signals required for control are represented by the given position and direction vector η d = [x d , y d , ψ d ] T and its firstη d and secondη d derivatives. It is assumed that all desired signals related to the ship position (x d , y d ) and heading ψ d are limited.
The control deviations related to the given position and direction vector η d and the velocity vector ν were defined as where ϑ 1 ∈ R 3 is the stabilizing function, which is the desired virtual control. Determining the two-sided derivative from Equation (12), and substituting relation (1) and that determined from Equation (13) into this derivative, we obtaiṅ The stabilizing function ϑ 1 was assumed as where K 1 = K T 1 > 0 is the diagonal positive definite gain matrix. The stabilizing function ϑ 1 , which is the desired virtual control for vector ν, was determined in relation to the Lyapunov function in the form V 1 = 0.5z T 1 z 2 . Substituting relation (15) into Equation (14) and using the relation R(ψ)R −1 (ψ) = I, we obtaiṅ Determining the two-sided derivative from Equation (13) and substituting the relation determined from Equation (5) into it, we obtaiṅ whereθ 1 , determined based on Equation (15), takes the forṁ With regard to the Lyapunov function in the form V 2 = V 1 + 0.5z T 2 Mz 2 , the following control law was adopted: where K 2 = K T 2 > 0 is the positively defined gain matrix, and the vectorb contains estimates of the parameters of the external bias term b describing external disturbances acting on the vessel. The disturbance observer for the bias vector b was constructed using the exponential convergent observer from [128] T is the estimate of the bias term, K 0 is the 3 × 3 positive definite symmetric observer gain matrix, and θ is the 3 × 1 intermediate auxiliary vector.
Considering the ship dynamics given by Formula (5) and the desired trajectory η d , which is smooth and limited, the information about all ship states x is provided. In this case, the control law is described by Formula (19). The law of adaptation of unknown parameters related to environmental disturbances is described by Formula (20) and has zero values as initial conditions. The entered design parameters These conditions mean that the entire closed control system is stable, and consequently the signals z 1 and z 2 have finite values.

Nonlinear PID Controller
The controller used to compare the obtained results of simulation tests was the nonlinear PID controller for DP systems, described in detail in [18]. The algorithm of this controller is given by the formula where R is the rotation matrix (2), S is the skew-symmetric matrix (4), M is the inertia matrix (7), C(ν) is the Coriolis and centripetal matrix (8), D is the matrix of hydrodynamic damping (9), and K P , K I and K D are the matrix gains of the nonlinear PID controller. In this algorithm, the position and orientation errors η e are determined in the inertial frame (X N , Y N ), while the velocity error ν e and the acceleration errorν e are determined in the body-fixed reference frame (X B , Y B ) associated with the moving ship The gain matrices of the nonlinear PID controller are positively defined:

Unconstrained Control Allocation
The vectors τ c or τ PID of the desired forces and moments determined by the DP controllers were divided by the allocation system into the commanded values for the actuators installed in the ship's hull (Figure 2). It is assumed that the ship is fitted with q propulsion devices located at positions where l i,x , l i,y are the distances of the propulsion devices from the origin of the coordinate system associated with the moving ship. These devices can provide the thrust force T i in the direction defined by the angle α i . The azimuth thrusters have the angle α i fixed in a certain direction and produce the following contributions to the generalized forces acting on the ship [120]: where L i = l i,x sin(α i ) − l i,y cos(α i ).
The sum of the generalized forces acting on the ship's hull from all installed propellers is given by the formula where In Formula (32), τ c ∈ R 3 is the input control vector, and u ∈ R q is the vector of the set forces for each actuator.
In control allocation, the optimization of a given quality indicator, such as the minimum energy expenditure for control, is often performed. The problem of optimal control allocation has been considered as a minimization problem with constrained equality [122] min u u T Wu in relation to where W ∈ R qxq is the weight matrix.
Considering the Lagrangian function, the formulated optimization problem (34) can be written as the unlimited minimization problem in the form where λ ∈ R 3 is the vector containing Lagrange multipliers. The Karush-Kuhn-Tucker (KKT) conditions provide the necessary conditions for the optimal solution After solving the system of Equations (36) and (37), the sought vector of Lagrange multipliers is obtained as Substituting Equation (38) to (36) gives When the matrix BW −1 B T is non-singular, the optimization problem (34) can be solved by finding the solution to the linear equation described by Formula (39).

Simulations
To demonstrate the effectiveness of the developed control algorithm, selected simulation tests were carried out with the dynamic positioning system.

Ship Model
The mathematical model of CyberShip I was chosen as the ship model for determining and testing the control system at low speed. This model was developed in the Department of Engineering Cybernetics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway. The physical model of this ship sails in the Marine Cybernetics Laboratory, NTNU (http://www.ntnu.edu/imt/lab/cybernetics (accessed on 10 September 2021)).
CyberShip I is a thruster-controlled model of an offshore supply vessel, made at a scale of 1:70. Its mass is m = 17.6 kg with length L = 1.19 m. The centre of gravity is located at x G = −0.04 m aft of the midship. This point was assumed to be the origin of the body-fixed coordinate system. In its general form, the mathematical model of CyberShip I is described by Formula (5). Based on hydrodynamic methods and system identification, model parameters for 3 degrees of CyberShip I freedom of motion were found [17,129,130]. The inertia matrix including zero-frequency hydrodynamic added inertia is as follows: The matrix of linear interactions related to hydrodynamic damping is defined as In real-time control systems, the mathematical model of the plant does not fully reflect its dynamics. To bring the simulation tests closer to the conditions in real-time control systems, new parameters of matrices M and D were determined for the model of ship dynamics. In Figure 2, the model of the plant is included in a block labeled "vessel dynamics". The principle of determining new parameter values was based on the values contained in the model described by matrices (40) and (41). The correction values ∆ were determined randomly and could at most amount to ± 50% of the parameter values included in these matrices. In this way, the following matrix values were determined: The model of Cybership I is equipped with four rpm-controlled thrusters with independently controllable azimuth angles α i . The thrusters are controlled by the rotational speed ω i . Figure 3 shows the location of the actuators installed on CyberShip I.
where B ∈ R 3×4 is the thruster configuration matrix described by Formula (33), in which while the components of the thrust vector of thrusters T ∈ R 4 are described by the formula where k 1T = k 2T = k 4T = 3.125 · 10 −3 , k 3T = 2.5 · 10 −4 . The thrusts generated by the azimuth thrusters are [131] T max = −T min = 0.8 0.8 0.1 0.8 T The inverse characteristics of the thrusters, determined based on Formula (46), have the form The limits imposed on the rotational speeds of thrusters were determined based on Formula (48) and related parameters as

Parameters of Control Systems
The mathematical model of the control system is shown in Figure 2. Here, the vessel model marked as "vessel dynamics" and described by Formula (5) contains the matrix coefficients M and D, given by Formulas (42) and (43), respectively. The coefficients of matrix C are determined from Formula (8). The tested control algorithms described by Formulas (19) and (22) also contain matrices M, C and D. In this case, the values of coefficients described by Formulas (40) and (41) were adopted. Other parameters of the controller determined by the backstepping method (19) are as follows: K 0 = diag{2, 2, 2}, K 1 = diag{0.05, 0.05, 0.05}, K 2 = diag{1, 1, 1}. The parameters of the gain matrices K P , K I and K D for the nonlinear PID controller (22) were determined based on the mathematical model of Cybership I, given by Formulas (40) and (41), using the Particle Swarm Optimization algorithm and assuming no external disturbances. In this way, the following parameters were obtained for the nonlinear PID controller: K P = diag{1. 25, 7.4, 9.95}, K I = diag{0.035, 0.343, 1.0}, K D = diag{3.57, 6.88, 6.06}.
The parameters of the reference model (10) were ζ i = 1, ω ni = 0.08, where i = 1, 2, 3. The maximum values in Formulas (19) and (22) were imposed as limits on the control signal generated by the DP controllers (51). In the control allocation system, the values of the weight matrix W were W = diag{1, 1, 10, 1}.
The  Figure 4, from which it can be seen that in steady states the disturbance estimates coincide. There are, however, two transients: the first after switching on the control and the second after starting to change the vessel position. The zoomed time segments representing transient states are shown in Figure 5. The first transient lasts about 8 s, while the second is longer and lasts about 80 s. In this latter transient state, the deviations from the real value are very small. Figure 6 shows that the proposed controller is able to follow the desired reference trajectory. The time histories showing the desired and actual ship positions (x, y) and courses ψ can follow the desired trajectory η d = [x d , y d , ψ d ] T with good precision. The deviations from the desired values are shown in Figure 7. It is noteworthy that greater values of deviations were recorded in the initial period of time and during the stabilization of the set exchange rate. This is mainly due to the low power of the actuators installed in the CyberShip I hull to generate the angular moment τ n . Figure 8 presents the time histories of ship velocity changes in the surge u, sway v and yaw r directions, while the next graphs show the time histories of other quantities occurring in the control system, such as the desired forces τ x , τ y and moment τ n , which are the output signals from the controllers (Figure 9) and the commanded rotational speeds of the thrusters installed in the vessel's hull ( Figure 10).

Case 2-Stochastic and Time-Varying Disturbances
In the next tested case, the disturbances affecting the ship with stochastic and timevarying forces and moment described by Formula (53) The simulation test results obtained for this case are shown in Figures 11-17. The time histories of environmental disturbances b and their estimatesb are shown in Figure 11, from which it can be seen that, in the steady states, the disturbance estimates coincide. There are, however, two transients: the first after switching on the control, and the second after starting to change the position of the vessel. The zoomed time segments representing the transient states are shown in Figure 11. The first transient lasts about 10 s, while in the second transient, small deviations are observed in estimatesb 1 andb 3 . There are no deviations in the estimateb 2 . The time histories of the desired and actual ship positions (x, y) and ship course ψ shown in Figure 13 follow the desired trajectory η d = [x d , y d , ψ d ] T with good precision. The deviations from the desired values are shown in Figure 14. Figure 15 presents the time histories of ship velocity changes in the surge u, sway v and yaw r directions. The next graphs show the time histories of other quantities occurring in the control system, such as the desired forces τ x , τ y and moment τ n , which are the output signals from the controllers (Figure 16), and the command rotational speeds of the thrusters installed in the vessel's hull ( Figure 17).    A quantitative comparison of the quality of the performance of the two considered controllers acting in the presence of constant and changing disturbances is given in Table 1, where x e = x d − x, y e = y d − y are the errors between the desired and current position, ψ e = ψ d -ψ is the difference between the desired and actual heading of the vessel and t f = 250 s. The results presented in Table 1 clearly show that the controller τ c (20) has better control quality than the nonlinear τ PID controller (22).

Discussion
Numerous institutions and universities are becoming increasingly interested in the development of control algorithms for autonomous marine surface vessels. One of the important tasks is the automation of the process of controlling the surface vessel's motion for the entire voyage, starting from the departure port and ending at the destination port. In this case, the desired route of the system consists of a number of different-type segments, and thus it may be necessary to use different controllers at different path stages. Ship navigation on the desired route defined in the above way requires the design of a control system that is capable of executing various tasks, such as ship undocking and docking, maneuvering in the port area, movement along the desired route with transit speed and stopping on the route. Such a solution was presented in [132].
The control algorithm presented in this article is planned to be used in a multioperational system to control the motion of a ship on those segments of the voyage route where the vessel will sail in dynamic positioning mode; i.e., in ports and very narrow navigation canals, on access routes to the port and in the maneuvers performed to reach the docking place.

Conclusions
The article presents a ship motion control system with a disturbance observer for the dynamic positioning of a fully actuated autonomous marine surface vessel in the presence of uncertain time-variant disturbances due to wind, waves and ocean currents. Both the Coriolis and centripetal matrix and the linear damping matrix are considered in the mathematical model of the vessel. The control strategy is introduced by the backstepping technique with a disturbance observer used to compensate for uncertainties associated with the disturbances. The simulation tests carried out on the model of a sea-going vessel have shown that the designed controller is effective in compensating for external disturbances. The proposed control system is characterized by a good quality of work both during the stabilization of the fixed position of the ship and when moving it to a new position.
The unconstrained control allocation was used to allocate the desired control to individual actuators. The applied allocation method allows for the distribution of the desired forces and moment determined by the DP controller into any number of thrusters installed in the ship's hull with a fixed thruster, which is a non-rotatable device, and its orientation angle α cannot be changed.  Data Availability Statement: Data available on request due to restrictions; e.g., privacy or ethical considerations.

Conflicts of Interest:
The authors declare no conflict of interest regarding the publication of this paper. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

MASS
Maritime Autonomous Surface Ship DP Dynamic Positioning FPSO Floating Production Storage and Offloading LQG Linear Quadratic-Gaussian NTNU Norwegian University of Science and Technology PID Proportional-Integral-Derivative