Path following for underactuated surface vessels with disturbance compensating predictive control

Abstract

It is common that underactuated surface vessels sailing on the sea suffer from strong external sea disturbances, such that the large roll motion can be probably caused resulting to the bad performance of path following. In order to realize the coordinate control of the rudder roll stabilization and path following, a robust controller with roll constraints is designed by combining predictive control with disturbance observer. Firstly, the rudder angle operating range is divided in the Serret–Frenet coordination frame, and then linear models corresponding to different equilibrium can be established with heading control and roll dynamics. Secondly, considering external disturbance and unmodeled error as the lumped disturbances, the disturbance observer can be utilized to achieve the real-time feedforward compensation leading to the improved system robustness. Thirdly, the output redefinition method is adopted to transfer the original system into a new minimum phase system. Based on the receding horizon and state prediction strategies of the predictive control, the analytical control law can be obtained, and the linear programming method is used to guarantee the roll constraint for the simplified computational burden. Lastly, the simulations have been carried out to show the good performance of the proposed algorithm with the heading control and roll reduction.

Introduction

Referring to the marine problems, path following has been a mainstream direction due to military and civilian necessities in recent years. Numerous control methods have been designed to tackle with kinds of vessels and detailed practical difficulties, such as feedback control, neural network control, and adaptive control technique. ^{1 –5} The sailing vessels are mostly underactuated, as the rudder is usually used as the control input to manipulate the heading, yaw velocity, and lateral displacement, and the surge velocity is maintained by the propeller. For this kind of manipulation strategy, there exist several issues which ought to be considered: the underactuated characteristic and the external disturbances in the sea. ⁶ It is worth mentioning that it is not unique for practical application, as in real robotics applications, the publications ^7,8 related to nonholonomic systems can reveal the same issue. When the vessel sails on the sea, it is inevitable to take the external disturbance into account. Among all vessel states of surge, sway, heave, yaw, pitch, and roll, the roll motion caused by stochastic wave disturbance is the most significant. The severe roll motion can cause the seasickness, reduce the performance of the shipping equipment and crew members, and even lead to the shipwreck. ⁹ The disturbances cannot be omitted in the control design, and cited publications above have taken into account the disturbance in the design process. In addition to disturbances, the tracking problem in the presence of sensor noise showed a more general condition for practical application. ¹⁰

As is well known, the rudder plays the role in the heading control. However, the roll moment generated by the rudder operation can be offset by the disturbance moment caused by the wave disturbance. Further, the respond period of roll motion is much longer than that of the yaw motion, and it makes the rudder roll stabilization feasible in the physical and control aspects. Compared with some other roll damping methods, without the equipped appendages, it has the advantages of low cost, small occupation space, and good roll reduction performance. ¹¹ As mentioned, owing to the chief role of the rudder in the heading control, the two goals of heading control and roll reduction cannot be considered separately, which implies that strong coupling exists. On account of the facts that the simple pursuit of roll reduction can deteriorate the path following performance and much more consideration about following the path can lead to excessive roll motion, it is urge to make the appropriate trade-off between these two tasks. In the study by Zhou and Wang, ¹² based on the different rudder frequency responses of roll and yaw, the H∞ control method was utilized for the realization of path following and rudder roll stabilization separately, and then the weighting functions can be chosen with the weighting function method. The path following can be accurately achieved with the external disturbance reduction, but the good trade-off performance had not been realized due to the transfer function-based separated control design method for path following and rudder roll stabilization. Using T-S fuzzy model, the nonlinear system of the yaw and roll motions can be obtained, and with the principle of parallel distribution compensation, a fuzzy control method was designed to improve the performance of roll reduction, but the path following performance had not been taken into account although it had already been affected by the roll reduction operation. ¹³ The advantages of fuzzy logic system in dealing with uncertainties have made it popular. The direct fuzzy control in Yu et al. ¹⁴ was proposed to compensate for effects caused by actuator saturation, but a new fuzzy sliding mode control method had achieved the same objective in Chu et al. ¹⁵ Actually, not only sliding mode control but also backstepping method can be integrated. Based on this idea, ¹⁶ the 3-D path following was realized. Furthermore, on the premise that the presented publications have pushed forward the theoretical research, many scholars and engineers in the shipping industry and related robotics area have put their best efforts to achieve big progress in practical implementation. ^17,18 Based on pure sliding model control and its combination with smooth control method, the effectiveness and applicability of this algorithm has been demonstrated in numerical simulation and practical experiment. ¹⁹ So far, these publications have concerned the tracking problem of single object, but with the rapid economic development, the control problem of multiple surface vessels have emerged. Many well-developed method and theories cannot be applied directly due to the peculiar researching difficulty, and then some bionics-based methods have grown and spread out rapidly. Based on swarm control, the collision avoidance and distribution for multiple vessels had been achieved. ^20,21

