Robust predictive control for anti-rolling path following of underactuated surface vessels using adaptive Kalman filter

Abstract

With the strong time-varying disturbances in wave fields, the large sway motion of underactuated surface vessels can be caused leading to the control performance reduction and large range of roll angle, such that a novel path following control with roll constraints is proposed. Firstly, the affine switching system model is constructed with coupled surge, sway, roll, and yaw. Then, with augmented state adaptive Kalman filter, the estimations of system state and disturbance can be obtained for the feedforward compensation. Considering the existing estimation error between real and estimated states, the state estimation has been adopted directly in the controller design for the performance improvements. With the convex optimization of linear matrix inequalities, the invariant ellipsoid constraint ensures the close-loop system uniformly bounded and stable. Lastly, the simulation results have verified the effectiveness of the designed control approach.

Introduction

Path following of surface vehicles has attracted the attention of the control community for many years, and in the last few years, the motion control problem of typical and asymmetric underactuated surface vessels (USVs) has been a research focus. ^{1 –3} At present, most marine surface vessels are only equipped with the propeller and rudder. The three degrees of freedom movement consisting of horizontal position and heading angle is controlled by the longitudinal propulsion of the propeller and steering moment of the rudder. As the number of controlled variables is more than the control actuators, this kind of control system is a typical underactuated system. ^{4 –6} Due to the wide application area and theoretical value, the problem how to handle such control system remains to be a researching focus, and several researches have been done such as the smooth transition from fully actuated to underactuated configuration. ⁷ Many intelligent control methods have also been successfully employed, such as fuzzy logic control. ^{8 –10} The desired path is usually a straight line or a waypoint path, which consists of piecewise straight lines with the curvature being zero. It is common for USVs to navigate at or nearly at constant velocity for a long distance. ¹¹ For a kind of nonlinear system with incomplete acceleration constraints, it is difficult to design the path following controller, especially in the navigation with strong time-varying disturbances from sea waves. It is due to a major error of path following or even strenuous roll motion caused by disturbances. To cope with the complex unknowns especially for disturbances, great achievements have been reached, and many effective methods have been proposed in recent years. The adaptive approximation-based control approach was developed, and the asymptotic convergence of output regulation could be achieved with the close-loop system uniformly ultimately bounded. ¹² Furthermore, the finite-time stabilizing observer and controller could be combined with the tracking problem, and the proposed method had been verified by simulation and experiment in Wang et al.. ¹³

The strenuous roll motion can not only affect the normal performance of the ship equipment, cargo fixing, and crew efficiency but also lead to the damage or even capsizing of surface vessels. To avoid such unsafe situation and ensure the safety of high assets, intelligent risk analysis and critical decision were rendered to diagnose the fault, evaluate the risk, and make the critical decision. ¹⁴ For most USVs, the rudder can produce the yaw moment to control the heading angle and can also generate the roll torque to achieve the roll stabilization at the same time. Compared with other roll stabilization facilities, rudder roll stabilization has the advantages of less investment, smaller space occupation, and efficient roll stabilization performance. For years, the path following problem of USVs has been addressed separately in the robustness and roll stabilization efficiency in many publications. For the path following problem of USVs with roll constraint, a robust constrained predictive control method was proposed utilizing adaptive Kalman filter and disturbance observer. The control design has been theoretically analyzed, and simulations and experiments have been demonstrated to illustrate the control efficiency. ¹⁵

For the model parameter uncertainties and external disturbances, the heuristic fuzzy control was proposed to achieve the robust adaptive following along the desired path. ¹⁶ For the ship motion system with nonlinear hydrodynamic coefficients, the fuzzy observer design method was proposed to achieve the rudder roll stabilization using T-S fuzzy model. ¹⁷ In Wang et al., ¹⁸ aiming at tracking an uncertain surface vessel with unknown disturbances, an adaptive universe-based fuzzy control was proposed with the combination of retractable fuzzy partitioning leading the achievement of the global asymptotic convergence. The internal model control was combined with the neural network, and the learning function of neural network was adopted to improve the system robustness. ¹⁹ However, the coordination control of path following and roll stabilization remains to be an open problem. The only consideration for fast convergence of path following error, which means the high gain controller, often leads to large roll angle and rudder saturation. On the other hand, the extreme pursuit for roll stabilization efficiency can significantly reduce the path following performance, rapidity, and maneuverability. Thus, it is a significant issue for USVs to make the effective compromise between the path following performance and roll stabilization to achieve the good following performance and rudder roll stabilization.

As is known, model predictive control (MPC) can carry out the online optimization updated by actual system states leading to the performance improvement and enhanced system robustness. Furthermore, the ability of handling the system constraints explicitly in advance is another featured strategy, such that it is widely utilized in the practical process industries. ²⁰ In Li et al., ²¹ based on online constrained optimization of MPC strategy, the roll constraint considered as the state constraint was ensured in the safety range, but the control system robustness had not been discussed. To deal with disturbances and constraints, the predictive controller was combined with disturbance observer based on the nominal model. The system stability and robustness were ensured, but the control algorithm had the computation burden due to online optimization twice respectively for control and feedforward compensation variables in every control period. ²² However, the achieved application studies are sort of limited. The assumption has been made that full motion states of USVs can be obtained accurately. In the practical application environment, not all states can be measured, while only partial states measurements can be obtained mixed with random noise. As MPC cannot effectively restrain the external time-varying disturbances using the nominal system model, the close-loop system performance can be deteriorated.

To handle the above existing problems, a novel direct constrained robust predictive control combined with augmented state adaptive Kalman filter is proposed for USVs to achieve path following and rudder roll stabilization in this article. The main contributions of this article include:

1. For the measurement of the partial motion states mixed with external disturbances, adaptive Kalman filter with the augmented state is adopted for the convenience for the control design.

2. Considering the existing estimation error between actual and estimated states, the real measurements are utilized directly in the performance index and control design, such that a direct constrained robust predictive control method combined with adaptive Kalman filter can be proposed for the achievements of the better performance of path following and rudder roll stabilization.

This article is organized as follows. The objective of this article is presented in the “Problem statement” section. In the “Adaptive Kalman filter with augmented state” section, adaptive Kalman filter with the augmented state is proposed, followed by a novel robust predictive control method in the “Robust predictive control with rudder roll stabilization” section. The “Simulation study” section presents simulations, and conclusions are presented in the final section.

Problem statement

