Design and Assessment of ADRC-Based Autopilot for Energy-Efficient Ship Steering

Abstract

This paper looks at a typical problem encountered in the process of designing an automatic ship’s course stabilisation system with the use of a relatively new methodology referred to as the Active Disturbance Rejection Control (ADRC). The main advantage of this approach over classic autopilots based on PID algorithms, still in the majority, is that it eliminates the tuning problem and, thus, ensures a much better average performance of the ship in various speed, loading, nautical and weather conditions during a voyage. All of these factors call for different and often dynamically variable autopilot parameters, which are difficult to assess, especially by the ship’s crew or owner. The original result of this article is that the required controller parameters are approximated based on some canonical model structure and analysis of the hydrodynamic properties of a wide class of ships. Another novelty is the use of a fully verified, realistic numerical hydrodynamic model of the ship as a simulation model as well as a basis for deriving a simplified model structure suitable for controller design. The preliminary results obtained indicate good performance of the proposed ADRC autopilot and provide prospects for its successful implementation on a real ship.

1. Introduction

A good design of a ship and her steering control system is crucial for the safety and efficiency, or fuel economy, of the operation (navigation) process. Sometimes, and to some extent, a good control system can even enhance the actual poor design of the ship’s hull, propeller and rudder, especially in heavy weather conditions.

Ship steering, realized either by automatic (autopilot) or manual (human) control, is twofold. First, it consists of course-keeping that is following a predetermined heading (most frequent) or track (as in modern integrated navigation systems). Second, it includes course-changing manoeuvres to reach the next straight-line section (leg) of the route according to the voyage plan or as a part of collision avoidance. In both of these tasks, the steering optimization by the proper selection of an autopilot structure/design and the tuning of its parameters can render better energy efficiency. For example, the problem with the often inappropriate tuning of the classic PID (Proportional–Integral–Derivative) regulators is the cause of significant energy (fuel) losses resulting from excessive ship’s yawing as well as excessive wear of the steering gear. Generally a properly tuned autopilot throughout the whole ship’s voyage is a crucial factor in the energy saving issue [1].

Real-life conditions are characterized by various types of uncertainties, such as inaccuracies in the system model, the occurrence of random processes such as winds, waves or currents, and other exogenous effects. Therefore, when designing the ship’s control system, it is necessary to use techniques that take into account the nonlinear effects of the process and the parametric uncertainty of the ship’s model [2].

There are many works on ship autopilots. Apart from the classic PIDs, there are more and more works on autopilots based on modern control algorithms, e.g., nonlinear [2] robust, adaptive [3,4] or intelligent [5]. Among those that deserve attention, there are also autopilots based on the ADRC algorithm [6].

However, their simulation verification (except classical sea trials) has usually been performed on simplified and not fully validated models [6]. Moreover, they often require the tuning of a number of parameters [7,8] (which often is only possible offline using AI-type algorithms) or some detailed knowledge about the ship, obtained, for example, from sea trials or identification procedures [9].

A comparative analysis of control performance obtained via the ADRC versus a traditional autopilot based on the PID algorithm is discussed in [10]. However, the data (ship model parameters) used there relate to the training ship Blue Lady, a physical model of the tanker made in 1:24 scale [11,12]. Therefore, the results obtained in the present research are practically more representative.

The problems of ship steering control design are wide and complex [13]. We can often define them in terms of classical optimisation, specifying an objective function as spread on the controller’s model structure and parameters and minimising/maximising it to obtain the ‘optimal’ design. Such an optimality appears, however, to be mostly local, for very specific conditions of the ship and environment, and the control often fails as it is very sensitive to model uncertainty and disturbances, as with PID controllers and other controllers designed using so-called optimal control techniques [13,14,15]. This is why other approaches are being sought, e.g., a disturbance observer-based control (DOBC) [13,16] and the relatively new ADRC method being investigated in the present work. According to its specific type, the rather fixed structure of an ADRC controller (additionally of partly non-parametric definition) cannot easily undergo an optimisation in a mathematical sense.

The stipulated energy efficiency for an ADRC controller, for all ship operational conditions, is treated hereafter qualitatively rather than quantitatively. This is also partly due to this being the preliminary stage of our ADRC development in the ship steering domain, and there being limited performance test results. At this moment, the ADRC is expected to be at least equivalent to the PID controller but without its ‘local’ features, and, hence, (as presenting with no need of tuning) can be deemed more effective over the whole ship’s operational lifetime or mission profile. For industrial implementation, as the next step to our commonly accepted MATLAB/Simulink-based studies in this paper, a more comprehensive and quantitative comparison of various controllers, and for various ships, is being scheduled in the future. This would involve different energy-efficiency indicators and hardware-in-the-loop (HIL) studies, especially on the full-mission ship-handling simulator at the authors’ premises (Maritime University of Szczecin, Poland), and full-scale trials.

The Active Disturbance Rejection Control (ADRC) method, originally proposed by Han in 1989 [17,18], is an innovative approach to the control of a plant whose dynamic properties are uncertain and affected by external disturbances. Han’s concept is based on the use of an extended state observer in the control loop (ESO).