For the process control with slow dynamics, model predictive control (MPC) is one of the effective methods to tackle the multi-objective optimal control with multivariable couplings and constraints, due to its receding horizon strategy and explicit consideration for constraints in the open loop optimization. Recently, besides the wide application in the industrial control, MPC has been appeared in the motion control with reduced computational burden and improved hardware condition. ^22,23 Based on the nonlinear model with roll dynamics, roll angle was considered as the state constraint, such that the safety range of the roll angle can be ensured with online constrained optimization of MPC property. The heading control and roll reduction had been achieved, but the time varying disturbances had not been discussed. ²⁴ To deal with constraints and disturbances, the nonlinear disturbance observer was adopted with the nominal model, such that the system stability and robustness can be guaranteed. ²⁵ In Zhang et al., ²⁶ a path following method was proposed with the combination of Kalman filter, disturbance observer, and robust constrained MPC. With partial measured states, the full states were estimated, and the control system robustness had been significantly improved due to the disturbance compensation in the control design. However, the computational burden was much large due to the online optimization for feedforward compensation and control with the application of linear matrix inequalities.

Inspired by Zhang et al., ²⁶ a robust path following control method is proposed in order to realize the coordinate control of path following and rudder roll stabilization. The main objectives of this study include:

With the disturbance compensating strategy, a robust MPC control is proposed for the nonlinear models of surface vessels with roll constraint in order to achieve the path following and rudder roll stabilization simultaneously.

Considering the non-minimum phase dynamics of the control system, ⁵ the output redefinition method is adopted to transfer the original system into a new minimum phase system.

For the computation reduction, analytical control law and linear programming method have been utilized.

This article is organized as follows. The objective of this article is briefly presented in the second section. In the third section, the disturbance compensating strategy is given, followed by a novel robust predictive control method in the fourth section. Several simulations have been carried out in the fifth section, and conclusions are in last section.

Problem statement

For path following control of underactuated surface vessels (USVs), two different methods have been addressed: (1) control design with error dynamics in the inertial coordination frame; (2) the adoption of the Serret–Frenet frame. The latter one is often utilized for the simplification of the control design. In Figure 1, the frame {SF} used for path following control is shown. Ω is the given desired path, e is the distance between the origins of {SF} and the body-fixed frame {B}, ${\psi }_{SF}$ is the desired tangential direction of the given path, and ψ is the heading angle.

Figure 1.

Illustration of path following in the Serret–Frenet frame.

The error dynamics of path following in the Serret–Frenet frame are given by

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

1

where $\overline{\psi }=\psi -{\psi }_{SF}$ is defined as the heading error. u, v, and r are the surge, sway, and yaw, respectively. κ is the curvature of the given path. For most path following issues in open water, the desired path is a straight line or a waypoint path consisted of piecewise straight lines with $\kappa =0$ . Assumed that the sway velocity is small enough to be neglected, which means $v\approx 0$ . Then, the error dynamics (1) can be simplified as

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

2

For the error dynamics, it is well understood that the control objective of path following is to drive path error e and heading angle error $\overline{\psi }$ to be zero. From the above analyses, the USV performance and stability can be severely affected by large roll motion, so the roll dynamics must be established to ensure the roll angle in the safety range. The coupled nonlinear dynamics with four degrees of freedom of surge u, sway v, roll p, and yaw r are given by ²²

$\begin{array}{l}m(\dot{u}-rv-x_G{r}^2+z_{G}rp)={\displaystyle \sum X}+{\tau }_x\\ m(\dot{v}+ur+x_G\dot{r}-z_G\dot{p})={\displaystyle \sum Y}+{\tau }_y\\ I_{xx}\dot{p}-I_{xz}\dot{r}-m{z}_G(ur+\dot{v})\text{=}{\displaystyle \sum K}+{\tau }_k\\ I_{zz}\dot{r}-I_{xz}\dot{p}+m{x}_G(ur+\dot{v})={\displaystyle \sum N}+{\tau }_n\end{array}$

3

where m is the vessel mass, $I_{xx}$ and $I_{zz}$ are the inertial moment about x- and z-axes, and $I_{xz}$ is the coupled inertial moment of x- and z-axes. x_G and z_G are the locations of the gravity center in the x-axis and z-axis. p is the roll angular velocity, and ϕ is the roll angle. The hydrodynamic forces X,Y and moments K,N are usually third-order Taylor serial polynomials of nonlinear hydrodynamic constants, and specific definitions can be referred in Fossen. ²⁷ ${\tau }_x$ , ${\tau }_y$ , ${\tau }_k$ , and ${\tau }_n$ are denoted as different hydrodynamic effects due to the control input.