In the open literature, path following of USVs has been addressed with two different approaches: one is to design the controller directly using the error dynamics of path following in the inertial frame and the other is to research on it in the Serret–Frenet frame instead to simplify the controller design in some ways. The frame {SF} used for path following control of USVs is shown in Figure 1. The origin of the frame {SF} is located at the closest point on the curve Ω from the origin of the body-fixed frame {B}. Ω is the given target path, e is defined as the distance between the origins of {SF} and {B}, ψ_SF is the tangential direction of the given path, and ψ is the heading angle of the vessel.

Figure 1.

Illustration of the coordinates in the Serret–Frenet frame.

The error dynamics of path following in the Serret–Frenet frame is given by ²³

$\begin{array}{l}\dot{e}=u\text{sin}\left(\overline{\psi }\right)+v\text{cos}\left(\overline{\psi }\right)\\ \dot{\overline{\psi }}=\dot{\psi }-{\dot{\psi }}_{\text{SF}}=\frac{\kappa }{1-e\kappa }\left(u\text{sin}\left(\overline{\psi }\right)-v\text{cos}\left(\overline{\psi }\right)\right)+r\end{array}$

1

where $\overline{\psi }=\psi -{\psi }_{\text{SF}}$ is defined as the heading error; u, v, and r are the surge, sway, and yaw velocity, respectively. κ is the curvature of the given path. For most path following problems for USVs in open water, the path is often a straight line or a waypoint path which consists of piecewise straight lines with the curvature κ = 0, and the sway speed is rather small enough to be ignored, such that it can be assumed v = 0, then the vessel forward speed is $U=\sqrt{u^2+v^2}\approx u$ . Therefore, the error dynamics (1) can be simplified as

$\begin{array}{l}\dot{e}=u\text{sin}\left(\overline{\psi }\right)+v\text{cos}\left(\overline{\psi }\right)\approx U\text{sin}\left(\overline{\psi }\right)\\ \dot{\overline{\psi }}=r\end{array}$

2

The desired path following can be achieved by driving path following error e and heading angle error $\overline{\psi }$ converging to be zero. The USV stability can be severely affected by large roll motion, so roll angle should be constrained in the safety range. The roll dynamics must be established, and the four degrees of freedom nonlinear equations of motion (surge u, sway v, roll p, and yaw r) are given by ^22,24

$\begin{array}{l}m\left(\dot{u}-vr-X_G{r}^2+Z_{G}pr\right)=X\\ m\left(\dot{v}+ur+X_G\dot{r}-Z_G\dot{p}\right)=Y\\ I_{XX}\dot{p}-I_{XZ}\dot{r}-m{Z}_G\left(\dot{v}+ur\right)=\mathit{K}-mgG{M}_T\varphi \\ I_{ZZ}\dot{r}-I_{XZ}\dot{p}+m{Z}_G\left(\dot{v}+ur\right)=N\end{array}$

3

where m is the ship mass, I_XX and I_XZ are the coupled moment of inertia between x,x- and x,z-axes, respectively. p is the roll velocity and ϕ is the roll angle; $G{M}_T$ denotes the metacentric height; X_G and Z_G are the locations of the gravity center in the x-axis and z-axis, respectively. The hydrodynamic forces X,Y and moments K,N are usually the three order Taylor series polynomials, the definitions in detail can be referred in the literature. ²⁵

Since the model (3) is a nonlinear system with strong couplings and uncertainties, it is difficult to design the controller directly. Due to the limited working range of the rudder $\left(|\delta |\le {35}^{\circ }\right)$ , the motion equations can be divided into three working areas by linearizing the original nonlinear model at different equilibrium points to facilitate the model-based design. The original nonlinear model is described by the switching model, and those three working areas are given by

${\Omega }_1=\left\{\delta |-{10}^{\circ }\le \delta \le {10}^{\circ }\right\},{\Omega }_2=\left\{\delta |{10}^{\circ }\le \delta \le {35}^{\circ }\right\},{\Omega }_3=\left\{\delta |-{35}^{\circ }\le \delta \le -{10}^{\circ }\right\}$

4

The simplified switching model has the following formulation as

$\begin{array}{l}\dot{\mathit{x}}={\mathit{A}}_i\mathit{x}+{\mathit{B}}_i\delta +{\mathit{B}}_{d}di=1,2,3\\ \mathit{y}=\mathit{H}\mathit{x}+v\\ x\in X\\ \delta \in U\end{array}$

5

If $\delta \in {\Omega }_i$ , the system states can be switched into the linear model i given by

$\begin{array}{l}{\mathit{A}}_i=\left[\begin{array}{cccccc}0& 0& 1& 0& 0& 0\\ 0& a_{i11}& a_{i12}& a_{i13}& a_{i14}& 0\\ 0& a_{i21}& a_{i22}& a_{i23}& a_{i24}& 0\\ 0& a_{i31}& a_{i32}& a_{i33}& a_{i34}& 0\\ 0& 0& 0& 1& 0& 0\\ U\beta & 0& 0& 0& 0& 0\end{array}\right] {\mathit{B}}_i=\left[\begin{array}{c}0\\ b_{i1}\\ b_{i2}\\ b_{i3}\\ 0\\ 0\end{array}\right] \\ \mathit{H}=\left[\begin{array}{cccccc}0& 0& 1& 0& 0& 0\end{array}\right]\end{array}$

6

where the nonlinear model (3) is linearized around the equilibrium points corresponding to δ = 0, δ = 15, and $\delta =-15$ , denoted as $[{\mathit{A}}_1,{\mathit{B}}_1]$ , $[{\mathit{A}}_2,{\mathit{B}}_2]$ , and $[{\mathit{A}}_3,{\mathit{B}}_3]$ , respectively. The system state is defined as $\mathit{x}={\left[x_1,x_2,x_3,x_4,x_5,x_6\right]}^{\text{T}}={\left[\overline{\psi },v,r,p,\varphi ,e\right]}^{\text{T}}$ and $\beta =\text{sin}\left(\overline{\psi }\right)/\overline{\psi }$ . The coefficients $a_{ijm}$ (i = 1, 2, 3, j = 1, 2, 3, m = 1, 2, 3, 4) and b_ij (i = 1, 2, 3, j = 1, 2, 3) are constants, obtained by the experiment and CFD software. ²⁰ The motion state of USVs includes two motions of the low and high frequencies. As the motion of high frequency results to a slight movement of the vessel and the position remains unchanged, only the motion of low frequency is considered for four degrees of freedom surge, sway, roll, and yaw. d is the moment caused by the disturbances and satisfies $|d(t)|\le d_{\text{max}}$ , $|\dot{d}(t)|\le d_{\epsilon }$ . Here yaw angle is chosen as the measurement with random noise, assume that the noise obeys normal distribution for white noise to facilitate the research. $X\triangleq \left\{x||x_5|\le {\varphi }_{\text{max}}\right\}$ and $U\triangleq \left\{\delta ||\delta |\le {\delta }_{\text{max}}\right\}$ are denoted as the constraints for roll angle and rudder angle, respectively.