The observer, apart from tracking the values of state variables, supports the estimation of a certain virtual variable consisting of combined external disturbances and components resulting from the uncertainty of the plant parameters and its nonlinearities.

The ADRC methodology has been gaining in popularity [19,20]. It should be noted here that, although based on the concept of ESO, it incorporates other ingredients, such as the transient profile generator, noise-tolerant tracking differentiator or special nonlinear feedback [17]. Regardless, it seems that the estimation of total disturbance and rejection of the same, discussed in this paper, are of key importance here.

Considering that complex control systems can be typically reduced to the first or second order, the ADRC facilitates the procedure through a joint treatment of external disturbances and uncertainty modelling. Since little knowledge is required of the process and the disturbance model, the ADRC can exercise effective control over a complex system in the absence of a detailed and accurate mathematical model. Owing to its simplicity, the technique is an attractive alternative to much more advanced model-based approaches, such as designs based on adaptive and/or robust control methodologies [4]. Thereby, practitioners looking for more efficient algorithms than the standard PID can consider using the ADRC.

The paper comprises six sections concluded with a discussion. In Section 2, the principal concept of the ADRC control design is explained, followed with a description of the mathematical model of the ship’s dynamics in Section 3. In Section 4, zero dynamics and a model suitable for the controller design are discussed. In Section 5, the authors explain how to design an autopilot using the ADRC technique. The final Section 6 and Section 7 include simulation tests and their results, followed with a discussion.

2. The Principal Concept of the ADRC Control

Let us explain the main idea behind the ADRC control technique. Rather than formal proof, we focus on an intuitive analysis of the system output stabilization. For simplicity and clarity, yet without losing generality, the following second order, SISO, nonlinear affine system should be considered

$$\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=f(x_1,x_2)+cu+d(t)\\ y=x_1\end{array}$$

(1)

where

$x_1,x_2$—state variables,

u—control variable,

y—system output,

f and c—known function and constant respectively,

d—known disturbance.

It should be borne in mind that these considerations apply as well to a general nonlinear model in the canonical (Brunovsky) form.

In order to stabilize the system output at the assumed set-point $y_d$, a simple feedback linearizing controller can be applied

$$u=\frac{-f(x_1,x_2)-d+K_p(y_d-x_1)+K_d{x}_2}{c}$$

(2)

where $y_d$ represents the desired output, while $K_p$ and $K_d$ are the controller gains.

To make the calculation more realistic, let us assume that the function f is unknown and the parameter c is uncertain, i.e.,: $c=c_0+\mathsf{\Delta }c$, where $c_0$ is its known estimate (cf. Section 3, Section 4 and Section 5), while $\mathsf{\Delta }c$ represents its range of uncertainty. Additionally, the system is exposed to the unknown exogenous load disturbances d(t).

On the basis of the above, the following is obtained

$$\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=f(x_1,x_2)+\mathsf{\Delta }cu+d+c_{0}u\end{array}$$

(3)

For the purpose of designing a stabilizing controller, all system uncertainties shall be grouped into a single expression $F(t)=f(x_1,x_2)+\mathsf{\Delta }cu+d$, which represents generalized or total disturbances.

In the next step, a subsequent state variable $x_3(t)=F(t)$ is introduced. Assuming that $\dot{F}(t)=G(t)$ where G(t) is unknown, we obtain a formally equivalent and always observable system

$$\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=x_3+c_{0}u\\ {\dot{x}}_3=G(t)\end{array}$$

(4)

If all the state vector components are available, the controller of the linear system above (4) (including the ‘virtual’ disturbance G(t)), can be expressed with the following formula (cf. (2))

$$u=\frac{-x_3+K_p(y_d-x_1)+K_d{x}_2}{c_0}$$

(5)

Note that with the controller expressed with (5), the disturbance $F(t)=x_3(t)$ is cancelled, and the remainder is a PD-type algorithm. This way, our system, written in the form (4), (5) has a completely known form

$$\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=u\end{array} \mathrm{where} u=K_p(y_d-x_1)+K_d{x}_2$$

(6)

Hence, its behaviour can be determined precisely via proper selection of the controller gains (e.g., by the pole-placement or LQR technique [21]) under laboratory conditions. Thus, in the further part we assume that $K_p$ and $K_d$ are known (cf. (38)). The key now is to have all the state vector components available.

The next and last step of the ADRC methodology is the design of an observer [21]. For this purpose, (4) is expressed in the matrix form. The observer is, obviously, designed for a disturbance-free system, i.e., a system expressed as follows

$$\underset{\dot{x}}{\underbrace{\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\end{array}\right]}}=\underset{A}{\underbrace{\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]}}\underset{x}{\underbrace{\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]}}+\underset{B}{\underbrace{\left[\begin{array}{c}0\\ c_0\\ 0\end{array}\right]}}u,\hspace{1em}y=\underset{C}{\underbrace{\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}}\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$$

(7)