Since the dynamic model (3) is a nonlinear system with uncertainties and strong couplings, it is difficult to be directly used in the control design. Due to the limited operating range of the control input produced by the rudder, it is assumed that there are three operating equilibrium points corresponding to different operation modes in the whole operating envelope for the USV. ²⁷ With the linearization of the original nonlinear model, the linear models with disturbances at different operating modes can be represented as

$\begin{array}{l}\dot{x}={\overline{\mathit{A}}}_{i}x+{\overline{\mathit{B}}}_i\delta +{\mathit{B}}_d\mathit{d},\delta \in {\Omega }_i\\ \mathit{x}\in \mathit{X}=\{\mathit{x}||{\mathit{x}}_5|\le {\varphi }_{\text{max}}\}\end{array}$

4

$\overline{{\mathit{A}}_i}=\left[\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& a_{i11}& a_{i12}& a_{i13}& a_{i14}\\ 0& a_{i21}& a_{i22}& a_{i23}& a_{i24}\\ 0& a_{i31}& a_{i32}& a_{i33}& a_{i11}\\ 0& 0& 0& 1& 0\end{array}\right]\text{, }{\overline{\mathit{B}}}_i=\left[\begin{array}{l}0\\ b_{i1}\\ b_{i2}\\ b_{i3}\\ 0\end{array}\right]$

5

where $\mathit{x}={[\overline{\psi },v,r,p,\varphi ]}^{\text{T}}$ is defined as the system state, and ${\overline{\mathit{A}}}_i,{\overline{\mathit{B}}}_i$ are the system state and input matrix of the ith subsystem. d is the lumped disturbances of the external wave disturbance and unmodeled error caused by piecewise linearization. δ is the rudder angle served as the control input, and $\mathit{X}=\{\mathit{x}||{\mathit{x}}_5|\le {\varphi }_{\text{max}}\}$ is the roll angle constraint. ${\Omega }_i$ is the operation area of the ith subsystem corresponding to the different operation range with $i\in \overline{N}=\{1,2,3\}$ , and the ith model and related parameters should be utilized if $\mathit{x}\in {\Omega }_i$ . The adjacent area between the ith and jth subsystems is denoted as ${\Omega }_{ij}$ , which is represented as

${\Omega }_{ij}=\{\mathit{x}(k)|\mathit{x}(k)\in {\Omega }_i,\mathit{x}(k+1)\in {\Omega }_j,i,j\in \overline{N}\}$

6

Based on the polyhedral modeling method, referring to the rudder angle range, the three subsets can be described by the polyhedrons as follows

$\begin{array}{l}{\Omega }_1=\{\delta |-{35}^{\circ }\le \delta \le -{10}^{\circ }\},\\ {\Omega }_2=\{\delta |{10}^{\circ }\le \delta \le {35}^{\circ }\},\\ {\Omega }_3=\{\delta |-{10}^{\circ }\le \delta \le {10}^{\circ }\}\end{array}$

7

where these three subsets are referred to two different path following conditions about the heading changing and heading keeping, respectively.

When the USV sails on the sea, the external sea waves can excite the large roll motion which can deteriorate the performance and stability. In order to verify the better disturbance compensation efficiency, it is obviously necessary to simulate the wave-induced motion accurately. The response amplitude operator-based method, namely response amplitude operator (RAO) can be utilized to realize the motion simulation. Conventionally, the vessel motion caused by hydrodynamic forces and moments can be simulated in the inertial coordination frame. Different from this method, the steady forward motion can also be simulated in the inertial seakeeping coordination frame, which coincides with the other common coordination frame, body-fixed coordination frame at the average time. ²⁸ Obviously, the two simulated motion results obtained from two different frames cannot directly be used to carry out some simple mathematical manipulations. As the roll motions obtained from the two frames are in the same direction, the roll motion can then be used to achieve the roll simulation in the system model, such that the performance and stability of the proposed method in this article can be demonstrated.

In the wave theories, the wave is considered to be consisted of the long-crested waves to facilitate the study, such that the wave amplitude can be assumed to be the sum of some harmonic waves. The assumption can be described as

$\zeta (t)={\displaystyle \sum _{i=1}^N{\zeta }_i\text{cos}(w_{i}t-{\epsilon }_i)}$

8

where $\zeta (t)$ is the wave amplitude. ${\zeta }_i$ and w_i are the amplitude and angular frequency of the ith harmonic wave, respectively. ${\epsilon }_i$ represents the random phase of the ith harmonic wave, and it is randomly distributed at the range $[0,2\pi ]$ .