Assuming that the simplified model (5) is selected to represent the original nonlinear system with the range $|\overline{\psi }|\le \frac{\pi }{2}$ , β can be contained in the compact set $\left\{\beta \text{|}\frac{\text{2}}{\pi }\le \beta \le 1\right\}$ , then each linear switching model in equation (5) can be described by the convex hull of two following linear systems

$\left[{\mathit{A}}_i, {\mathit{B}}_i\right]\in \text{Co}\left\{\left[{\mathit{A}}_{i1}, {\mathit{B}}_i\right],\left[{\mathit{A}}_{i2}, {\mathit{B}}_i\right]\right\}, i=1,2,3$

7

where ${\mathit{A}}_{i1}$ and ${\mathit{A}}_{i2}$ are the system state matrices when $\beta =\frac{2}{\pi }$ and $\beta =1$ , respectively.

Therefore, the problem discussed in the article can be described as follows. On the basis of the switching models formulated in equation (4), with the consideration of measurement noise and external disturbances caused by wind, wave, and current, design adaptive Kalman filter with augmented state to eliminate the measurement noise and estimate the full motion state and disturbance to make the feedforward compensation $u_d(k)$ . With the input constraint, design direct constrained robust predictive control and obtain the compound controller $u(k)=u_{\text{mpc}}(k)+u_d(k)$ to achieve the desired path following, ensure roll angle in a safe constrained range, and restrain the effect of random disturbance d. So that the control system has a good performance and robustness.

Adaptive Kalman filter with augmented state

In the sophisticated sea conditions, USVs can be affected by the random disturbance d, system noise w, and measurement noise v, which makes it difficult to obtain the system state accurately or have a good performance of path following. Factually, the wave disturbance d can be considered as the Gaussian colored noise. The random disturbance d was considered as the white noise, and there existed the approximation error after adopting Kalman filtering method in the study by Yu. ²⁶ The colored noise, transferred into white noise, was estimated by shaping filter but this method had the computation burden. Due to the motion characteristic of low frequency, it can be approximately regarded as $\dot{d}\approx 0$ in the state estimation period. Taking the disturbance d as the augmented state, adaptive Kalman filtering method is adopted to make the optimal estimations of system state and disturbance to reduce the effects caused by random disturbance and measurement noise, and make the direct feedforward compensation for the estimated disturbance, so as to improve the system robustness and reduce the computation complexity of feedback control $u_{\text{mpc}}$ .

For the further simplification, with roll disturbance force and moment omitted, the yaw disturbance moment is mainly considered. After the first-order Eulerian discretization of the augmented dimensional switching model (5), the equations can be formulated as

$\begin{array}{l}\overline{\mathit{x}}\left(k+1\right)={\overline{\mathit{A}}}_i\overline{\mathit{x}}\left(k\right)+{\overline{\mathit{B}}}_i\delta \left(k\right)+w\left(k\right)\\ {\overline{\mathit{y}}}_1\left(k\right)=\overline{\mathit{H}}\overline{\mathit{x}}\left(k\right)+v\left(k\right)\end{array}$

8

where $\overline{\mathit{x}}\left(k+1\right)=\left[\begin{array}{l}\widehat{x}\left(k+1\right)\\ \widehat{d}\left(k+1\right)\end{array}\right]={\left[\begin{array}{ccccccc}\widehat{\overline{\psi }}& \widehat{v}& \widehat{r}& \widehat{p}& \widehat{\varphi }& \widehat{e}& \widehat{d}\end{array}\right]}^{\text{T}}$ , T is the sampling period.

$\begin{array}{l}{\overline{\mathit{A}}}_i\in \text{Co}\left\{{\mathit{I}}_{7*7}+T\left[\begin{array}{ccccccc}0& 0& 1& 0& 0& 0& 0\\ 0& a_{i11}& a_{i12}& a_{i13}& a_{i14}& 0& 0\\ 0& a_{i21}& a_{i22}& a_{i23}& a_{i24}& 0& 1\\ 0& a_{i31}& a_{i32}& a_{i33}& a_{i34}& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ \frac{2U}{\pi }& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right],{\mathit{I}}_{7*7}+T\left[\begin{array}{ccccccc}0& 0& 1& 0& 0& 0& 0\\ 0& a_{i11}& a_{i12}& a_{i13}& a_{i14}& 0& 0\\ 0& a_{i21}& a_{i22}& a_{i23}& a_{i24}& 0& 1\\ 0& a_{i31}& a_{i32}& a_{i33}& a_{i34}& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ U& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right]\right\}\\ {\overline{\mathit{B}}}_i=T{\left[\begin{array}{ccccccc}0& b_{i1}& b_{i2}& b_{i3}& 0& 0& 0\end{array}\right]}^{\text{T}},i=1,2,3\end{array}$

$\overline{\mathit{H}}$ is the augmented measurement matrix.

The detailed process of filtering and estimations for state and disturbance is as follows

$\overline{\mathit{x}}\left(k|k-1\right)=\overline{\mathit{A}}\overline{\mathit{x}}\left(k-1\right)+\overline{\mathit{B}}u\left(k\right)$

9

$\overline{\mathit{P}}\left(k|k-1\right)=\overline{\mathit{A}}\overline{\mathit{P}}\left(k-1|k-1\right){\overline{\mathit{A}}}^{\text{T}}+\overline{\mathit{Q}}\left(k\right)$

10

$\overline{\mathit{K}}\left(k\right)=\frac{\overline{\mathit{P}}\left(k|k-1\right){\overline{\mathit{H}}}^{\text{T}}}{\overline{\mathit{H}}\overline{\mathit{P}}\left(k|k-1\right){\overline{\mathit{H}}}^{\text{T}}+\mathit{R}\left(k\right)}$

11

$\overline{\mathit{x}}(k|k)=\overline{\mathit{x}}(k|k-1)+\overline{\mathit{K}}(k)({\overline{y}}_1(k)-\overline{\mathit{H}}\overline{\mathit{x}}(k|k-1))$

12