In the formalism above, F(t) is perceived by the state observer as a constant state variable $x_3(t)$ (cf. the third equation, where ${\dot{x}}_3=0$). Hence, the reconstruction of exact state can take place only in the steady-state $x^s={\left[x_1^s,\hspace{0.17em}{x}_2^s,\hspace{0.17em}{x}_3^s\right]}^T={\left[y_d,\hspace{0.17em}0\hspace{0.17em},f(y_d,0)+\mathsf{\Delta }c{u}^s+d\right]}^T$, ${\dot{x}}^s=0$.

There is a permanent transient process for non-constant disturbances (such as in the event of a step change of the set-point $y_d$), where the observer state component ${\widehat{x}}_3(t)$ follows up the varying total disturbance $x_3(t)$. Implementation of the controller (5) will obviously require a replacement of the original unmeasurable states $x_i$ with their estimates ${\widehat{x}}_i$ (i = 2,3) using the principle of certainty equivalence.

To summarise the concept of the ADRC controller, let us first write down the general formula for the observer in Kalman form

$$\dot{\widehat{x}}=A\widehat{x}+Bu+L(y-C\widehat{x})$$

(7a)

where L is the observer gain matrix (cf. Section 5).

Now we can obtain its Laplace transform form (for scalar inputs)

$$\widehat{X}(s)=H(s)\left[\begin{array}{c}U(s)\\ Y(s)\end{array}\right]$$

(7b)

where the matrix transfer function is

$$H(s)={\left[\begin{array}{c}{\left(sI\text{-}A+LC\right)}^{-1}B\\ {\left(sI\text{-}A+LC\right)}^{-1}L\end{array}\right]}^T$$

(7c)

and I denotes 3 × 3 identity matrix.

Finally, by using (5) we can write the formula for the whole ADRC controller

$$u=\frac{-{\widehat{x}}_3+K_p(y_d-x_1)+K_d{\widehat{x}}_2}{c_0}$$

(7d)

where the state estimates are taken from (7b).

A natural attempt to formally compare ADRC to a PID controller encounters certain difficulties and methodological challenges. There are some similarities, e.g., ADRC contains PD gains. On the other hand, ADRC does not require direct integral action (which is a weakness of PID). Moreover (as shown in the article), ADRC tuning is global rather than local and, thus, seriously simplified.

3. Ship Manoeuvring Dynamics

3.1. Nonlinear Model

The ship’s motion in the horizontal plane is represented by two vectors of state variables: the position vector [x₀, y₀, ψ]^T, where (x₀, y₀) and ψ denote earth-based linear position coordinates and angular position (heading), accordingly, and the velocity vector in the ship’s body axes ν = [u, v, r]^T, where (u, v) and r are linear velocity components (longitudinal and lateral, or surge and sway) and angular velocity (rate of turn or yaw). The origin is taken at the intersection of the midship section and centreplane, with the ship’s centre of gravity coinciding with this point (usually, at max., a few % of the ship’s length for merchant ships, and of negligible effect in practical manoeuvring). The mentioned former vector is related to the latter one by the kinematic law, whereas the latter is governed by the Newtonian laws of dynamics for a rigid body in fluid (where we adopt Fossen [22] notation)

$$\left[\begin{array}{c}\begin{array}{c}{\dot{x}}_0\\ {\dot{y}}_0\end{array}\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}u\cos \psi -v\sin \psi \\ u\sin \psi +v\cos \psi \end{array}\\ r\end{array}\right]$$

(8)

$$\mathit{M}\dot{\mathit{v}}+\mathit{C}\left(\mathit{v}\right)\mathit{v}+\mathit{D}\left(\mathit{v}\right)\mathit{v}=\mathit{\tau }$$

(9)

where M—mass/inertia matrix (including added masses), C—centripetal matrix (also including added masses), D—damping matrix (velocity-related hydrodynamic hull forces), τ—control (external force) vector

$$\mathit{M}=\left[\begin{array}{ccc}m-X_{\dot{u}}& 0& 0\\ 0& m-Y_{\dot{v}}& 0\\ 0& 0& J_z-N_{\dot{r}}\end{array}\right]$$

(10)

$$\mathit{C}=\left[\begin{array}{ccc}0& 0& -\left(m-Y_{\dot{v}}\right)v\\ 0& 0& \left(m-X_{\dot{u}}\right)u\\ \left(m-Y_{\dot{v}}\right)v& -\left(m-X_{\dot{u}}\right)u& 0\end{array}\right]$$

(11)

$$\mathit{D}=\left[\begin{array}{ccc}-d_{11}\left(\mathit{v}\right)& 0& 0\\ 0& -d_{22}\left(\mathit{v}\right)& -d_{23}\left(\mathit{v}\right)\\ 0& -d_{32}\left(\mathit{v}\right)& -d_{33}\left(\mathit{v}\right)\end{array}\right]$$

(12)

$$\mathit{\tau }=\left[\begin{array}{c}{\tau }_X\\ {\tau }_Y\\ {\tau }_N\end{array}\right]$$

(13)

In matrix D, we use the following specific expansion of second order (standard form within our research centre)