The wave energy spectrum is the key for the wave issues, the wave spectrum recommended by International Towering Tank Conference is used in this article, and it is formulated as

$S(w)=\frac{0.0081g^2}{w^2}\text{exp}(-\frac{3.11}{{{\rm H}}_{1/3}^2{w}^5})$

9

where S is the power density, g is the gravity, and w is the angular frequency of the wave. ${{\rm H}}_{1/3}$ is denoted as the significant wave height, and it is the parameter used to represent the sea conditions. With the adopted wave spectrum, the amplitude of the harmonic wave can be described as

${\zeta }_i=\sqrt{2S(w_i)\Delta w}$

10

where $\Delta w$ is the fixed range of the angular frequencies in the wave spectrum. Substituting (10) into (8), the wave amplitude can then be obtained.

As RAO is the linear transfer function from wave to vessel motion, there exists the linear property, such that with the wave slope related RAO, the vessel motion can be formulated as

${\eta }_j(t)={\displaystyle \sum _{i=1}^N{k}_i{\zeta }_i}|G_{i,j}|\text{cos}(w_{i}t+{\epsilon }_i+\angle G_{i,j})$

11

where ${\eta }_j$ is the amplitude of the jth motion. $k_i=2\pi /\lambda =w_i^2/g$ is the wave number, and λ is the wave length. $|G_{i,j}|$ and $\angle G_{i,j}$ are the amplitude and phase of jth RAO motion.

After the wave analysis, considering the advantages of robust MPC in tacking with constrained multivariable optimization problems, it is adopted to achieve the coordinate control of path following and rudder roll stabilization. Then, the analyzed problem in this article can be described as: For the path following model (2) with ship dynamics (4), taking the wave disturbances on the sea into account, an integrated controller $\delta (k)={\delta }_{mpc}(k)+{\delta }_{dc}(k)$ is designed, such that the desired heading angle can be followed with the constrained roll angle and restrained disturbances in order to improve the system stability and robustness. For the designed controller, ${\delta }_{mpc}(k)$ and ${\delta }_{dc}(k)$ are referring to control inputs produced by nominal MPC and disturbance compensation, respectively.

Disturbance compensating design

As the time-varying disturbances and unmodeled errors exist, the control input ${\delta }_{mpc}(k)$ obtained based on the nominal system model prediction cannot eliminate the steady-state error in path following, such that it can lead to the poor system robustness. Referring to the separation principle, the disturbance estimation can be independent to the control design, and then the disturbance estimation is firstly implemented.

In order to deal with the existed disturbance in the control system, the disturbance observer is derived and developed. ¹⁵ With the designed disturbance observer, the real-time disturbance in the system can be estimated with input and output states of the control system, and then the feedforward compensation can be carried out to eliminate the effects caused by unknown disturbances, unmodeled dynamics, and other factors. Besides, the disturbance observer still has the advantages of simple structure and low computational complexity, such that the disturbance compensation mechanism has been widely used in many areas. Referring to Yang et al., ²⁹ the linear disturbance observer is designed as follows

$\begin{array}{l}\dot{z}=-\mathit{L}{\mathit{B}}_d(z+\mathit{L}\mathit{x})-\mathit{L}({\overline{\mathit{A}}}_i\mathit{x}+{\overline{\mathit{B}}}_i{\delta }_{mpc})\\ \widehat{d}=z+\mathit{L}\mathit{x}\end{array}$

12

where z is the auxiliary state, $\mathit{L}=\text{diag}\{l_1,l_2,l_3,l_4,l_5\}$ is the gain matrix of the linear observer, and $\widehat{d}$ is the estimation of the lumped disturbance d. The structure diagram of the disturbance observer is shown in Figure 2.

Figure 2.

Structure diagram of disturbance observer.

The estimation error of the lumped disturbance is defined as

$e_d=d-\widehat{d}$

13

For the accurate estimation of the lumped disturbance, the disturbance estimation $\widehat{d}$ tends to equal to the lumped disturbance d. Equivalently, the estimation error defined as (9) tends to be zero. The derivative of the estimation error is given by

${\dot{e}}_d=\dot{d}-\dot{\widehat{d}}=\dot{d}-\mathit{L}{\mathit{B}}_d{e}_d$

14

According to the ship kinematics, the motion state of the USV is the composition of the low-frequency and high-frequency motions. Among the two kinds of motions, the high-frequency motion can only result to the slight motion with the average position unchanged, such that the high-frequency motion can be neglected, and only the low-frequency motion is considered. Then, it is naturally assumed that the lumped disturbance remains unchanged during the sampling periods with $\stackrel{.}{d}\approx 0$ . If the observer gain L is properly chosen to make $-\mathit{L}{\mathit{B}}_d$ the Hurwitz matrix, the estimation error can then tend to be zero.