$\widehat{d}\left(k\right)=\overline{C}\overline{\mathit{x}}\left(k|k\right)$

13

$\overline{\mathit{P}}\left(k|k\right)=\left(\mathit{I}-\overline{\mathit{K}}\left(k\right)\overline{\mathit{H}}\right)\overline{\mathit{P}}\left(k|k-1\right)$

14

where $\overline{C}=\left[\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& 1\end{array}\right]$ .

The disturbance d, estimated in equation (13), mainly affects the yaw dynamics, so the feedforward compensation u_d is formulated as

$u_d(k)=-\frac{\widehat{d}(k)}{b_{i2}}$

15

Due to the sophisticated wave condition, the system and measurement noise cannot be determined in advance. The adaptive Sage–Husa filtering method, based on time-varying disturbance statistics, is adopted to estimate the motion state. The variance matrices of system and measurement noise can be estimated and corrected constantly to prevent the filter divergence caused by the significant error of estimated state. Based on Sage–Husa adaptive filter, not only the state can be estimated using the measurements but also real-time statistics for different noise can be done. The filtering method still has the advantages of simple form, high filtering precision, and fast-response speed. To obtain the real-time estimations of system and measurement noise matrices $\overline{\mathit{Q}},\overline{\mathit{R}}$ online and prevent the filtering divergence caused by the noise uncertainties, forgetting factor b is introduced to limit the memory length of the filter, and the effect of new measurements on the current estimation is added. The formulations are given by ²⁷ :

$\begin{array}{l}\overline{\mathit{R}}(k)=\left[1-z\left(k-1\right)\right]\mathit{R}\left(k-1\right)+z\left(k-1\right)\left[\mathit{v}\left(k\right){\mathit{v}}^{\text{T}}\left(k\right)-\overline{\mathit{H}}\overline{\mathit{P}}\left(k|k-1\right){\overline{\mathit{H}}}^{\text{T}}\right]\\ \overline{\mathit{Q}}(k)=\left[1-z\left(k-1\right)\right]\overline{\mathit{Q}}\left(k-1\right)+z\left(k-1\right)\left[\overline{\mathit{K}}\left(k\right)\mathit{v}\left(k\right){\mathit{v}}^{\text{T}}\left(k\right){\overline{\mathit{K}}}^{\text{T}}\left(k\right)+\overline{\mathit{P}}\left(k\right)\right]\end{array}$

16

where $z(k-1)=\frac{1-b}{1-b^k}$ is denoted as the weighting coefficient, the value range of forgetting factor b is usually $0.95\le b\le 0.99$ , and $\mathit{v}(k)=\overline{\mathit{y}}(k)-\overline{\mathit{H}}\overline{\mathit{x}}(k)$ is defined as the new information sequence.

Robust predictive control with rudder roll stabilization

The strong disturbances caused by wind, wave, and current can lead to the large sway motion of USVs, then roll angle can be out of the safety constrained range. Thus, it is necessary to consider the fast convergence of following error and safety range of roll angle. With the advantage of dealing with state and input constraints, predictive control method with rudder roll stabilization is proposed in this section.

From the previous section, it can be concluded that there exists the observation error between the real and observation motion states after the state estimation of adaptive Kalman filtering method and disturbance compensation. The estimation error still exists in the disturbance compensation. Therefore, adaptive Kalman filtering should be taken into consideration dealing with the state and input constraints using predictive control. Considering the integration design of estimation and control, a direct constrained robust predictive control is proposed based on the state estimation and mixed ${\mathit{H}}_2/{\mathit{H}}_{\infty }$ approach method. The estimation state is used directly in the control design and performance index to present the good path following performance with constrained roll angle and restrained disturbance.

After the feedforward compensation adopted in the model (5), the discrete predictive model with the disturbance can be formulated as follows ^28,29

$\begin{array}{l}\widehat{\mathit{x}}\left(k+1\right)={\mathit{A}}_i\widehat{\mathit{x}}\left(k\right)+{\mathit{B}}_i{u}_{\text{mpc}}\left(k\right)+{\mathit{B}}_{\mathit{d}}{d}_1\left(k\right)+\mathit{K}\left(k\right)[{\overline{\mathit{y}}}_1\left(k\right)-\mathit{H}\widehat{\mathit{x}}\left(k\right)]\\ \mathit{y}\left(k\right)=\left[\begin{array}{c}\mathit{H}\widehat{\mathit{x}}\left(k\right)\\ \mathit{R}{u}_{\text{mpc}}\left(k\right)\end{array}\right]\end{array}$

17

where $\mathit{y}(k)$ is the system output, $\widehat{\mathit{x}}\left(k\right)=\widehat{\mathit{x}}\left(k|k\right)$ represents the system state measured at time k, $x\left(k+i\right)=x\left(k+i|k\right)$ refers to the state at time k + i predicted based on the measurements at time k, and $u\left(k\right)=u\left(k|k\right)$ refers to the control move to be implemented at time k. ${\mathit{A}}_i=\mathit{M}{\overline{\mathit{A}}}_i{\mathit{M}}^{\text{T}}$ , ${\mathit{B}}_i=\mathit{M}{\overline{\mathit{B}}}_i$ , $\mathit{M}=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\end{array}\right]$ , $i=1,2,3$ . $\mathit{K}=\mathit{M}\overline{\mathit{K}}$ is the Kalman filter gain, the submatrix of the previous gain $\overline{\mathit{K}}$ in the augmented system. $d_1\left(k\right)=d\left(k\right)-\widehat{d}\left(k\right)$ and $\mathit{e}\left(k\right)=x\left(k\right)-\widehat{\mathit{x}}\left(k\right)$ are the estimation errors of disturbance compensation and state, respectively, and both of them are Euclidean norm-bounded. ³⁰ Equivalently, ${\Vert d_1\Vert }_2^2={\displaystyle \sum _{k=0}^{\infty }{d}_1^{\text{T}}\left(k\right)d_1\left(k\right)}\le {\varpi }^2,\varpi >0$ and ${\Vert \mathit{e}\Vert }_2^2={\displaystyle \sum _{k=0}^{\infty }{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)}\le {\xi }^2,\xi >0$ .

With a constant state-feedback control law $u_{\text{mpc}}\left(k\right)=\mathit{F}\widehat{\mathit{x}}\left(k\right)$ , the discrete predictive model (17) can be transformed into

