Path following of under-actuated ships based on backstepping and predictive control method

Abstract

To deal with the sideslip angle caused by the current disturbances or transverse motion for path following of under-actuated ships, a nonlinear observer established by an exponential function is introduced in the backstepping approach which converts the path following into heading control. Then, the model predictive control (MPC) method is used as a heading controller, addressing the rudder optimization. A linear extended state observer technology was exploited to estimate yaw rate, external disturbances, and internal uncertainties, which could avoid measuring the high-order state used in the MPC controller and promote the accuracy of the MPC internal model. Moreover, an inverse tangent function is applied to develop a new method for switching the reference heading angle to reduce rudder amplitude when the ship is choosing the next waypoint. Finally, the validity and reliability of the design method were verified through comparative computer simulation experiments.

Introduction

In recent years, the path-following control of underactuated ships has become an important research topic in nonholonomic systems. ¹ Fossen applied the Line of Sight (LOS) guidance to a straight linear path following, where the algorithm was converted to heading control. ² Then, the integral LOS was presented in Lekkas and Fossen ³ and Chen et al. ⁴ to offset the sideslip angle due to current disturbances or transverse motion. To avoid the sideslip angle, the adaptive method was proposed in Huang and Fan. ⁵ However, the above algorithms addressed the small or slow sideslip angle only, the linear extended state observer (LESO) was presented in Liu et al. ⁶ and Wang et al. ⁷ for estimating the sideslip angle. However, external disturbances, model uncertainties and unknown parameters should also be considered. An iterative sliding mode control with reinforced learning was proposed in Shen and Dai ⁸ and Shen et al. ⁹ where the environmental disturbances and internal uncertainties can be ignored. Then, the neuron adaptive iterative sliding mode was presented in Shen and Jing ¹⁰ to address the uncertainty of parameters and unknown disturbances. In Wang et al., ¹¹ a multiobjective genetic algorithm was used to address the restrictions caused by the model perturbation and external interferences. Robust and adaptive control was applied in Zwierzewicz ¹² and Padideh et al. ¹³ to offset the modeling inaccuracies and external interferences. The adaptive technique was also applied in Vo et al. ¹⁴ and Hu et al. ¹⁵ to address the nonlinear hydrodynamic term or external disturbances. Then, radial basis function (RBF) neural networks were applied in Shen et al. ¹⁶ and Wan et al. ¹⁷ to offset model uncertainties. In Zeng et al., ¹⁸ a simplified RBF neural network was proposed. In Yu et al. ¹⁹ and Zhang et al.^,20 a minimum learning parameter (MLP) neural network approach was used to decrease the calculation burden. The MLP algorithm was also applied in Mu et al. ²¹ and Shen and Wang, ²² which could tackle the uncertainty of dynamics and external interferences. Most RBF or MLP weights or adaptive laws were updated by combining with the tracking errors, but when the tracking error changes rapidly, the accuracy of the tracking is reduced. In Ding et al., ²³ an improved observer was proposed to estimate unknown interferences. A nonlinear observer was proposed in Zheng and Feroskhan ²⁴ where the estimation of interferences is more accurate. In addition, an active disturbance rejection control (ADRC) method, presented by Han based on proportional integral derivative control, Han ²⁵ has been successfully used in many control areas where the system has unknown parts. In Huang and Fan^,5 the ADRC technology was proposed with LOS law for the path following, and the key of ADRC, namely extended state observer (ESO), was used to estimate the disturbances. Then, the ADRC was extended by Gao ²⁶ to linear ADRC whose key is LESO, due to the difficulty of adjusting the ESO parameters. In Li et al., ²⁷ the LESO was applied to solve the uncertainty of dynamics and external interferences. The LESO was also used in Li et al., ²⁸ it promotes robustness to wave to address the control input constraints, the auxiliary system was used in Fan et al. ²⁹ and Wang et al. ³⁰ In Liu et al., ³¹ the constraints could be handled by the bounded nonlinear function in the iterative sliding mode controller. In Nie and Lin, ³² a novel simple attitude-tracking compensator was introduced to solve the input saturation. There is a standard optimization method for complex constrained systems, ³³ namely model predictive control (MPC), which can directly incorporate dynamic constraints to address the control input constraints easily. ³⁴ In Yao et al., ³⁵ the MPC method was used to obtain an optimal expected angular velocity satisfying the constraint. The linear MPC was used in Hu et al., ³⁶ where input range and incremental constraints were considered. However, linear MPC may introduce inaccuracy when the external disturbances are large. In Liu et al., ³⁷ the MPC algorithm with improved least squares support vector machines were presented for promoting the controller robustness under disturbances or uncertainties. Both linear MPC and nonlinear MPC were proposed in Zheng et al. ³⁸ to address input constraints. However, the inaccuracy of the internal model and the problem of measuring high-order state always increase the difficulties of using the MPC algorithm.