$$\{\begin{array}{l}{d}_{11}\left(\mathit{v}\right)=X_{\left|u\right|u}\left|u\right|\hfill \\ d_{22}\left(\mathit{v}\right)=Y_{\left|u\right|v}\left|u\right|+Y_{\left|v\right|v}\left|v\right|+Y_{\left|r\right|v}\left|r\right|\hfill \\ \begin{array}{l}{d}_{23}\left(\mathit{v}\right)=Y_{\left|u\right|r}\left|u\right|+Y_{\left|v\right|r}\left|v\right|+Y_{\left|r\right|r}\left|r\right|\hfill \\ d_{32}\left(\mathit{v}\right)=N_{\left|u\right|v}\left|u\right|+N_{\left|v\right|v}\left|v\right|+N_{\left|r\right|v}\left|r\right|\hfill \\ d_{33}\left(\mathit{v}\right)=N_{\left|u\right|r}\left|u\right|+N_{\left|v\right|r}\left|v\right|+N_{\left|r\right|r}\left|r\right|\hfill \end{array}\hfill \end{array}$$

(14)

For our purposes, the component τ_X in (13) is related only to the propulsive force of the propeller, proportional roughly to its revolutions squared, τ_X = (1 − t_P)·X_nn·n², where t_P—thrust deduction, while the other two components represent the rudder lateral force and moment, τ_Y = Y_uuδ·u²δ and τ_N = N_uuδ·u²δ, respectively (δ < 0 for starboard rudder, in (rad)).

For the purpose of the present research we have adapted the nonlinear model of the slightly unstable large tanker Esso Bernicia,

Table 1. Ship data (ESSO 193 kDWT tanker).

Parameter	Value
Length overall (LOA), m	323.7
Length between perpendiculars (L, LBP), m	304.8
Beam (B), m	47.17
Draught (summer) (T), m	18.46
L/B	6.46
B/T	2.56
Depth (H), m	23.7
Displacement (mass), t	225,500
Block coefficient (cB)	0.83
M/E power, kW	22,400
Propeller rpm, rpm	80
Service speed, kn	16

, that is well published [23] and validated in full scale against widely available and very extensive sea trials. An unstable ship is more challenging for the design and assessment of the autopilot. The selected ship is one of few considered as benchmarks in the area of ship manoeuvring and steering simulation. The original publication [23], however, does not show the exact equations for the structure of forces comprising the listed hydrodynamic derivatives, directly referring to [24]. However, the latter source is too general and has caused a lot of interpretation confusions and, thus, application or implementation problems. A better documentation of the model can be found in [22] that reproduces the original simulation results of [23], but only after correcting the power of the derivatives in the rudder flow speed (‘c’) (also erroneous in [24]) from two to one, which should now look even more natural.

After ensuring that the full model of [22,23] is properly working, and limiting it to the deep water condition, of primary interest in the present preliminary study, we attempted the model simplification to reach the above stated structure with a proper level of adequacy. Our simple standard structure allows to better understand and control the manoeuvring model, as well as to assess the steering simulation results. Among others, the side effect of the propeller, leading to asymmetry in the port and starboard turning circle, was omitted (also of rather negligible effect in the turning circle and zig-zag tests according to the original parameters). In the surge equation, the rudder force and hull resistance coupling with sway velocity also had to be rejected. The other actions are commented on in the annotations (see the footer) to

Table 2. Model parameters and hydrodynamic derivatives (full scale, deep water case).

Name	Value (–)	Scaling Multiplier (in Terms of L, see Table 1, to Be Applied before Entering Values into Equations)
$m$	1 (unity)	-
$J_z$	0.05 1	L2
$X_{\dot{u}}$	−0.05	-
$X_{\left|u\right|u}$	−0.0377	L−1
$X_{nn}$	1.984 × 10−5 2	L
$t_P$	0.22	-
$Y_{\dot{v}}$	−1.02	-
$Y_{\left|u\right|v}$	−1.205	L−1
$Y_{\left|v\right|v}$	−2.4	L−1
$Y_{\left|u\right|r}$	0.298 3	-
$Y_{uu\delta }$	0.2704 4	L−1
$N_{\dot{r}}$	−0.0732	L2
$N_{\left|u\right|v}$	0.519 5	-
$N_{\left|u\right|r}$	−0.231	L
$N_{\left|v\right|r}$	−0.300	L
$N_{uu\delta }$	−0.1274 4	-

¹ Together with $N_{\dot{r}}$ equal to original inertia moment coefficient (with added inertia), i.e., =0.1232. ² Reduced against the original (due to rejecting u-related factors in propeller thrust) by ca. 44% to reach the equivalent effect in balancing the hull resistance at service speed. ³ Considered as viscous contribution; when combined (as negative) with $\left(1-X_{\dot{u}}\right)$, it produces the original value 0.752 for the ‘effective’ centripetal term $\left(1-Y_{\left|u\right|r}\right)$ in [23]. ⁴ Increased by 30% against the original (as inherently related to the mentioned rudder flow average speed ‘c’) to partly account for: (1) the propeller slipstream (the effect of rpm), (2) the omitted original effect of coupling with the drift angle. The displayed value proves somewhat optimal to simulate turning circle and zig-zag as close as possible to those of original model. Both the original rudder force and moment derivatives were here equally affected. ⁵ Extracted from the original −0.451 and to be considered viscous contribution, since the Munk moment derivative $\left(Y_{\dot{v}}-X_{\dot{u}}\right)$ = −0.970 is treated independently and added to our tabulated value.

, which also gathers all of the parameters of our simplified nonlinear model. Instead of a total of 41 derivatives [23], we arrived at 14 derivatives.

The ‘bis’ superscript, originally used in the symbols of derivatives in [23], is omitted in our presentation, but the necessary scaling/multiplication per the last column of Table 2 should be exercised. The damping derivatives in (14), if not listed in Table 2, equal zero.

We use the unity mass, that means ‘1′ has to be substituted for m in (10) and (11). Since the inertia moment and hydrodynamic forces have already been non-dimensionalised (divided), among others, by the ship’s displacement (mass), providing the derivatives of [23], we can operate our algorithm (by default, assigned to the absolute quantities) almost directly with the non-dimensional parameters. In this case, speaking in general, there is no need to divide force by mass in order to obtain acceleration (for change of velocity) because the non-dimensional forces (and moments), and their components, as combined with the ship’s velocities and some linear dimensions, already represent appropriate linear accelerations (and angular accelerations).

Selected standard manoeuvres of our ship model, to prove its performance, are presented in Appendix A.

3.2. Rudder Dynamics

Ignored in the controller design process, the steering gear model will be part of the overall ship model to ensure that the best possible representation of the real ship dynamics is obtained in the simulations (according to [22])

$$\{\begin{array}{l}\left|\delta \right|\le {\delta }_{max} \hfill \\ \left|\dot{\delta }\right|\le {\dot{\delta }}_{max}\hfill \\ \mathrm{if} \left({\delta }_{ord}-\delta \right)/T_{SG}\ge {\dot{\delta }}_{max} \Rightarrow \dot{\delta }={\dot{\delta }}_{max} \nRightarrow \dot{\delta }=\left({\delta }_{ord}-\delta \right)/T_{SG}\hfill \end{array}$$

(15)

where ${\delta }_{max}=35°$; ${\dot{\delta }}_{max}=2.33°/\mathrm{s}$; $T_{SG}=1 \mathrm{s}$.

3.3. Linearised Model

Assuming that velocity u is not changing considerably during rudder manoeuvres, that is valid, to say, for mild manoeuvres, i.e., with low to medium rudder angles (depending upon inherent nonlinearities in the dynamics of a particular ship) we can linearise the dynamics. Within this condition, the low rudder angles should also be followed by small hull drift angles and dimensionless yaw velocities (corresponding to track curvatures), which are the main factors in hull hydrodynamic forces.

Under the above restraints, (u ≈ const, du/dt = 0), the longitudinal ship dynamic equation now disappears. The other two equations in sway and yaw can be classically and straightforwardly written, with pure derivatives of v and r kept on the left-hand side (and the right-hand side components thus divided by relevant mass-related items), as a set of coupled linear equations of first order (or second order vs. heading ψ, where r = dψ/dt) as follows

$$\{\begin{array}{c}\dot{v}=a_{11}\left(u\right)·v+a_{12}\left(u\right)·r+b_1\left(u\right)·\delta \\ \dot{r}=a_{21}\left(u\right)·v+a_{22}\left(u\right)·r+b_2\left(u\right)·\delta \end{array} \mathrm{or} \left[\begin{array}{c}\dot{v}\\ \dot{r}\end{array}\right]={\mathit{A}}_{\mathbf{1}}\left(u\right)\left[\begin{array}{c}v\\ r\end{array}\right]+{\mathit{B}}_{\mathbf{1}}\left(u\right)\delta$$

(16)

where

$${\mathit{A}}_{\mathbf{1}}\left(u\right)=\left[\begin{array}{cc}{Y}_{\left|u\right|v}·L^{-1}/\left(1-Y_{\dot{v}}\right)& \left(Y_{\left|u\right|r}-\left(1-X_{\dot{u}}\right)\right)/\left(1-Y_{\dot{v}}\right)\\ \left(N_{\left|u\right|v}+\left(Y_{\dot{v}}-X_{\dot{u}}\right)\right)·L^{-2}/\left(J_z-N_{\dot{r}}\right)& N_{\left|u\right|r}·L^{-1}/\left(J_z-N_{\dot{r}}\right)\end{array}\right]·u$$

(17)

$${\mathit{B}}_{\mathbf{1}}\left(u\right)=\left[\begin{array}{c}{Y}_{uu\delta }·L^{-1}/\left(1-Y_{\dot{v}}\right)\\ N_{uu\delta }·L^{-2}/\left(J_z-N_{\dot{r}}\right)\end{array}\right]·u^2$$

(18)

In the above, with regard to the full equations of sway and yaw of the previous Section 3.1, just two items have been withdrawn, namely $Y_{\left|u\right|v}$ and $N_{\left|v\right|r}$, as producing null partial derivatives versus v and r around the point (u ≠ 0, v = r = 0). One should be aware that all the parameters (mass-related and hydrodynamic derivatives) in Equations (17) and (18) assume the pure numerical values per the second column of Table 2. The mentioned scaling (multiplication by proper power of L) is not required as this has already been incorporated in (17) and (18).

It is remarkable that the reduced, linear ship model (16)–(18) is still unstable on course, although the linearisation (dependent upon a method) often ‘stabilises’ the ship, especially if we identify the linear model from the zig-zag test.

3.4. Rough Estimate of Rudder Moment Derivative b₂ for Arbitrary Ship

Since when designing the autopilot algorithm we need at least a rough estimate of the parameter b₂ in Equation (16) (see also Section 5 and Section 6)), in this section we derive an approximate formula for b₂, valid for a wide class of merchant ships.

The yaw equation taken from (9), while limiting to the rudder moment (denoted by τ_N) only, i.e., after disregarding matrices C and D, reads

$$\left(J_z-N_{\dot{r}}\right)\dot{r}={\tau }_N$$

(19)

and represents the (additive) contribution of the rudder to the total yaw response. As such, this equation itself also fully describes the ship’s initial motion from rest when hull hydrodynamic forces are not yet developed.

Having adopted the usual, well-known (in ship hydrodynamics) structure of all terms in (19), one obtains the ship’s inertia moment $J_z$, added inertia $N_{\dot{r}}$ and rudder moment ${\tau }_N$, accordingly, as follows

$$J_z=m{L}^2{\left(r_z^{\prime }\right)}^2$$

(20)

$$m=\rho c_{B}LB{T}_{ }$$

(21)

$$N_{\dot{r}}=-k_{66}{J}_z$$

(22)

$${\tau }_N=x_R\left(1+a_H\right)0.5\rho A_R{v}_R^2{c}_{yR}$$

(23)

$$A_R=A_R^{\prime }LT$$

(24)

where: m—displacement (mass); $\rho$—water density; $r_z^{\prime }$—dimensionless radius of gyration; $k_{66}$—dimensionless coefficient of added inertia; $x_R$—rudder effective abscissa (herein geometrically assumed, = −L/2); $a_H$—dimensionless hull-rudder interaction coefficient (augmentation of effective rudder lateral force due to hull presence); $A_R$ and $A_R^{\prime }$—rudder area and its dimensionless coefficient (as ratio of LT), respectively; $v_R$—rudder flow velocity, i.e., reference velocity in the complicated flow in behind propeller/hull conditions for defining the rudder (dimensionless) force coefficients; $c_{yR}$—rudder lateral force coefficient, herein equal to the lift coefficient $c_L$, due to the assumption of the ship’s centreplane straight-line movement).

In the latter, no rudder local drift exists and the rudder angle of attack (‘α’) means rudder angle δ.

Furthermore, we can write down

$$v_R=u\left(1-w\right)\sqrt{1+c_{Th}}$$

(25)

$$c_{Th}=\frac{8}{\pi }\frac{k_T}{J^2}$$

(26)

$$c_{yR}=c_L\left(\alpha ,c_{Th}\right)\approx \delta {\frac{\partial c_L}{\partial \alpha }|}_{c_{Th}=0}$$

(27)

where: u—ship’s longitudinal speed (as input); displacement (mass); w—wake fraction (dimensionless); $c_{Th}$—propeller thrust loading coefficient, as based on standard propeller thrust coefficient (‘k_T’) and advance ratio (‘J’) and which is constant for the steady-state motion with different propeller rpm.

Herein, the propeller jet velocity is assumed for the reference velocity. However, due to the availability of data and roughness of the approach, we take the free-stream (uniform flow) hydrodynamic characteristics of the rudder in terms of the c_L data. This corresponds to null $c_{Th}$. Note that for other, non-zero values of $c_{Th}$, the ordinates of c_L are always lower.

A straightforward transformation of all the Equations (19)–(27) produces the following simple formula

$$\dot{r}=b_2\left(u\right)·\delta$$

(28)

with

$$b_2\left(u\right)=-\frac{0.25\left(1+a_H\right)A_R^{\prime }{u}^2{\left(1-w\right)}^2\left(1+c_{Th}\right)\frac{\partial c_L}{\partial \alpha }}{{\left(r_z^{\prime }\right)}^2{c}_{B}LB\left(1+k_{66}\right)}$$

(29)

In (29), some items (especially those of hydrodynamic nature) may not be known accurately at the lowest effort. Therefore, the authors recommend using the following data of

Table 3. Recommended (rough) values of dimensionless parameters in b₂ formula.

Parameter	Value
$a_H$	0.5
$A_R^{\prime }$	0.02
$w$	0.4
$c_{Th}$	2
$\frac{\partial c_L}{\partial \alpha }$	2.3
$r_z^{\prime }$	0.24
$k_{66}$	1

for any ship type (among seagoing merchant ships) and size. These values, together with three geometrical parameters of our tanker in Table 1, namely c_B, L and B, return b₂ for the speed u = 8.23 m/s equal to −9.18 × 10⁻⁴, as compared to −7.54 × 10⁻⁴ that can be obtained exactly from the parameters of the above ship nonlinear model by means of (16) and (18).

For additional guidance and assessment of Table 3, particularly while modelling a class of ships or examining the limits of b₂, the reader may refer to the wider literature, for example to:

• Reference [25] for $a_H$, where even higher values are permitted. Such a significance of $a_H$ has also been proved in the authors’ experience gained through various identification attempts of mathematical ship models. The empirical meaning of this parameter also tends to correct the theoretically assumed $\frac{\partial c_L}{\partial \alpha }$, see below, since they both stand beside each other in (23) as a product, taking into account (27);
• Reference [26] for $A_R^{\prime }$. However, according to the authors’ database, the provided magnitude 0.02 can be somewhat considered as the upper limit for many 1-screw seagoing merchant ships (e.g., tankers, bulk carriers and containerships), where the value 0.015 can often be found in existing ships;
• Reference [27] for $\frac{\partial c_L}{\partial \alpha }$. The listed value 2.3 comes from the upper boundary of c_L (=1.2), which for the usual rudder aspect (height-to-chord length) ratio of ca. 1.5–1.6 occurs at rudder angle 30° [27];
• Reference [28] for $r_z^{\prime }$, which reports a range of values 0.05 to 0.062 for ${\left(r_z^{\prime }\right)}^2$. However, this can be easily determined from the actual ship’s longitudinal weight distribution. With such a method, the authors confirmed the data of [28] to be partly correct. Additionally, note that the potential error in the gyration radius is often repaired by the empirical assumption (or full-scale identification) of $k_{66}$. The latter is hardly determinable for typical ship hulls, although there are a lot of exact solutions for simple bodies or various regression expressions. In the deep water condition, both the ship’s moment of inertia ($J_z$) and her added inertia ($N_{\dot{r}}$) roughly seem to be of the same order and, hence, $k_{66}=1$ in Table 3.

The wake fraction and the propeller hydrodynamic characteristics, as well as the propeller selection methods (in terms of rpm, pitch, diameter), are well established in the literature. The practical range of w is approximately 0.3 to 0.45 and depends on whether a ship is slender (e.g., container carrier) or full (e.g., tanker, bulk carrier). Although being partly independent, the propeller thrust loading coefficient c_Th varies accordingly between about 1.5 and 3.

4. Zero Dynamics and the Final Controller Design Model

Complex control systems can be typically reduced to the first or second order. Since the ADRC’s order is selected taking into consideration the relative degree [29] of the plant, lumping many non-trackable terms under total disturbance facilitates the simplification.

Applying the input–output (I/O) linearization method, the model is transformed (cf. (16)) as follows

$$\{\begin{array}{l}\dot{v}=a_{11}(u)v+a_{12}(u)r+b_1(u)\delta \\ \dot{r}=a_{21}(u)v+a_{22}(u)r+b_2(u)\delta \\ \dot{\psi }=r\end{array}$$

(30)

with the output

$$y=\psi$$

(31)

The resultant form is suitable for direct application of the ADRC technique in accordance with the theory outlined in Section 2. However, in order to avoid the here typical formalism of Lie derivatives, the output (31) is differentiated twice with respect to time, to obtain the equations

$$\begin{array}{l}{\dot{y}}_1=y_2\\ {\dot{y}}_2=f(x)+g(x)\delta \end{array}$$

(32)

where $y=y_1=\psi$, $\dot{y}=y_2=r$.

$$\begin{array}{l}f(x)=a_{21}(u)v+a_{22}(u)r\\ g(x)=b_2\end{array}$$

(33)

The partial transformation of the coordinates from x to y is defined by the equations

$$y_1=\psi ; y_2=r$$

(34)

The transformation generates a new system (32) of the second order, representing the external dynamics. Since the original system (30), (31) is of the third order, the missing part of the dynamics must be complemented with the internal (zero) dynamics (of the first order).

An analysis of the system zero dynamics can be performed through examination of the original system (30), (31) under the assumption that its output equals zero $y(t)\equiv 0$. Thereby, we obtain the restricted motion of the system [29], confined to the set

$$Z^{*}=\left\{x:h(x)=L_{f}h(x)=0\right\}=\left\{x:r=0\right\}$$

(35)

where $y=h(x)$

The motion of the original system (30), (31) on $Z^{*}$, with the input $\delta =\alpha (x)=-f(x)/g(x)$ (compare (33)), represents the asymptotically stable zero dynamics

$$\dot{v}=a_{0}v\hspace{1em}\mathrm{where}\hspace{1em}{a}_0=(a_{11}-a_{21}\frac{b_1}{b_2})<0$$

(36)

Thus, for the purpose of the controller design, the model (30), (31) can be represented by means of the system of second order (32), written as in (1)

$$\begin{array}{l}\{\begin{array}{l}\dot{\psi }=r\\ \dot{r}=a_{11}(u)v+a_{22}(u)r+b_2\delta \end{array}\\ y=\psi \end{array}$$

(37)

5. Autopilot Design Using the ADRC Technique

The ship course-keeping controller will be designed using the ADRC technique discussed above. The ADRC concept is based on the assumption that the object model structure as well as its parameters are unavailable, i.e., the parameters and exogenous disturbances that affect the process are not known. The only known data are the roughly estimated rudder efficiency coefficient $b_2$ of (37) (see (29)). Based on the framework provided by the second-order model (vs. heading) of (37), we can simply apply the theory given in Section 2, in particular Formula (5), defining the controller, and (7) specifying the model for observer design. To maintain consistency with Section 2, denotation of the rudder efficiency coefficient $b_2$ has been changed to $c_0$. The ADRC controller gains, $K_p,K_d$, have been calculated by the pole-placement method (see (6), (38)), and the observer gain matrix $\mathit{L}={\left[\begin{array}{ccc}1.0484& 0.44955& 0.001\end{array}\right]}^T$ has been tuned using the standard LQR procedure.

6. Simulation Results

Based on the mathematical model of the ship (with the steering gear attached) presented in the previous Sections and the general idea of the ADRC system, simulation tests were carried out using the MATLAB/Simulink environment.

The block diagram of the ADRC-based heading angle control system is shown in Figure 1.

The following symbols are used: ψ—heading angle, ε—error signal, α—steering gear command, δ—rudder deflection angle, r—angular speed and d—disturbances.

At the first step, it is necessary to determine the parameters of the ADRC system, i.e., $c_0$, $K_p$ and $K_d$.

The approximate value of $c_0$ is obtained from the theoretical study (see Section 3.4) where $c_0=b_2=$−9.18 × 10⁻⁴ (1/s²). $K_p$ and $K_d$ values are determined via the pole-placement method [30] based on the canonical system (6) as follows

$$K_p=\frac{36}{T_{settle}^2}, \hspace{1em} K_d=\frac{2}{T_{settle}}$$

(38)

where T_settle denotes, to some extent, an arbitrarily assumed settling time of the control system (6).

The second step consists of checking the static and dynamic properties of the ADRC system with parameters calculated as above. The results of this experiment are shown in Figure 2.

In the third stage of the research the robustness of the ADRC system to disturbances was examined.

Two kinds of disturbances were selected during the test: constant disturbance of 25 degrees (in terms of the counteracting rudder deflection), which represents the impact of the wind or sea current, and a colourised white noise corresponding to the impact of waves. Results of the simulation with constant disturbance are shown in Figure 3.

The noteworthy fact is that the maximal heading angle deviation is about 5 degrees. The steady-state value of the heading angle (not entirely shown in this graph) is equal to zero.

The obtained settling time $T_{settle}$ is calculated for the case when the deviation of the course angle is equal to or less than one degree and amounts to about 800 s.

During the last test, a white noise generator with the sampling time and power parameter set as 0.5 s and 0.0001, respectively, was connected to a band-pass filter, described by the following transfer function [10,31]

$$G(s)=\frac{10.8s}{s^2+0.672s+2.6729}$$

(39)

The obtained simulation results are presented in Figure 4 and Figure 5. The first shows variations of disturbance d(t) as coloured-noise, expressed in degrees, whereas in Figure 5 the heading angle vs. time is shown.

As can be seen, the course angle deviation values obtained are practically negligible. However, the robust response to random disturbances requires high performance of the steering gear machinery.

7. Discussion and Conclusions

If we compare the above performance of ADRC to PID with laboratory, well-tuned controller parameters (in our case resolved to: K_p = 3.81, K_d = 3.95, K_i = 0.00547), then the transients, for the step disturbance, are shown in Figure 6.

Differences in the response between ADRC and PID seem to be practically negligible, so that ADRC, in the conditions of our ship simulation model, proved to be similar in performance to PID.

However, it should be strongly emphasized that such a good tuning of the PID controller is not achievable in practice. From this point of view, and the fact of the simple, global tuning of the ADRC controller, there is no doubt that this controller performs much better.

Unfortunately, due to the research and development on ADRC for ship steering being at an early stage, and hence a significant lack of other reference data, especially for a more or less direct comparison and assessment of our ADRC results, we are not able to put them in a much wider context.

The ADRC autopilot tests carried out are very preliminary, but the results are already promising. The proposed system with the ADRC algorithm is characterized not only by robustness against the influence of constant and random disturbances but the first priority is the ease of its tuning, which is relatively simple.

We just need to know the parameter $b_2$ $(c_0)$ and the gains $K_p,K_d$. However, for the estimation of $b_2$, we proposed the formula (29), valid for a wide class of merchant vessels, while the gains $K_p,K_d$ can be easily determined by the general canonical model (6). This comes down to the choice of settling time $T_{settle}$ (38) that can be rather arbitrarily performed under laboratory conditions. As shown in Figure 2, due to the robustness of this parameter, its nominal selection (the settling time 150 s) should ensure sufficiently good control performance for a real ship.

There is no doubt that the presented system requires further research. Important issues are, for example, its functioning with the introduced dead-band in order to reduce the intensity of the steering gear operation, or the issue of assessing its performance using integral quality indices (cf. [10]).

Nevertheless, the final evaluation of the system would require tests on a real ship, which, due to the fact that the simulation tests were carried out on a verified, full ship model, provides a chance that the implementation will work correctly.