$\begin{array}{l}\widehat{\mathit{x}}\left(k+1\right)=\mathit{\varphi }\widehat{\mathit{x}}\left(k\right)+{\mathit{B}}_d{d}_1\left(k\right)+\mathit{K}\left(k\right)\mathit{H}\mathit{e}\left(k\right)\\ y\left(k\right)=\left[\begin{array}{c}\mathit{H}\widehat{\mathit{x}}\left(k\right)\\ \mathit{R}{u}_{\text{mpc}}\left(k\right)\end{array}\right]\end{array}$

18

where $\mathit{\varphi }\left(k\right)={\mathit{A}}_i\left(k\right)+{\mathit{B}}_i\left(k\right)\mathit{F}$ .

For a good path following performance and robust stability, the mixed $H_2/H_{\infty }$ approach method is proposed, and the indexes are defined as follows:

• (1) The ${\mathit{H}}_2$ constraint: the system output $y(k)$ satisfies

${\Vert y\Vert }_2^2={\displaystyle \sum _{k=0}^{\infty }{y}^{\text{T}}\left(k\right)y\left(k\right)}={\displaystyle \sum _{k=0}^{\infty }\left[{\widehat{\mathit{x}}}^{\text{T}}\left(k\right){\mathit{Q}}_1\widehat{\mathit{x}}\left(k\right)+u_{\text{mpc}}^{\text{T}}\left(k\right){\mathit{R}}_{\text{1}}{u}_{\text{mpc}}\left(k\right)\right]}\le {\alpha }^2$

19

where ${\mathit{Q}}_1={\mathit{H}}^{\text{T}}\mathit{H}$ and ${\mathit{R}}_1={\mathit{R}}^{\text{T}}\mathit{R}$ refer to weighting matrices of system state and input, respectively.

• (2) The ${\mathit{H}}_{\infty }$ constraint: the transfer function $\mathit{G}$ from the external disturbance to the system output satisfies

${\Vert \mathit{G}\Vert }_{\infty }\le \gamma$

20

Consider a Lyapunov function $V\left(\widehat{\mathit{x}}\right)={\widehat{\mathit{x}}}^{\text{T}}\mathit{P}\widehat{\mathit{x}}$ , $\mathit{P}>0$ , then

$\begin{array}{l}V\left(\widehat{\mathit{x}}\left(k+1\right)\right)-V\left(\widehat{\mathit{x}}\left(k\right)\right)\\ ={\left[\mathit{\varphi }\left(k\right)\widehat{\mathit{x}}\left(k\right)+{\mathit{B}}_d{d}_1\left(k\right)+\mathit{K}(k)\mathit{H}\mathit{e}\left(k\right)\right]}^{\text{T}}\mathit{P}\left[\mathit{\varphi }\left(k\right)\widehat{\mathit{x}}\left(k\right)+{\mathit{B}}_d{d}_1\left(k\right)+\mathit{K}\left(k\right)\mathit{H}\mathit{e}\left(k\right)\right]-{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)\mathit{P}\widehat{\mathit{x}}\left(k\right)\\ =\left[\begin{array}{ccc}{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)& d_1^{\text{T}}\left(k\right)& {\mathit{e}}^{\text{T}}\left(k\right)\end{array}\right]\mathit{N}\left(k\right){\left[\begin{array}{ccc}{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)& d_1^{\text{T}}\left(k\right)& {\mathit{e}}^T\left(k\right)\end{array}\right]}^{\text{T}}\\ +{\gamma }^2\left[d_1^{\text{T}}\left(k\right)d_1\left(k\right)+{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)\right]-y^{\text{T}}\left(k\right)y\left(k\right)\end{array}$

21

where