In order to implement the disturbance compensation mechanism, the stability analysis of the designed disturbance observer is necessary. The Lyapunov function of the linear disturbance observer is chosen as

$V(e_d)=\frac{1}{2}{e}_d^T{e}_d$

Using the estimation error derivative (14), the derivative of the selected Lyapunov function can be obtained as

$\dot{V}=e_d^T{\dot{e}}_d=e_d^T(\dot{d}-\dot{\widehat{d}})=-e_d^T\mathit{L}{\mathit{B}}_d{e}_d$

15

If there exists $\mathit{L}{\mathit{B}}_d>0$ , then $\dot{V}=-e_d^T\mathit{L}{\mathit{B}}_d{e}_d<0$ can be guaranteed resulting to the convergence of the designed disturbance observer. The detailed theoretical analyses can be referred in Li and Sun. ²³ As the nominal MPC is a discrete model-based optimization method, a discrete form of the disturbance observer (12) is needed in order to maintain the consistency. Using the first-order Euler discretization method, the discrete disturbance observer can be formulated as

$\begin{array}{l}z(k+1)={\mathit{A}}^{\prime }z(k)+{{\mathit{B}}^{\prime }}_1\mathit{x}(k)+{{\mathit{B}}^{\prime }}_2{\delta }_{mpc}(k)\\ \widehat{d}(k)=z(k)+\mathit{L}x(k)\end{array}$

16

where ${\mathit{A}}^{\prime }=-\mathit{L}{\mathit{B}}_{d}T+I$ , ${{\mathit{B}}^{\prime }}_1=-\mathit{L}({\mathit{B}}_d\mathit{L}+\overline{\mathit{A}})T$ , and ${{\mathit{B}}^{\prime }}_2=-\mathit{L}\mathit{B}T$ . T is denoted as the sampling time, and I represents the identity matrix with the appropriate dimension 5 here.

As the control matrix B _d is not the square matrix due to the underactuated characteristics of the linearized model (4), the feedforward compensation is only considered for the yaw moment, and the time-varying compensation is designed as

${\delta }_d=-\frac{{\widehat{d}}_c}{b_{i2}}$

17

Constrained MPC control design

With the wave disturbance on the sea, the disturbance observer has been designed to estimate the lumped disturbance. Further, in order to realize the coordinate control of heading control and roll rudder stabilization, the constrained MPC method combined with the disturbance observer is proposed to ensure the system stability and performance with the safety range of roll angle. As is known, MPC is an effective method to deal with the multivariable problem with the state and input constraints.

To facilitate the control design, the discrete form of the nominal system model can be formulated as

$\mathit{x}(k+1+i|k)={\mathit{A}}_i\mathit{x}(k+i|k)+{\mathit{B}}_i\delta (k+i|k),\delta \in {\Omega }_i$

18

where $\mathit{x}(k+i|k)$ and $\delta (k+i|k)$ are the state and control move predicted at time $k+i$ based on the measurements at time k. A _i and B _i are discrete system states matrix and control input matrix respectively corresponding to the operation range. With the first order Euler discretization, it can be obtained that ${\mathit{A}}_i={\overline{\mathit{A}}}_{i}T+I$ and ${\mathit{B}}_i={\overline{\mathit{B}}}_{i}T$ .

With the discrete nominal model (18), the optimization in this article can be described as