Most of the existing observers are unable to resolve external disturbances, model uncertainties and unknown parameters. Considering abovementioned research, this paper combines the backstepping method with the MPC algorithm for path following. Model predictive control has the advantages of good control and robustness, which can effectively overcome process uncertainty, nonlinearity and parallelism, and can easily handle various constraints in the process controlled and manipulated variables. Motivated by Wang et al., ⁷ the nonlinear observer is designed, via the exponential function, for estimating the sideslip angle. The MPC is exploited as the heading controller, addressing the rudder constraints and optimization. Furthermore, the LESO is used to estimate the yaw rate, external interferences as well as model uncertainties, which could avoid measuring high order state and promote the accuracy of the internal model. Moreover, an inverse tangent function is applied to construct a switching reference heading angle, such that the rudder amplitude can be reduced when the ship is selecting the next waypoint.

The main contents of the article are as follows:

In part 2, the control objective is presented.

In part 3, a controller for ship path following based on backstepping and MPC method was designed.

In part 4, compared simulation analyses were carried out.

In part 5, the concluding remarks are given.

Control objective

The ship motion parameters are shown in Figure 1.

Figure 1.

The ship motion parameters.

where u denotes surge velocity, v denotes sway velocity, r denotes yaw rate, ψ is heading angle, U  =  (u²  +  v²)^1/2 is ground velocity, and β = arctan(v/u) denotes sideslip angle, which is difficult to calculate because v cannot be measured directly. Then, the ship path-following model is:

$$\{\begin{array}{c}\dot{x}=\sqrt{u^2+v^2}\cos (\psi +\beta )\\ \dot{y}=\sqrt{u^2+v^2}\sin (\psi +\beta )\\ \dot{\psi }=r\\ \dot{r}=-{\displaystyle \frac{1}{T}r-{\displaystyle \frac{a}{T}{r}^3+{\displaystyle \frac{K}{T}\delta +f}}}\end{array}$$

(1)

where K and T denote maneuverability index, a denotes a positive coefficient, δ denotes rudder angle, and f denotes total unknowns. The ship cross-track error is shown in Figure 2.

Figure 2.

The ship cross-track error.

where y_d denotes a reference path desired by the waypoint, y_e denotes cross-track error, θ denotes straight line angle, and R denotes turning range radius. Then, the control objective is to design δ that is ship following the reference path, namely satisfy y_e  =  0.

Path-following controller

In the path-following control system, the article uses the backstepping method to control the reference heading angle and the MPC method to control ship's heading. Furthermore, the LESO is used to estimate the yaw rate, external interferences as well as model uncertainties, which could avoid measuring high order state and promote the accuracy of the internal model. Then, the control design is shown in Figure 3.

Figure 3.

The control design diagram.

The backstepping algorithm for reference heading

In this section, the backstepping algorithm was used to design reference heading angle. It allows complex nonlinear systems to be decomposed into subsystems of no more than the order of the system, and then partial Lyapunov functions and intermediate virtual control quantities are designed for each subsystem until the entire control law is completed. The backstepping method is applied to produce the initial reference heading angle ψ_d1 and a nonlinear observer was designed to estimate sideslip angle β. This design philosophy ensures a stable convergence of the system and ensures stable tracking of the heading, reducing system chatter. According to Figure 2, we have

$$y_e={\displaystyle \frac{Ax+By+C}{\sqrt{A^2+B^2}}={\displaystyle \frac{x\tan \theta -y+C}{\sqrt{{\tan }^2\theta +1}}}}$$

(2)

where A, B, and C denote different factors. Then,

$$\begin{array}{rl}{\dot{y}}_e=& (\dot{x}\tan \theta -\dot{y})/\sqrt{{\tan }^2\theta +1}\\ =& \tan \theta \sqrt{u^2+v^2}(cos(\psi +\beta )-sin(\psi +\beta ))/\sqrt{{\tan }^2\theta +1}\\ =& \sqrt{u^2+v^2}sin(\theta -\psi -\beta )\end{array}$$

(3)

Hence, the initial reference heading angle ψ_d1 can be designed

$${\psi }_{\mathrm{d}1}=k_1\tanh (k_2{y}_e)+\theta -\hat{\beta }$$

(4)

where k₁ and k₂ are the positive parameters, $\hat{\beta }$ denote the estimation of β, and considering estimation error $\tilde{\beta }=\hat{\beta }-\beta$ is small, namely $|\tilde{\beta }|<\mathrm{\Delta }\beta$ , Δβ is a small positive constant.

Proof:

Choosing the Lyapunov function V₁  =  (1/2)y_e² and then,

$$\begin{array}{rl}{\dot{V}}_1& =y_e{\dot{y}}_e\\ =& -\sqrt{u^2+v^2}{y}_e\sin (k_{1}tanh(k_2{y}_e)+\beta -\hat{\beta })\\ =& -U{y}_e\sin (k_{1}tanh(k_2{y}_e)-\tilde{\beta })\\ \approx & -U{y}_e\sin (k_{1}tanh(k_2{y}_e))\le 0\end{array}$$

(5)

Thus, the designed ψ_d1 could guarantee y_e is bounded and converges to 0. The next waypoints P _{k  + 2} should be selected when ship position P_{(x, y)} satisfies

$$(x-x_{k+1}{)}^2+(y-y_{k+1}{)}^2\le R^2$$

(6)

Then, to estimate β and motivated by a linear observer in Wang et al., ⁷ a nonlinear sideslip angle observer is designed. Firstly, the kinematics is rewritten as

$$\begin{array}{rl}\dot{y}=& u\sin \psi +v\cos \psi \\ =& u\sin \psi +\omega \end{array}$$

(7)

where $\omega =v\cos \psi$ . The state observer will be built to estimate ω.

$$\{\begin{array}{c}\dot{\hat{y}}=u\sin \psi +\hat{\omega }-k_{y1}{h}_{\hat{y}-y}\\ \dot{\hat{\omega }}=-k_{y2}{h}_{\hat{y}-y}\end{array}$$

(8)

where k_y1 and k_y2 are the positive parameters, ŷ denotes estimation of y, $\hat{\omega }$ denote estimation of ω, h is a bounded nonlinear function based on the exponential, and built as

$$h_{\epsilon }=\{\begin{array}{c}1-e^{-k_h\epsilon },\epsilon \ge 0\\ e^{k_h\epsilon }-1,\epsilon <0\end{array}$$

(9)

where k_h is a positive parameter. We can see that, h∈(−1, 1) is an increasing odd function.

Proof:

We select the positive coefficients η and choose the Lyapunov function $V_1=({(\hat{y}-y)}^2+\eta {(\hat{\omega }-\omega )}^2)/2$ . Differentiating V₂ with respect to time, we have

$$\{\begin{array}{c}{\dot{V}}_2=(\hat{y}-y)(\dot{\hat{y}}-\dot{y})+\eta (\hat{\omega }-\omega )(\dot{\hat{\omega }}-\dot{\omega })\\ =-k_{y1}(\hat{y}-y)h_{\hat{y}-y}+z_1{z}_2+z_3\\ z_1=\hat{\omega }-\omega \\ z_2=(\hat{y}-y)-\eta k_{y2}{h}_{\hat{y}-y}\\ z_3=-\eta \dot{\omega }(\hat{\omega }-\omega )\end{array}$$

(10)

Considering the characteristic of h and choosing the suitable parameter k_y2 < 1/η, z₂ and (ŷ−y) have the same sign. If (ŷ−y) > 0 and z₂ > 0, then $\dot{\hat{\omega }}<0$ , and $\hat{\omega }$ will decrease to $\hat{\omega }<\omega$ or z₁ < 0, hence z₁z₂ < 0. Meanwhile, if (ŷ−y) < 0 and z₂ < 0, then $\dot{\hat{\omega }}>0$ , and $\hat{\omega }$ will increase to $\hat{\omega }>\omega$ or z₁ > 0, we still have z₁z₂ < 0. The assumption, z₃ is small compared with the other parts, may be satisfied because ω is bounded and time-varying slowly. So, there is a small positive constant σ, leading to z₁z₂  +  z₃≤ σ. Thus, we have

$$\begin{array}{rl}{\dot{V}}_2=& -k_{y1}(\hat{y}-y)h_{\hat{y}-y}+z_1{z}_2+z_3\\ \le & -k_{y1}(\hat{y}-y)h_{\hat{y}-y}+\sigma \le 0\end{array}$$

(11)

Then, the estimated sideslip angle $\hat{\beta }$ can be computed by $\hat{\beta }=\arctan (\hat{\omega }/(u\cos \psi ))$ .

Model predictive control method for heading controller

For the MPC's path-following problem, a suitable predictive control method not only helps to ensure the stability of the control system but also reduces the amount of computation and increases the calculation speed. An improved MPC method is proposed to design a heading controller. Considering the large control input may be designed when the reference value changes suddenly and greatly, a new switching method is proposed there. We introduce an inverse tangent function to design the switching reference heading angle

$${\psi }_{\mathrm{d}}=\{\begin{array}{c}{\psi }_{\mathrm{d}1},|{\psi }_{\mathrm{d}1}-\psi |\le {\psi }_{max}\\ \psi +k_{\psi }\arctan ({\psi }_{\mathrm{d}1}-\psi ),|{\psi }_{\mathrm{d}1}-\psi |>{\psi }_{max}\end{array}$$

(12)

where ψ_d is new reference heading angle, ψ_max is the switching range, and k_ψ is a positive parameter. The MPC is used to design the bow controller so that the bow follows the reference bow. Using the Norrbin model as a prediction model, the system state is discretized and predicted based on the Euler iterative method as follows. We use formula (13) to discrete and predict the system states.

$$\{\begin{array}{c}\overline{\psi }(k+1)=\psi (k)+T_c{\dot{f}}_{\psi }(k)\\ \overline{\psi }(k+1)=\hat{r}(k)+T_c{\dot{f}}_r(k)\end{array}$$

(13)

where ȓ(k) denotes the estimation of r(k), and T_c is the predictive sampling time:

$$\{\begin{array}{cc}{\dot{f}}_{\psi }(k)& =\hat{r}(k)\\ {\dot{f}}_r(k)& =-{\displaystyle \frac{1}{T}\hat{r}(k)-{\displaystyle \frac{a}{T}(\hat{r}(k){)}^3+{\displaystyle \frac{K}{T}\delta +\hat{f}}}}\end{array}$$

(14)

where $\hat{f}$ denote the estimation of total unknowns. The next heading angle is:

$$\overline{\psi }(k+N_p)=\psi (k)+\sum _{j=0}^{N_p-1}(T_c{\dot{f}}_{\psi }(k+j))$$

(15)

The error ē_ψ(k  +  j), j  =  1, 2, …, N_P denotes propeller moment. The value of N_p and ${\psi }_{\mathrm{d}}$ are related as in equation (16):

$${\overline{e}}_{\psi }(k+N_p)=\overline{\psi }(k+N_p)-{\psi }_{\mathrm{d}}(k+N_p)$$

(16)

A function with a rudder angle is created from the prediction error and the current error:

$$\underset{\genfrac{}{}{0ex}{}{{\delta }_{min}\le \delta \le {\delta }_{max}}{\mathrm{\Delta }{\delta }_{min}\le \mathrm{\Delta }\delta \le \mathrm{\Delta }{\delta }_{max}}}{min}{f}^{*}({\overline{e}}_{\psi },\delta )=\sum _{j=0}^{N_p}Q\overline{e}{(k+j)}^2+\sum _{i=1}^{Nc}G\delta {(k+i)}^2$$

(17)

where Q and G are weights used to balance the relationship between control performance and rudder angle performance. Q is used to control the path deviation and G controls the rudder angle. The greater the path deviation, the greater the rudder angle required by the ship. Therefore, the larger the Q value, the larger the G value required. Δδ denotes rudder rate, N_c denotes control horizon, and N_c < N_p. By setting the constraints δ_min≤ δ ≤ δ_max and Δδ_min≤ Δδ ≤ Δδ_max,. The stability of MPC has been proved in Marjanovic and Lennox. ³⁹

Linear extended state observer algorithm for evaluating total unknowns and yaw rate

The traditional LESO technique has more solutions for external disturbances and model uncertainties, but less for control input constraints and system optimization. In this paper, a path-following predictive control algorithm combined with LESO is proposed to deal with the rudder angle constraint and optimization problem. The traditional MPC with a single structure is improved by using the LESO technique to approximate the uncertainties and external disturbances and to compensate for the prediction model with feedback. In addition, the longitudinal and lateral velocities and the yaw rate are estimated using a nonlinear observer and LESO respectively, avoiding the problem that the required system state values for the controller are not easily measured.

In path-following control, it is vital to deal with external disturbances and to calculate the yaw rate. In this section, the LESO is used to evaluate total unknowns f and yaw rate r. Considering equation (1), the improved LESO can be written as follows:

$$\{\begin{array}{c}\dot{\hat{\psi }}=\hat{r}-l_1(\hat{\psi }-\psi )\\ \dot{\hat{r}}=-{\displaystyle \frac{1}{T}\hat{r}-{\displaystyle \frac{a}{T}{\hat{r}}^3+{\displaystyle \frac{K}{T}\delta +\hat{f}-l_2(\hat{\psi }-\psi )}}}\\ \dot{\hat{f}}=-l_3(\hat{\psi }-\psi )\end{array}$$

(18)

where l₁, l₂, and l₃ denote positive values. $\hat{\psi }$ , ȓ, and $\hat{f}$ denote estimations of ψ, r, and f, respectively.

Proof:

Firstly, the observer estimation errors are obtained

$$\{\begin{array}{c}\dot{\tilde{\psi }}=\tilde{r}-l_1\tilde{\psi }\\ \dot{\tilde{r}}=-{\displaystyle \frac{1}{T}\tilde{r}-{\displaystyle \frac{a}{T}({\hat{r}}^3-r^3)+\tilde{f}-l_2\tilde{\psi }}}\\ \dot{\tilde{f}}=-l_3\tilde{\psi }-\dot{f}\end{array}$$

(19)

where $\tilde{B}=\hat{B}-B$ (B denotes ψ, r and f) are estimation errors. We select the following Lyapunov function:

$$V_3={\displaystyle \frac{1}{2}{\tilde{\psi }}^2+{\displaystyle \frac{1}{2}\chi {\tilde{r}}^2+{\displaystyle \frac{1}{2}\zeta {\tilde{f}}^2}}}$$

(20)

where χ and ζ denote positive values. Then,

$$\{\begin{array}{c}{\dot{V}}_3=\tilde{\psi }\dot{\tilde{\psi }}+\chi \tilde{r}\dot{\tilde{r}}+\zeta \tilde{f}\dot{\tilde{f}}\\ =\tilde{\psi }(\tilde{r}-l_1\tilde{\psi })+\chi \tilde{r}\left(-{\displaystyle \frac{1}{T}\tilde{r}-{\displaystyle \frac{a}{T}({\hat{r}}^3-r^3)+\tilde{f}-l_2\tilde{\psi }}}\right)\\ +\zeta \tilde{f}(-l_3\tilde{\psi }-\dot{f})\\ =-l_1{\tilde{\psi }}^2-{\displaystyle \frac{1}{T}\chi {\tilde{r}}^2-{\displaystyle \frac{a}{T}\chi \tilde{r}({\hat{r}}^3-r^3)}}\\ +(1-\chi l_2)\tilde{\psi }\tilde{r}+\tilde{f}(\chi \tilde{r}-\zeta l_3\tilde{\psi })-\zeta \tilde{f}\dot{f}\\ =-l_1{\tilde{\psi }}^2-{\displaystyle \frac{1}{T}\chi {\tilde{r}}^2-{\displaystyle \frac{a}{T}\chi \tilde{r}({\hat{r}}^3-r^3)+w_1+w_2+w_3}}\\ w_1=(1-\chi l_2)\tilde{\psi }\tilde{r}\\ w_2=\tilde{f}(\chi \tilde{r}-\zeta l_3\tilde{\psi })\\ w_3=-\zeta \tilde{f}\dot{f}\end{array}$$

(21)

Then, if l₂  =  1/χ, then w₁  =  0. According to formula (19), we have $\tilde{r}=\dot{\tilde{\psi }}+l_1\tilde{\psi }$ , next

$$\begin{array}{rl}{w}_2=& \tilde{f}(\chi \tilde{r}-\zeta l_3\tilde{\psi })=\tilde{f}(\chi (\dot{\tilde{\psi }}+l_1\tilde{\psi })-\zeta l_3\tilde{\psi })\\ =& \chi \tilde{f}\dot{\tilde{\psi }}+\tilde{f}\tilde{\psi }(\chi l_1-\zeta l_3)\end{array}$$

(22)

Considering assumption II, so there is a small constant μ which satisfies w₁  +  w₂  +  w₃≤ μ. Then, we have

$$\begin{array}{rl}{\dot{V}}_3=& -l_1{\tilde{\psi }}^2-{\displaystyle \frac{1}{T}\chi {\tilde{r}}^2-{\displaystyle \frac{a}{T}\chi \tilde{r}({\hat{r}}^3-r^3)+w_1+w_2+w_3}}\\ \le & -l_1{\tilde{\psi }}^2-{\displaystyle \frac{1}{T}\chi {\tilde{r}}^2-{\displaystyle \frac{a}{T}\chi (\hat{r}-r)({\hat{r}}^3-r^3)+\mu }}\\ \le & 0\end{array}$$

(23)

Simulation analysis

The problem of achieving a ship's path following is a very complex process. It requires the sophisticated cooperation of several devices such as rudders and propellers. In the simulation experiments, the parameters must not only be designed to match the actual conditions of the ship but also to ensure the stability of the control system.

Simulation model

The ship mathematical model group (MMG) model is applied in the simulation ship motion model finally, expressed as Zhang et al.^40,41:

$$\{\begin{array}{c}\dot{x}=\sqrt{u^2+v^2}\cos (\psi +\beta )\\ \dot{y}=\sqrt{u^2+v^2}\sin (\psi +\beta )\\ \dot{\psi }=r\\ \dot{u}=[(m+m_y)vr+X_H+X_P+X_R+X_W+X_{Wave}\\ +(m_x-m_y)V_c\sin ({\psi }_c-\psi )r]/(m+m_x)\\ \dot{v}=[(m+m_x)ur+Y_H+Y_P+Y_R+Y_W+Y_{Wave}\\ -(m_x-m_y)V_c\cos ({\psi }_c-\psi )r]/(m+m_y)\\ \dot{r}=(N_H+N_P+N_R+N_W+N_{Wave})/(I_{zz}+J_{zz})\end{array}$$

(24)

where m denotes ship mass, m_x and m_y denote additional masses, X_H, Y_H, and N_H denote ship hull forces (moment), X_P, Y_P, and N_P denote propeller forces (moment), X_W, Y_W, and N_W denote wind forces (moment), X_Wave, Y_Wave, and N_Wave denote wave forces (moment). ψ_c and V_c denote current direction and velocity, respectively. I_ZZ denotes moment of inertia of the ship around the vertical axis, J_ZZ denotes additional moment of inertia, and X_R, Y_R, and N_R denote rudder forces (moment) given by

$$\{\begin{array}{c}{X}_R=-(1-t_R)F_N\sin \delta \\ Y_R=-(1+a_H)F_N\cos \delta \\ N_R=-(x_R+a_H{x}_H)F_N\cos \delta \end{array}$$

(25)

where t_R denotes rudder resistance deductions, α_H denotes ratio of lateral force, x_H denotes distance from lateral force center to gravity, F_N denotes rudder positive pressure.

The M.V “YULONG” is used to prove the validity of the proposed method in the simulation. The main parameters of the ship are shown in Table 1.

Table 1.

Parameters of M.V “YU LONG.”

Length L/m	126
Width B/m	20.8
Draft d/m	8.0
Displacement ▽/t	14,625
Block coefficient Cb	0.681
RPM r/min	100
Index K	216
Index T	0.478

The dynamics are given by MMG model (3), and the wind (moment) X_W, Y_W and N_W are given by Wang et al. ⁴² :

$$\{\begin{array}{c}{X}_W=-0.5{\rho }_a{A}_f{U}_R^2{C}_{wx}({\alpha }_R)\\ Y_W=0.5{\rho }_a{A}_s{U}_R^2{C}_{wy}({\alpha }_R)\\ N_W=-0.5{\rho }_a{A}_s{L}_{oa}{U}_R^2{C}_{wn}({\alpha }_R)\end{array}$$

(26)

where ρ_α denotes air density, α_R denotes wind angle of attack, U_R denotes relative wind speed, A_f and A_s denote frontal and lateral projected areas, respectively. L_oα denotes length overall, C_wx(α_R), C_wy(α_R), and C_wn(α_R) denote nondimensional wind coefficients, respectively.

The waves (moments) X_Wave, Y_Wave and N_Wave are given by Qian et al. ⁴³ :

$$\{\begin{array}{c}{X}_{Wave}=0.5\rho L{\alpha }_s^{2}cos\chi C_{Xw}(\lambda )\\ Y_{Wave}=0.5\rho L{\alpha }_s^2\sin \chi C_{Yw}(\lambda )\\ N_{Wave}=0.5\rho L{\alpha }_s^2\sin \chi C_{Nw}(\lambda )\end{array}$$

(27)

where λ denotes wavelength, χ denotes encounter angle, ρ denotes seawater density, α_s denotes wave amplitude, C_Xw(λ), C_Yw(λ), and C_Nw(λ) denote coefficients related to wavelength λ, respectively.

Simulation results

Case 1. To illustrate the ability of avoiding the great rudder amplitude when ship moving at the turning waypoint by introducing the proposed switching reference heading method, we choose three reference heading angles for the compared simulation, namely: (1) ψ_d1 without switching, (2) ψ_ds1 with the switching (k_ψ = 0.7), (3) ψ_ds2 with the switching (k_ψ = 0.5). The others control parameters: k₁ = π/4, k₂  =  0.003, k_y1  =  0.11, k_y2  =  0.009, k_h  =  20, ψ_max  =  1.5°, T_c  =  1.5 s, N_p  =  24, N_c  =  1, l₁  =  1.2, l₂  =  0.48, l₃  =  0.064. The desired waypoint: P₁(0, 0), P₂(4000, 0), P₃(7000, 3000), P₄(8000, 7000), and P₅(9500, 7000). The initial states: u  =  7.0 m/s, v  =  0, r  =  0, ψ = 0, (x₀, y₀)  =  (0, 100 m). The values of external disturbances are shown in Table 2:

Table 2.

The values of external disturbances.

Wind speed m/s	10
Wind direction °	45° + 30°sin(0.02t)
Current speed m/s	1
Current direction°	45° + 10°sin(0.003t)
Wavelength m	83
Wave encounter angle °	ψ + 135°−30°sin(0.02t)
Wave amplitude m	3

The simulation results based on the different reference heading angles are shown in Figures 4 to 8. Figures 4 to 6 represent the path, expected reference heading, actual heading, rudder angle and rate based on three different reference heading angles (without switching ψ_d1, switching ψ_ds1 and ψ_ds2). In Figure 4, we can see that, all y₁, y_s1, and y_s2 can follow the desired path y_d successfully under the environmental disturbances, and y_s2 could follow the next path smoothly when ship is moving at the turning point. Figure 5 shows that, the jumps of ψ_ds1 and ψ_ds2 (based on switching) are less than ψ_d1 (without switching) at the turning point, and the actual heading ψ_s1 and ψ_s2 are more smoother. In Figure 6, δ₁, δ_s1, and δ_s2 are variable due to the time-varying interferences such as wave, wind and current. However, δ_s1 and δ_s2 are more smaller and smoother than δ₁, and δ_s2 could avoid the large rudder angle effectively when ship is selecting the next waypoint. All the rudder rates Δδ₁, Δδ_s1, and Δδ_s2 are less than 6°/s, namely can be constrained within the saturation. Figure 7 illustrates the validity of the presented sideslip observer, and the sideslip angle β is approximated successfully. Figure 8 shows that the yaw rate r and total unknowns f could be estimated accurately by the LESO.

Figure 4.

The path following based on different reference headings.

Figure 5.

The different reference headings and actual headings.

Figure 6.

The rudder angle and rate.

Figure 7.

The sideslip angle.

Figure 8.

The yaw rate and total unknowns.

Case 2. To illustrate the effectiveness of addressing the rudder constraints and optimization via the proposed MPC algorithm, we compare the designed MPC controller with the PD controller proposed in Li et al. ²⁷ and the sliding mode method (SMC: the sliding surface is $s={\dot{\psi }}_e+0.1\tanh (0.12{\psi }_e)$ and $\dot{s}=-0.1s$ ). The desired waypoint: P₁(0, 0), P₂(4000, 0), P₃(7000, 3000), P₄(10,000, 3000), P₅(12,000, 1500), and P₆(13,500, 1500). The initial states: u  =  7.2 m/s, v  =  0, r  =  0, ψ = 0, (x₀, y₀)  =  (0, −100 m). The current set is −45° + 10°sin(0.003t), other parameters are same as the Case 1.

In Figure 9, y_MPC, y_PD, and y_SMC can achieve the path following in the presence of the time-varying interferences, demonstrating the effectiveness of addressing the internal uncertain and external disturbances by three controllers. Figure 10 shows that the oscillations of δ_MPC are less than δ_SMC, and the performance of δ_MPC is much better than δ_PD.

Figure 9.

The path following based on different controllers.

Figure 10.

The rudder angle and rate based on different controllers.

In general, the heading angle and rudder angle are constantly fluctuating to counteract the effects of time-varying disturbances. The rudder angle varies particularly sharply due to first-order high-frequency wave disturbances. The oscillations that occur in the control signals in the article are affected by the accuracy of the model predictions, but the main cause is the influence of external disturbances, and such oscillations are consistent with nautical reality.

Case 3. To prove the availability of handling the sideslip angle based on the designed nonlinear observer, we compare the built nonlinear observer with the linear observer used in Wang et al. ⁷ The desired waypoint: P₁(0, 0), P₂(3500, 0), P₃(6000, 2000), and P₄(8000, 2000). The initial states: u  =  7.0 m/s, v  =  0, r  =  0, ψ = 0 (x₀, y₀)  =  (0, 200 m). Other parameters are the same as Case 1. To compare the performance of the observer explicitly, we select a new current set ψ_c with a jump point. The new set is given

$$\{\begin{array}{c}{\psi }_c={30}^{\circ }+{30}^{\circ }\sin (0.025t),t<300\\ {60}^{\circ }+{30}^{\circ }\sin (0.025t),300\le t<600\\ {90}^{\circ }+{30}^{\circ }\sin (0.025t),600\le t<900\\ {120}^{\circ }+{30}^{\circ }\sin (0.025t),900\le t<1200\\ {60}^{\circ }+{30}^{\circ }\sin (0.025t),1200\le t<1500\\ {20}^{\circ }+{30}^{\circ }\sin (0.025t),t\ge 1500\end{array}$$

(28)

The simulation results based on the different sideslip observers are shown in Figures 11 and 12, and the subscript N represents the nonlinear observer, meanwhile, L represents the linear observer. In Figure 11, both the y_N and y_L can follow the reference path yd effectively in the presence of the interferences with jump point. Figure 12 illustrates that the sideslip angle β_N or β_L can be estimated successfully by the nonlinear or linear observer, respectively. However, the estimation error β_eN is smaller than β_eL, especially at the jump point of disturbances and the turning point of ship motion.

Figure 11.

The path following based on different observers.

Figure 12.

The sideslip angle and estimation errors.

In conclusion, the above simulation results show that the set side slip angle control algorithm can make the ship return to the set path stably and smoothly. The control algorithm designed in this paper can effectively solve the problem of model uncertainty and external disturbances. To achieve an accurate following of the path, the ship needs to keep steering in the presence of external disturbances, resulting in a certain amplitude of oscillation in both heading and rudder angles. However, the MPC oscillation amplitude is smaller and shorter at the beginning, and the variation is smoother, which shows the MPC algorithm's ability to handle rudder angle constraints and optimization.

Conclusions

This paper combines the backstepping algorithm with MPC methodology for path following of ships with the rudder constraints and optimization. A nonlinear observer is designed to compensate the sideslip angle, and the LESO technology is used to estimate the yaw rate and total unknowns. Furthermore, a new switching method is presented to reduce the rudder amplitude at the turning point. Finally, the compared simulation shows that the proposed controller can drive ship to follow the desired path effectively with a smaller and smoother rudder angle in the presence of time-varying interferences such as wind, wave, and current. The sideslip angle, yaw rate, and the total unknowns can be determined accurately. Moreover, the rudder behavior would be promoted successfully by the presented switching reference heading method.

The main researches work in the future are: (1) Applying the controller designed in the article to large-scale ship maneuvering simulators. (2) Apply the improved algorithm to the actual ship. (3) Research collaborative control algorithms for multiple ships to achieve consistency control, formation control, obstacle avoidance control, etc.