$\mathit{N}\left(k\right)=\left[\begin{array}{ccc}\chi & {\mathit{\varphi }}^{\text{T}}\left(k\right)\mathit{P}{\mathit{B}}_d& {\mathit{\varphi }}^{\text{T}}\left(k\right)\mathit{P}\mathit{K}\left(k\right)\mathit{H}\\ {\mathit{B}}_d^{\text{T}}\mathit{P}\mathit{\varphi }\left(k\right)& {\mathit{B}}_d^{\text{T}}\mathit{P}{\mathit{B}}_d-{\gamma }^2\mathit{I}& {\mathit{B}}_d^{\text{T}}\mathit{P}\mathit{K}\left(k\right)\mathit{H}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\left(k\right)\mathit{P}\mathit{\varphi }\left(k\right)& {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\left(k\right)\mathit{P}{\mathit{B}}_d& {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\left(k\right)\mathit{P}\mathit{K}\left(k\right)\mathit{H}-{\gamma }^2\mathit{I}\end{array}\right]$

and $\chi ={\mathit{\varphi }}^{\text{T}}\left(k\right)\mathit{P}\mathit{\varphi }\left(k\right)-\mathit{P}+{\mathit{H}}^{\text{T}}\mathit{H}+{\mathit{F}}^{\text{T}}{\mathit{R}}^{\text{T}}\mathit{R}\mathit{F}$

Due to the system stability, it can be assumed that ${{\displaystyle \text{lim}}}_{k\to \infty }\widehat{\mathit{x}}(k)=0$ , and then

$\begin{array}{c}{\displaystyle \sum _{k=0}^{\infty }\left[V\left(\widehat{\mathit{x}}\left(k+1\right)\right)-V\left(\widehat{\mathit{x}}\left(k\right)\right)\right]}={\displaystyle \sum _{k=0}^{\infty }\left[{\widehat{\mathit{x}}}^{\text{T}}\left(k+1\right)\mathit{P}\widehat{\mathit{x}}\left(k+1\right)-{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)\mathit{P}\widehat{\mathit{x}}\left(k\right)\right]}\\ =-{\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)\end{array}$

22

Adding the both sides of equation (21) from $k=0$ to $k\to \infty$ , then using equations (19) and (22), and carrying out a simple manipulation, we have

$\begin{array}{c}{\Vert y\left(k\right)\Vert }_2^2={\displaystyle \sum _{k=0}^{\infty }\left[y^{\text{T}}\left(k\right)y\left(k\right)\right]}={\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)+{\displaystyle \sum _{k=0}^{\infty }\left[{\gamma }^2{d}_1^{\text{T}}\left(k\right)d_1\left(k\right)+{\gamma }^2{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)\right]}\\ +{\displaystyle \sum _{k=0}^{\infty }\left[\begin{array}{c}{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)d_1^{\text{T}}\left(k\right){\mathit{e}}^{\text{T}}\left(k\right)\end{array}\right]\mathit{N}\left(k\right){\left[\begin{array}{c}{\widehat{\mathit{x}}}^{\text{T}}\left(k\right)d_1^{\text{T}}\left(k\right){\mathit{e}}^{\text{T}}\left(k\right)\end{array}\right]}^{\text{T}}}\end{array}$

23

Referring the definition of the ${\mathit{H}}_{\infty }$ constraint, it is the transfer function from the disturbance to the system output. Thus, with the consideration of the system state error and the assumption $\widehat{\mathit{x}}(0)=0$ , equation (20) is equivalent to

${\displaystyle \sum _{k=0}^{\infty }\left[{\mathit{y}}^{\text{T}}\left(k\right)\mathit{y}\left(k\right)-{\gamma }^2{d}_1^{\text{T}}\left(k\right)d_1\left(k\right)-{\gamma }^2{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)\right]}\le 0$

24

From equation (23), we can conclude that if and only if $\mathit{N}(k)\le 0$ such that equation (24) is satisfied. Equivalently, if $\mathit{N}(k)\le 0$ exists, equation (24) is satisfied and so the constraint (20) is.

Using Schur complement, the condition $\mathit{N}(k)\le 0$ is equivalent to

$\left[\begin{array}{cccccc}\mathit{Q}& 0& 0& \mathit{Q}{\mathit{A}}_i^{\text{T}}+{\mathit{Y}}^{\text{T}}{\mathit{B}}_i^{\text{T}}& \mathit{Q}{\mathit{H}}^T& {\mathit{Y}}^{\text{T}}{\mathit{R}}^{\text{T}}\\ 0& {\gamma }^2{\alpha }^2\mathit{I}& 0& {\alpha }^2{\mathit{B}}_d^{\text{T}}& 0& 0\\ 0& 0& {\gamma }^2{\alpha }^2\mathit{I}& {\alpha }^2{\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}& 0& 0\\ {\mathit{A}}_i\mathit{Q}+{\mathit{B}}_i\mathit{Y}& {\alpha }^2{\mathit{B}}_d& {\alpha }^2\mathit{K}\mathit{H}& \mathit{Q}& 0& 0\\ \mathit{H}\mathit{Q}& 0& 0& 0& {\alpha }^2\mathit{I}& 0\\ \mathit{R}\mathit{Y}& 0& 0& 0& 0& {\alpha }^2\mathit{I}\end{array}\right]\ge 0$

25

where $\mathit{Q}={\alpha }^2{\mathit{P}}^{-1}$ , $\mathit{F}=\mathit{Y}{\mathit{Q}}^{-1}$ , $i=1,2,3$ .

On the premise that equation (25) is satisfied, equation (23) can be written as

$\begin{array}{c}{\displaystyle \sum _{k=0}^{\infty }\left[y^{\text{T}}(k)y(k)\right]}\le {\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)+{\displaystyle \sum _{k=0}^{\infty }\left[{\gamma }^2{d}_1^{\text{T}}\left(k\right)d_1\left(k\right)+{\gamma }^2{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)\right]}\\ ={\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)+{\gamma }^2{\Vert d_1\Vert }_2^2+{\gamma }^2{\Vert \mathit{e}\Vert }_2^2\\ \le {\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)+{\gamma }^2{\varpi }^2+{\gamma }^2{\xi }^2\end{array}$

26

If the ${\mathit{H}}_2$ constraint (19) is satisfied, we have

${\displaystyle \sum _{k=0}^{\infty }\left[{\mathit{y}}^{\text{T}}\left(k\right)\mathit{y}\left(k\right)\right]} \le {\widehat{\mathit{x}}}^{\text{T}}(0)\mathit{P}\widehat{\mathit{x}}(0)+{\gamma }^2{\varpi }^2+{\gamma }^2{\xi }^2\le {\alpha }^2$

27

Using Schur complement, equation (27) is equivalent to

$\left[\begin{array}{ccc}1& {\widehat{\mathit{x}}}^{\text{T}}(0)& {\left[{\gamma }^2\left({\varpi }^2+{\xi }^2\right)\right]}^{\text{T}}\\ \widehat{\mathit{x}}(0)& \mathit{Q}& 0\\ {\gamma }^2\left({\varpi }^2+{\xi }^2\right)& 0& {\gamma }^2{\alpha }^2\left({\varpi }^2+{\xi }^2\right)\end{array}\right]\ge 0$

28

• (3) The invariant ellipsoid constraint

It can be concluded from equation (21) that ${\gamma }^2{d}_1^{\text{T}}\left(k\right)d_1\left(k\right)+{\gamma }^2{\mathit{e}}^{\text{T}}\left(k\right)\mathit{e}\left(k\right)-{\mathit{y}}^{\text{T}}\left(k\right)y\left(k\right)<0$ cannot be guaranteed if only equation (25) is satisfied, so that the considered Lyapunov function adopted cannot be guaranteed to decrease monotonically. Based on the effect of the disturbance $d_1$ , the constraint of invariant ellipsoid sets is proposed to ensure the state $\widehat{\mathit{x}}(k+1)$ in the set if $\widehat{\mathit{x}}(k)$ is. Equivalently

$\widehat{\mathit{x}}{\left(k+1\right)}^{\text{T}}{\mathit{Q}}^{-1}\left(k+1\right)\widehat{\mathit{x}}\left(k+1\right)\le \widehat{\mathit{x}}{\left(k\right)}^T{\mathit{Q}}^{-1}\left(k\right)\widehat{\mathit{x}}\left(k\right)\le 1$

Thus, the invariant ellipsoid constraint is formulated as

$\mathit{Q}\left(k+1\right)\le \mathit{Q}\left(k\right)$

29

• (4) The input constraint: Physical limitations inherent in the rudder invariably impose hard constraints on the manipulated variable. Considering the previous designed feedforward compensation u_d, the input constraint for predictive control move $u_{\text{mpc}}$ can be represented by

$\begin{array}{l} u_{n\text{max}}=u_{n\text{max}}-|u_d|\\ {\Vert u_{\text{mpc}}\left(k+i\right)\Vert }_2^2\le u_{n\text{max}}^2,i\ge 0\end{array}$

30

Due to the previous state-feedback control, we have $u_{\text{mpc}}=\mathit{F}\widehat{\mathit{x}}=\mathit{Y}{\mathit{Q}}^{-1}\widehat{\mathit{x}}$ , and then equation (30) can be written as

$\underset{i\ge 0}{\text{max}}{\Vert u_{\text{mpc}}\left(k+i\right)\Vert }_2^2=\underset{i\ge 0}{\text{max}}{\Vert Y{\mathit{Q}}^{-1}\widehat{\mathit{x}}\left(k+i\right)\Vert }_2^2\le \underset{z\in \xi }{\text{max}}{\Vert \mathit{Y}{\mathit{Q}}^{-1}z\Vert }_2^2={\lambda }_{\text{max}}\left({\mathit{Q}}^{-1/2}{\mathit{Y}}^{\text{T}}\mathit{Y}{\mathit{Q}}^{-1/2}\right)$

31

Using Schur complement, equation (31) is equivalent to

$\left[\begin{array}{cc}{u}_{n\text{max}}^2& Y\\ Y^{\text{T}}& \mathit{Q}\end{array}\right]\ge 0$

32

• (5) The roll angle constraint

The large sway motion of USVs can be caused by external strong disturbances resulting to the bad control performance and stability. So it is necessary to take the roll angle constraint into account, which can be represented by

${\Vert \varphi \left(k+i\right)\Vert }_2^2\le {\varphi }_{\text{max}}^2,i\ge 1$

33

Selecting $\varphi \left(k+i\right)={\mathit{H}}_1\widehat{\mathit{x}}\left(k+i\right)$ and ${\mathit{H}}_1=[0 0 00 1 0]$ , it means that the roll angle is the system output ${\widehat{\mathit{x}}}_5$ , then equation (33) is equivalent to

${\widehat{\mathit{x}}}^{\text{T}}\left(k+1\right){\mathit{H}}_1^{\text{T}}{\mathit{H}}_1\widehat{\mathit{x}}\left(k+1\right)\le {\varphi }_{\text{max}}^2$

34

For the left side of equation (34), we have

$\begin{array}{c}{\widehat{\mathit{x}}}^{\text{T}}\left(k+1\right){\mathit{H}}_1^{\text{T}}{\mathit{H}}_1\widehat{\mathit{x}}\left(k+1\right)={\left[\begin{array}{c}{\mathit{Q}}^{-1/2}x\left(k\right)\\ d_1\left(k\right)\\ \mathit{e}\left(k\right)\end{array}\right]}^{\text{T}}\left[\begin{array}{c}{\mathit{Q}}^{-1/2}{\mathit{\varphi }}^{\text{T}}\left(k\right)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\end{array}\right]{\mathit{H}}_1^{\text{T}}{\mathit{H}}_1{\left[\begin{array}{c}{\mathit{Q}}^{-1/2}{\mathit{\varphi }}^{\text{T}}(k)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\end{array}\right]}^{\text{T}}\left[\begin{array}{c}{\mathit{Q}}^{-1/2}x\left(k\right)\\ d_1\left(k\right)\\ \mathit{e}\left(k\right)\end{array}\right]\\ \le \left(1+{\varpi }^2+{\xi }^2\right){\lambda }_{\text{max}}\left({\mathit{H}}_1{\left[\begin{array}{c}\mathit{Q}{\mathit{\varphi }}^{\text{T}}\left(k\right)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\end{array}\right]}^{\text{T}}\left[\begin{array}{ccc}{\mathit{Q}}^{-1}& & \\ & \mathit{I}& \\ & & \mathit{I}\end{array}\right]\left[\begin{array}{c}\mathit{Q}{\mathit{\varphi }}^{\text{T}}\left(k\right)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^T{\mathit{K}}^T\end{array}\right]{\mathit{H}}_1^{\text{T}}\right)\end{array}$

35

where ${\lambda }_{\text{max}}(\ast )$ is the largest eigenvalue.

For $\forall k>0$ , there exists ${\widehat{\mathit{x}}}^{\text{T}}\left(k\right){\mathit{Q}}^{-1}\widehat{\mathit{x}}\left(k\right)\le 1$ . Then, we have

${\Vert \left[\begin{array}{c}{\mathit{Q}}^{-1/2}\widehat{\mathit{x}}\left(k\right)\\ d_1\left(k\right)\\ \mathit{e}\left(k\right)\end{array}\right]\Vert }^2={\left[\begin{array}{c}{\mathit{Q}}^{-1/2}\widehat{\mathit{x}}\left(k\right)\\ d_1\left(k\right)\\ \mathit{e}\left(k\right)\end{array}\right]}^{\text{T}}\left[\begin{array}{c}{\mathit{Q}}^{-1/2}\widehat{\mathit{x}}\left(k\right)\\ d_1\left(k\right)\\ \mathit{e}\left(k\right)\end{array}\right]\le \left(1+{\varpi }^2+{\xi }^2\right)$

If equation (33) or (34) is satisfied, the roll angle can be constrained in the given safety range. If there exists

${\widehat{\mathit{x}}}_5^2\mathit{I}-\left(1+{\varpi }^2+{\xi }^2\right)\left({\mathit{H}}_1{\left[\begin{array}{c}\mathit{Q}{\mathit{\varphi }}^{\text{T}}\left(k\right)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\end{array}\right]}^{\text{T}}\left[\begin{array}{ccc}{\mathit{Q}}^{-1}& & \\ & \mathit{I}& \\ & & \mathit{I}\end{array}\right]\left[\begin{array}{c}\mathit{Q}{\mathit{\varphi }}^{\text{T}}\left(k\right)\\ {\mathit{B}}_d^{\text{T}}\\ {\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}\end{array}\right]{\mathit{H}}_1^{\text{T}}\right)\ge 0$

36

and then equation (34) is satisfied.

Using Schur complement, equation (36) can be written as

$\left[\begin{array}{cc}\frac{\mathit{Q}}{\left(1+{\varpi }^2+{\xi }^2\right)}& \mathit{Q}{\mathit{A}}_i^{\text{T}}{\mathit{H}}_1^{\text{T}}+Y^{\text{T}}{\mathit{B}}_i^{\text{T}}{\mathit{H}}_1^{\text{T}}\\ {\mathit{H}}_1{\mathit{A}}_i\mathit{Q}+{\mathit{H}}_1{\mathit{B}}_{i}Y& {\widehat{\mathit{x}}}_5^2\mathit{I}-\left(1+{\varpi }^2+{\xi }^2\right){\mathit{H}}_1{\mathit{B}}_d{\mathit{B}}_d^{\text{T}}{\mathit{H}}_1^{\text{T}}+{\mathit{H}}_1\mathit{K}\mathit{H}{\mathit{H}}^{\text{T}}{\mathit{K}}^{\text{T}}{\mathit{H}}_1^{\text{T}}\end{array}\right]\ge 0$

37

where $\mathit{Q}={\alpha }^2{\mathit{P}}_{ }^{-1}$ , $\mathit{F}=Y{\mathit{Q}}_{ }^{-1}, i=1,2,3$ .

Remark

If ${\mathit{B}}_d=\mathit{K}=0$ corresponding to omitting the disturbance and estimation error, equation (30) can be reduced to the similar formulation in Olalla et al. ³¹

Theorem 1

With the constraints described in the linear matrix inequalities (LMIs) (15), (25), (28), (29), (32), and (37), there exists the control law $u=u_{\text{mpc}}+u_d$ that ensures the asymptotic robust stability of the close-loop system (18).

Proof

From the above derivations, LMI (15) is the control move for the disturbance compensation. The monotonic decrease of $V(\widehat{\mathit{x}}\left(k\right))$ is ensured by LMIs (25) and (29). The constraints for performance index, input, and roll angle are ensured by LMIs (28), (32), and (37) respectively. On the premise of the theorem feasibility at time k, assuming that the optimal solution is $\left\{{\mathit{Q}}_0,\dots ,{\mathit{Q}}_{\mathit{N}-1},Y_0,\dots ,{\mathit{F}}_0,\dots ,{\mathit{F}}_{\mathit{N}-1},{\alpha }^{\ast }\left(k\right)\right\}$ . From equations (25) and (29), it can be concluded that there exists the feasible solution $\{{\mathit{Q}}_1,\cdots ,{\mathit{Q}}_{\mathit{N}-1},{\mathit{Q}}_{\mathit{N}-1},Y_1,\cdots ,Y_{\mathit{N}-1},Y_{\mathit{N}-1},{\mathit{F}}_1,\cdots ,{\mathit{F}}_{\mathit{N}-1},{\mathit{F}}_{\mathit{N}-1},{\alpha }^{\ast }(k)\}$ with ${\alpha }^{\ast }\left(k\right)=\alpha \left(k+1\right)$ at time $k+1$ . From the optimization theory, there exists ${\alpha }^{\ast }\left(k\right)\ge \alpha \left(k+1\right)\ge {\alpha }^{\ast }\left(k+1\right)$ .

In addition, from equations (21) and (26), without the disturbance and estimation error, the monotonic decrease of $V(\widehat{\mathit{x}}\left(k\right))$ , which is equivalent to $V(\widehat{\mathit{x}}\left(k+1\right))<V(\widehat{\mathit{x}}\left(k\right))$ , can still be satisfied using the invariant ellipsoid constraint (30). Thus, the close-loop system (18) is asymptotically robust stable.

Simulation study

To demonstrate the effectiveness of the proposed algorithm, several simulations have been done. The model parameters $\left[{\mathit{A}}_i,{\mathit{B}}_i\right],i=1,2,3$ used for model ship, which can be referred in the study by Wang et al. ¹⁸ for details. The surge speed is constant at U = 7 m/s. The parameters of controller and cost function are selected as follows: ${\mathit{B}}_d=0.02*{\left[\begin{array}{cccccc}0& 0& 1& 0& 0& 0\end{array}\right]}^{\text{T}}$ , $\mathit{H}=\text{eye}(6)$ , $\mathit{R}=0.02$ , and $\gamma =2$ . The control period is set to T = 0.5 s. The initial state is $\left(\overline{\psi },v,r,p,\varphi ,\mathit{e}\right)=\left(-45/57.3,0,0,0,0,20\right)$ . The constraints of roll angle and rudder angle are $|\varphi |<5^{\circ }$ and $|\delta |\le {35}^{\circ }$ , respectively. Assuming that the measurement noise is the stationary Gaussian white noise with the mean value about 0 and variance about 0.01. Considering the non-dimensional disturbance $d(t)=0.02\left(\text{sin}\left(0.1t\right)+\text{cos}\left(0.1t+\frac{\pi }{4}\right)+\text{cos}\left(0.1t+\frac{\pi }{6}\right)\right)$ , the simulation results are shown in Figures 2 to 10.

Figure 2.

Heading angle error.

Figure 3.

Tracking error.

Figure 4.

Roll angle.

Figure 5.

Heading angle error.

Figure 6.

Tracking error.

Figure 7.

Roll angle.

Figure 8.

Rudder angle.

Figure 9.

Estimated disturbance.

Figure 10.

Path following comparison.

Figures 2 to 4 show the heading angle error, tracking error, and roll angle under the condition with or without the proposed adaptive Kalman filter. The measurement noise has almost eliminated showing the filtering efficiency. Considering the existence of external disturbance, Figures 5 to 9 compare the control performance using the proposed algorithm and typical MPC, respectively. Figures 5 and 6 are the comparison of heading angle error and tracking error. Compared with the typical algorithm, using the disturbance compensation of the proposed method in this article, both heading angle error and tracking error have shown the faster convergence speed and smaller convergence range. In Figure 7, the convergence range of roll angle obtained using the disturbance compensation is apparently smaller with the state constraint. Figure 8 is the comparison of rudder angle. With the disturbance compensation, the curve of rudder angle is smoother and has a smaller fluctuation range with the input constraint. The smaller fluctuation range can result to an application significance, as the violent fluctuation affects the service life of the rudder and reduce the vessel safety. Figure 9 compares the disturbance estimation with actual disturbance. The estimated value is in a small fluctuation range of the actual value, showing a good estimation performance.

In Figure 10, path following comparison is presented, and it can be seen that the vessel can reach the desired path in a shorter time and has a smaller fluctuation range even with the disturbance, demonstrating the better performance of the proposed algorithm.

From the above figures, it can be concluded that adaptive Kalman filtering method can obtain the system disturbance and facilitate the control design, such that the predictive control with the disturbance compensation can achieve the faster convergence speed with the smaller convergence ranges of roll and rudder angles. The smaller the convergence range means the more safety condition and better performance of restraining the disturbance.

Conclusions

For the path following issue of USVs, the large sway motion and rudder saturation caused by time-varying disturbances can lead to the constraints and robustness problem. Considering the control performance, stability, and safety, a direct constrained predictive control with roll constraints is proposed for path following based on augmented Kalman filtering method. Utilizing the switching modeling method, the model complexity has been reduced facilitating the control design. The adaptive Kalman filter with the augmented state is adopted for the estimations of state and stochastic disturbance to eliminate the measurement noise and obtain the disturbance compensation. With the proposed predictive control, all considered constraints can be transformed into the convex optimization of LMIs using Schur complement. The effectiveness of the proposed control method has been illustrated by simulations.