$\begin{array}{l}\text{min}J\\ s.t.\{\begin{array}{l}\mathit{x}(k+1+i|k)={\mathit{A}}_i\mathit{x}(k+i|k)+{\mathit{B}}_i\delta (k+i|k)\\ \mathit{x}\in \mathit{X}=\{\mathit{x}||\varphi |\le {\varphi }_{\text{max}}\}\\ \delta \in \mathit{U}=\{\delta ||\delta |\le {\delta }_{i\text{max}}\}\end{array}\end{array}$

19

where $\mathit{x}\in \mathit{X}=\{\mathit{x}||\varphi |\le {\varphi }_{\text{max}}\}$ refers to the state constraint of roll angle, and $\delta \in \mathit{U}=\{\delta ||\delta |\le {\delta }_{i\text{max}}\}$ is the control input constraint of rudder angle. The objective function J with the quadratic form can be defined as ³⁰

$J={\displaystyle \sum _{i=1}^{N_p}{[y_d(k+i)-y(k+i)]}^{\text{T}}\mathit{Q}[y_d(k+i)-y(k+i)]}+{\displaystyle \sum _{i=1}^{N_c}\mathit{R}{[\delta (k+i-1|k)]}^2}$

20

where $y_d=\mathit{C}{\mathit{x}}_d$ is the desired system output, and x _d is the desired system state with the coordination of roll angle and path following error. N_p and N_C are the prediction and control horizons respectively, generally subject to $N_C\le N_P$ . The positive semi-definite matrix Q and positive definite R are the gain matrices of system state and control input with $\mathit{Q}\ge 0$ and $\mathit{R}\ge 0$ .

With the designed disturbance observer and nominal model optimization description, the coordinate control design of rudder roll stabilization and path following can then be realized in the following processes.

System analysis of minimum phase characteristics

The non-minimum characteristic of the control system in this article is caused by the coordination control of heading control and roll reduction only with rudder served as the control input. Different from the general control system, the control system with the non-minimum phase has the unstable internal dynamics resulting to the increased difficulty for the control design. In order to tackle with such control systems, the output redefinition method can be adopted to ensure the stability of the internal system dynamics by choosing a new output and approaching the original system.

The output matrix C in the objective function (20) has the significant influence on the system stability. If all system states information is available, equivalently $\mathit{C}=[\begin{array}{cccc}1& 1& 1& \begin{array}{cc}1& 1\end{array}\end{array}]$ , the root locus of the system transfer function is shown in Figure 3. From Figure 3, it can be seen that in the right half plane, there exists the zero-pole meaning that the control system is of non-minimum phase. With the original output matrix, the designed predictive control law can make the system root enter the right half plane resulting to the divergence and instability of the control system.

Figure 3.

Root locus of non-minimum phase system.

Due to the limited capacity for the non-minimum phase system using the predictive control method, the output matrix is redefined as $\mathit{C}=[\begin{array}{cccc}\begin{array}{cc}1& 0\end{array}& 1& 0& 0\end{array}]$ meaning the integration of the measured system states $\overline{\psi }$ and r. The system root locus is shown in Figure 4. The system root is in the left half plane and the system is stable, such that the designed predictive control can ensure the system stability.

Figure 4.

Root locus of minimum phase system.

Nominal predictive control law design

Based on the nominal state space model (18), the system state prediction at the prediction horizon N_p can be obtained as

$\mathit{x}(k+N_p|k)={\mathit{A}}_i^{N_p}\mathit{x}(k)+{\mathit{A}}_i^{N_p-1}\mathit{B}\delta (k)+\dots +{\mathit{A}}_i^{N_p-N_c}\mathit{B}\delta (k+N_c-1)$

21

In order to facilitate the description, some sets are denoted as

$\mathit{W}={[\begin{array}{cccc}{y}_d(k+1)& y_d(k+2)& \cdots & y_d(k+N_P)\end{array}]}^T\mathit{F}={[\begin{array}{cccc}\mathit{C}{\mathit{A}}_i& \mathit{C}{\mathit{A}}_i{}^2& \cdots & \mathit{C}{\mathit{A}}_i{}^{N_P}\end{array}]}^T$

$\mathit{Y}={[\begin{array}{cccc}y(k+1)& y(k+2)& \cdots & y(k+N_P)\end{array}]}^T\mathit{U}={[\delta \begin{array}{cccc}(k)& \delta (k+1)& \cdots & \delta (k+N_c-1)\end{array}]}^T$

$\overline{\mathit{Q}}=\left[\begin{array}{cccc}\mathit{Q}& & & \\ & \mathit{Q}& & \\ & & \ddots & \\ & & & \mathit{Q}\end{array}\right]\overline{\mathit{R}}=\left[\begin{array}{cccc}\mathit{R}& & & \\ & \mathit{R}& & \\ & & \ddots & \\ & & & \mathit{R}\end{array}\right]\mathrm{\Gamma }=\left[\begin{array}{cccc}\mathit{C}{\mathit{B}}_i& 0& \cdots & 0\\ \mathit{C}{\mathit{A}}_i{\mathit{B}}_i& \mathit{C}{\mathit{B}}_i& \cdots & 0\\ ⋮& ⋮& \ddots & 0\\ \mathit{C}{\mathit{A}}_i{}^{N_P-1}{\mathit{B}}_i& \mathit{C}{\mathit{A}}_i{}^{N_P-2}{\mathit{B}}_i& \cdots & \mathit{C}{\mathit{A}}_i{}^{N_P-N_{\mathit{C}}}{\mathit{B}}_i\end{array}\right]$

With the simplified denotations, the cost function and the state prediction can then be rewritten as

$\begin{array}{l}J={(\mathit{W}-\mathit{Y})}^T\overline{\mathit{Q}}(\mathit{W}-\mathit{Y})+{\mathit{U}}^T\overline{\mathit{R}}\mathit{U}\\ \mathit{Y}=\mathit{F}\mathit{x}+\mathit{\Gamma }\mathit{U}\end{array}$

22

For the unconstrained optimization problem, the necessary condition $\frac{\partial J}{\partial U}=0$ for the minimization of the objective function can be utilized to obtain

$\mathit{U}={({\mathit{\Gamma }}^T\overline{\mathit{Q}}\mathit{\Gamma }+\overline{\mathit{R}})}^{-1}{\mathit{\Gamma }}^T\overline{\mathit{Q}}(\mathit{W}-F\mathit{x})$

23

With the receding horizon optimization strategy, the current control (the first row of the optimal control sequence) can be applied to act on the system, and it can be obtained as

${\delta }_{mpc}(k)=[\mathit{I}\text{ 0 0 }\dots \hspace{0.17em}\text{0]}\mathit{U}={\mathit{K}}_1-{\mathit{K}}_2\mathit{x}(\mathit{k})$

24

where ${\mathit{K}}_1={({\mathit{\Gamma }}^T\overline{\mathit{Q}}\mathit{\Gamma }+\overline{\mathit{R}})}^{-1}{\mathit{\Gamma }}^T\overline{\mathit{Q}}\mathit{W}$ and ${\mathit{K}}_2={({\mathit{\Gamma }}^T\overline{\mathit{Q}}\mathit{\Gamma }+\overline{\mathit{R}})}^{-1}{\mathit{\Gamma }}^T\overline{\mathit{Q}}F$ .

Constraints handling

Aiming at the reduced complexity and calculation burden caused by the constraints optimization using sequential quadratic programming, the linear programming method is adopted here introducing a new variable α, and the constraints handling problem can be formulated as

$\begin{array}{l}\underset{{\delta }_{mpc}(k)}{\text{min}}|\alpha |\\ s.t.\{\begin{array}{l}\alpha ={\delta }_{mpc}(k)-{\mathit{K}}_1+{\mathit{K}}_2\mathit{x}(k)\\ \varphi \in \mathit{X}\\ \delta (k)\in \mathit{U}\end{array}\end{array}$

25

where $\delta (k)={\delta }_{mpc}(k)+{\delta }_{dc}(k)$ is denoted as the integrated control with the nominal predictive control and the disturbance compensation. ${\delta }_{dc}(k)={\delta }_d(k)T$ , and T is the sampling period.

From the optimization algorithm (25), it can be seen that the model choice is determined by the current control. In order to avoid the system vibration caused but the frequent model switching, a soft model switching approach is utilized, and the historical errors are taken into account. Then, at time k, it can be formulated as

$e_{yi}(k)=y(k)-\widehat{y}(k),i=1,2,3$

26

where $e_{yi}(k)$ is the residual between the system output measured from the actual vessel and output prediction with the ith model at kth instant. The performance objective can be formulated to indicate the matching degree of the USV model and utilized predictive model at the instant time, which is defined as

$J_y={\displaystyle \sum _{h=1}^{\mathit{L}}{\theta }_h}{e}_y(k-h)$

27

where h is historical window length, and ${\theta }_h$ is the weighting factor. The smaller J_y equals to the bigger matching degree, such that the predictive model corresponding to the smallest $J_y(k)$ is chosen to be implemented at the current time control. At the next time $k+1$ , the predictive model is re-chosen referring to the recalculation of $J_y(k+1)$ .

Simulations

The parameters in the switching models of the USV can be referred in the literature, ²⁴ and then the matrices of system state and input in (4) can be selected as

${\overline{\mathit{A}}}_1={\overline{\mathit{A}}}_3=\left[\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& -0.02276& -2.7910& -0.009211& -0.1169\\ 0& -0.0009168& -0.1068& 0.009949& 0\\ 0& 0.002032& -0.3058& -0.01982& -0.04486\\ 0& 0& 0& 1& 0\end{array}\right],$

${\overline{\mathit{A}}}_2=\left[\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& -0.04062& -0.1899& -0.06664& -0.09348\\ 0& 0.0001167& -0.1468& 0.007198& -0.000828\\ 0& 0.002594& -0.3051& -0.01434& -0.04471\\ 0& 0& 0& 1& 0\end{array}\right],$

${\overline{\mathit{B}}}_1={\overline{\mathit{B}}}_3={\left[\begin{array}{cccc}\begin{array}{cc}0& -0.05699\end{array}& 0.002838& 0.004081& 0\end{array}\right]}^T,$

and ${\overline{\mathit{B}}}_2={\left[\begin{array}{cccc}\begin{array}{cc}0& -0.04575\end{array}& 0.002279& 0.003277& 0\end{array}\right]}^T$ . Assumed that the initial system state is $\mathit{x}=({25}^{\circ },0,0,0,0^{\circ })$ . The surge speed is constant at $u=7$ m/s, and the encounter angle is ${90}^{\circ }$ . The range constraints of roll and rudder angle are defined as $\mathit{X}=\left\{\mathit{x}||\varphi |\le {15}^{\circ }\right\}$ and $\mathit{U}=\left\{\delta ||\delta |\le {35}^{\circ }\right\}$ , respectively. The sampling period is 0.08 s. Due to the adopted single step prediction, the weighting matrices of the cost function are selected as $\mathit{Q}=400$ and $\mathit{R}=0.8$ , respectively. In order to evaluate the robustness and feasibility, the original nonlinear model is utilized in this section. Furthermore, due to the similar optimization strategy designed for rudder roll stabilization, the adaptive LQG-based control method ⁵ has been utilized for the comparison of the control performance in different conditions. The simulation results are as follows:

1. If the system measurement matrix is chosen as $\mathit{C}=[\begin{array}{cccc}\begin{array}{cc}1& 1\end{array}& 1& 1& 1\end{array}]$ , it can be obtained from the analyses in the “System analysis of minimum phase characteristics” subsection that the non-minimum phase characteristic of the control system can lead to the system divergence. It can be seen from Figure 5 that the heading angle errors are seriously divergent with different rates using both control methods, and this navigation condition can severely deteriorate the control system performance, even cause the safety accidents on the seaway.

2. In order to facilitate the study, the original system can be transformed into an equivalent one with the non-minimum phase characteristic by the output matrix redefinition $\mathit{C}=[\begin{array}{cccc}\begin{array}{cc}1& 0\end{array}& 1& 0& 0\end{array}]$ . Without the external disturbance, the simulation results for the performance comparison are shown in Figures 6 to 8. It can be seen that the heading angle error is convergent with the satisfactory ranges of roll angle and rudder angle. Using the predictive control proposed in this article, the heading error can convergent to the stable state faster, and overshoot of roll and rudder angles can become smaller resulting to the better control performance and longer service life of rudder. With the comparison between Figure 5 and Figure 6, the effectiveness of the output redefinition method can be illustrated.

3. With the significant wave height ${{\rm H}}_{1/3}=0.8\hspace{0.17em}\text{m}$ , the simulated wave amplitude has been shown in Figure 9 based on the wave spectrum theory, and comparison of heading angle error with or without the disturbance compensation mechanism is depicted in Figure 10. From Figure 10, the simulation results without and with the disturbance compensation are represented with the blue line and red line, respectively. From the blue line without the disturbance compensation, it can be seen that the existence of the wave disturbance can generate the significant error of the heading angle illustrating the disability of nominal predictive control in tackling with the time-varying disturbance. Compared with the red line, with the implementation of the disturbance observer, the steady error of heading angle tracking has been significantly reduced. It can be seen clearly from Figures 11 and 12 that the oscillations in the curves of roll and rudder angles have been decreased remarkably decreased revealing the necessity validation for the combination of disturbance compensation in present of external wave disturbances.

Figure 5.

Heading angle error.

Figure 6.

Heading angle error.

Figure 7.

Roll angle.

Figure 8.

Rudder angle.

Figure 9.

Random wave generated by the wave energy spectrum.

Figure 10.

Heading angle error comparison without/with the disturbance compensation.

Figure 11.

Roll angle comparison without/with the disturbance compensation.

Figure 12.

Rudder angle comparison without/with the disturbance compensation.

Conclusion

Aiming at the coordination problem of path following and rudder roll stabilization, the integrated constrained MPC control is proposed combined with the disturbance compensation.

According to rudder range limitation, the linear model is established, and the disturbance observer is adopted to estimate and compensate the lumped disturbance such that the capacity of nominal MPC can be enhanced to tackle time-varying disturbance. The non-minimum phase characteristic is analyzed, and a new minimum phase system approaching the original system is produced using the output redefinition method. For the sake of reduced calculation burden, the linear programming is utilized to tackle with the constrained optimization instead of quadratic programming. The simulations have been carried out to demonstrate the error convergence of the heading angle and path following with the constraints.

Although the proposed method offers an available realization of the coordination of the rudder roll stabilization and path following, the practical application still remains to be a big problem due to the complex sea conditions, and in the further works, the performance and stability of the proposed method will be tested in the experimental research.