Longitudinal motion control for high-speed trimaran based on computational fluid dynamics and predictive control

Abstract

This article investigates a model predictive control (MPC) with disturbance observer (DOB) for a trimaran longitudinal motion control. Firstly, to design the trimaran longitudinal motion stability controller, the mathematical model of the trimaran requires to be obtained. The hydrodynamic coefficients in the mathematical model are obtained by the computational fluid dynamics simulation. Secondly, a T-foil with the fixed attack angle is selected as an antipitching appendage. It is verified that the T-foil is effective in restraining the longitudinal motion of the trimaran through numerical simulation. Lastly, to enhance the ability of the T-foil for restraining the severe longitudinal motion, a controller based on the MPC method with DOB is designed to control the attack angle. The effect of the proposed algorithm is verified by theoretical analysis and simulation.

Introduction

Trimaran produces heave and pitch motions when sailing at high speed, which will adversely affect the navigation of the trimaran. The problem of how to restrain motion of a trimaran has caused great attention among engineers and technicians in different countries. For ship stabilization, the commonly used method is to install antipitching appendages on the bottom of the ship. Many appendages are used for ship stability like bilge keel, fin, T-foil, flap, and so on. ¹ In these appendages mentioned above, T-foil has the advantages of low resistance and small impact on the speed. In this study, T-foil is selected as the antipitching appendage for a high-speed trimaran.

The fixed T-foil installed on the trimaran has a certain antipitching capability. To improve stability performance, a control system is designed to change the attack angle of the appendage. Consequently, the choice of a reasonable control algorithm is critical to the control system. T-foil and flaps on the fast ferry were effective for suppressing the longitudinal motion and the worst vertical acceleration is proposed as the PD (Proportional, Differential) controller index. ^2,3 The PD controller has a simple principle, convenient application, and strong robustness, but the control accuracy is not very high. Multivariable classic control was used to decrease the motion sickness incidence of a high-speed ferry. ^4,5 The essence of the multivariable class control is also a PD controller (first-order filter and second-order filter). Multivariable class controller has a better effect on ship stability than PD controller, nonetheless, it still has the disadvantage of PD controller. Using a multivariable robust control algorithm for a ferry stabilization, the result proved the effectiveness of advanced algorithm and could be used in actual marine environment. ⁶ In addition, adaptive fuzzy control in tracking an uncertain marine vehicle with unknown disturbances was used. ⁷ The method has a good control effect, even so, it is easy to misoperate when the system breaks down. Sliding mode control was used to design a ship steering autopilot. ⁸ In Wang et al.’s study, ^{9 –13} a finite-time observer-based control scheme was performed, and external disturbance could be exactly attenuated with a good performance of the tracking control. In Zhang et al.’s study, ¹⁴ based on the model predictive controller, path-following of underactuated surface vessels was achieved. A translation rotation cascade control approach and a learning control framework were represented for unmanned surface vehicles. ^15,16 In Wang et al.’s study ¹⁷ a global asymptotic output tracking of nonlinear second-order systems with power integrators was proposed. These control algorithms have been widely used in the marine theory and engineering practice. In this article, considering the constraint of the T-foil, the method of model predictive control (MPC) combined with disturbance observer (DOB) is used for the trimaran longitudinal stability controller design, handling the constraint problem of the T-foil.

In this article, T-foil is selected as the antipitching appendage. As a result, the constraints of attack angle and the state of the trimaran must be considered. As is known, MPC can handle the constraint problem and convert the constraint conditions into quadratic programming while optimizing the corresponding performance index functions. ¹⁸ The control law contains the constraint condition and achieves the real-time control of the system under the premise of satisfying the constraint condition. Furthermore, the longitudinal motion control of a high-speed trimaran also expands the application of MPC. Due to the complexity of the marine environment, the trimaran sails in the wave can be encountered much unknown disturbance. Accordingly, a DOB is required to design. ^{19 –21} To tackle the problem of constraint and disturbance, the synthesis of MPC with DOB is proposed in this article.

In recent years, the computational fluid dynamics (CFD) simulation combined with experiment has been widely used. The numerical simulation results were consistent with the experiment. ^{22 –25} In this article, the mature CFD method is used to obtain the trimaran hydrodynamic coefficients.

The main contributions of this article are presented as follows:

A scale geometric model of the trimaran is built based on the trimaran profile and parameters. Using the CFD simulation to calculate the hydrodynamic coefficients of the trimaran, the antipitching appendage is designed, the lift force provided by the T-foil is calculated, and the slope of the lift force curve is determined by the polynomial fitting.

The model predictive controller with DOB is designed. The effective set method is used to constrain the amplitude of the attack angle of T-foil in actual application.

The structure of this article is as follows: the second section represents the model of the trimaran. In the third section, a T-foil is designed and approximately obtained the slope of the T-foil through the polynomial fitting. Then, an MPC with a DOB controller is designed and presented in the fourth section. The fifth section proves the effect of the proposed control algorithm through simulation.

Problem statement

For the control problem of the trimaran longitudinal motion, this article uses an active-controlled T-foil as the antipitching appendage. Firstly, a mathematical model of the trimaran is established. Consider the trimaran as an ordinary rigid body when the trimaran sails in marine. The forces and moments are produced as follows: buoyancy of the trimaran, the lift force, and moment generated by the trimaran and the active T-foil. The disturbance force and moment are generated by the wave, wind, and current. Through the dynamic theorem of the rigid body and the moment of momentum, the mathematical model for the longitudinal motion is shown as ²⁶

$\{\begin{array}{l}Z\left({\ddot{x}}_3,{\dot{x}}_3,x_3,{\ddot{x}}_5,{\dot{x}}_5,x_5,\right)=\\ m({\ddot{x}}_3+v{\dot{x}}_5)-F_T-2F_H-F_b\text{cos}(x_5)-mg\text{cos}(x_5)\\ M\left({\ddot{x}}_3,{\dot{x}}_3,x_3,{\ddot{x}}_5,{\dot{x}}_5,x_5,\right)=\\ I{\ddot{x}}_5+(l_T-x_G)F_T+2(x_H-x_G)F_H+(x_b-x_G)F_b\text{cos}(x_5)\end{array}$

1

where m is the mass of the trimaran, I is longitudinal inertia of the trimaran, x ₃ is heave motion, x ₅ is pitch motion, v is the velocity of the trimaran, F_T is the force generated by the T-foil, F_H is lift force of the trimaran, F_b is the trimaran buoyancy, F ₃ and M ₅ are the wave disturbance force and moment, $|l_T|,|x_G|,|x_H|,|x_b|$ are the T-foil gravity center, the trimaran gravity center trimaran, the distance of lift point, and the buoyancy point to the trimaran center, respectively.

The mathematical model in equation (1) describes the force situation of the trimaran during navigation. Typically, to design the controller, the mathematical model requires to be linearized as ²⁷

$\begin{array}{l}(m+a_{33}){\ddot{x}}_3+b_{33}{\dot{x}}_3+c_{33}{x}_3+a_{35}{\ddot{x}}_5+b_{35}{\dot{x}}_5+c_{35}{x}_5=F_T+F_3\\ (I+a_{55}){\ddot{x}}_5+b_{55}{\dot{x}}_5+c_{55}{x}_5+a_{53}{\ddot{x}}_3+b_{53}{\dot{x}}_3+c_{53}{x}_3=M_T+M_5\end{array}$

2

where $a_{ij}$ is the added mass, $b_{ij}$ is the damping coefficient, $c_{ij}$ is the restoring coefficient, and M_T is the moment of the T-foil.

Trimaran hydrodynamic

The trimaran parameters are listed in Table 1. For the trimaran studied in this article, the CFD simulation is used to obtain the hydrodynamic coefficients.

Table 1.

Parameters of the trimaran.

Name	Parameter
Main hull length (m)	3.343
Main hull beam (m)	0.240
Main hull waterline (m)	0.122
Main hull depth (m)	0.249
Side hull length (m)	1.157
Side hull beam (m)	0.064
Side hull waterline (m)	0.043
Side hull depth (m)	0.171
Side hull gas (m)	0.288
Mass (kg)	43.14
Center of gravity	(−1.622, 0, 0.1864)

Theory of ship computational fluid dynamics

Ship CFD is based on three equations of fluid mechanics (continuity equation, momentum equation, and energy equation). The Reynolds-averaged Navier–Stokes (N-S)[Please approve the abbreviation “N-S” identified for the expansion “Navier–Stokes” is appropriate.] method is the most widely acknowledged method for describing turbulent flows. Normally, ship CFD requires to solve the N-S equation. The complete N-S equation is presented as ²³

$\frac{\partial (\rho {\overline{u}}_i)}{\partial x_i}=0$

3

$\frac{\partial (\rho {\overline{u}}_i)}{\partial t}+\frac{\partial }{\partial x_j}(\rho {\overline{u}}_i{\overline{u}}_j+\rho \overline{{u^{\prime }}_i{u^{\prime }}_j})=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial {\overline{\tau }}_{ij}}{\partial x_j}$

4

where ${\overline{\tau }}_{ij}$ are the mean viscous stress tensor components, as shown in equation (5). p is the mean pressure, ${\overline{u}}_i$ is the averaged Cartesian components of the velocity vector, $\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}$ is the Reynolds stresses, ρ is the fluid density

${\overline{\tau }}_{ij}=\mu \left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$

5

where μ is the dynamic viscosity.

Simulation process

The geometry model of the trimaran is built by a 3-D modeling software. To obtain the hydrodynamic coefficients in equation (2), the calculation domain is required. Figure 1 shows the calculation domain, the domain includes six faces, and the type of faces is provided in Table 2.

Figure 1.

Calculation domain of the trimaran.

Table 2.

The boundary conditions.

Boundary condition	Type
Inlet	Velocity inlet
Outlet	Pressure outlet
Side	Symmetry
Bottom	Velocity inlet
Top	Velocity inlet

Then, the grid is necessary for hydrodynamic coefficients of the trimaran calculation. Figure 2 shows four grid conditions. The total number of the grid is 50–600 k corresponding to grid 1 to grid 4, and the calculation result shows that the number of grids has little effect on the results. In general, the denser the grid, the more accurate the calculation results. But, the denser the grid, the larger the calculation. To accelerate the simulation process, grid 3 is selected as the calculation grid in this article and shown in Figure 3(a) and (b).

Figure 2.

Four grid conditions: (a) grid 1, (b) grid 2, (c) grid 3, and (d) grid 4.

Figure 3.

(a) and (b) Distribution of the grid 3.

The working conditions are given in Table 3. The $k\text{-}\epsilon$ turbulence model is used to calculate the sailing state of the trimaran in different directions after building the calculation domain and meshing. The sailing states of the trimaran at the direction of 45°, 90°, 135°, and 180° are shown in Figure 4.

Table 3.

Working conditions.

Name	Value
Water density	$1.0\times {10}^3{\text{kg/m}}^3$
Speed of the trimaran	20 kn, 30 kn, 40 kn
Wave direction	45°, 90°, 135°, and 180°

Figure 4.

The trimaran sailing at different wave directions through CFD simulation: (a) 45°, (b) 90°, (c) 135°, and (d) 180°. CFD: computational fluid dynamics.

This article focuses on the trimaran sails in the heading wave. Therefore, the wave direction is set to 180°. From the CFD calculation above, the added mass and the damping coefficients are obtained and shown in Figures 5 and 6, respectively. The restoring coefficients are not changed with the speed and the encounter frequency, and Table 4 provides the restoring coefficients.

Figure 5.

(a)–(d) Added mass.

Figure 6.

(a)–(d) Damping coefficients.

Table 4.

The restoring coefficients.

Name	Parameter
$c_{33}(\text{N/m})$	6450.14
$c_{35}(\text{N})$	1080.1
$c_{53}(\text{N})$	1080.1
$c_{55}(\text{N/m})$	3920.7

Trimaran motion equation

Heave and pitch motions of the trimaran are mainly influenced by external factors, such as wind, wave, and current. For the trimaran longitudinal stability controller design, equation (2) is converted to a space state equation and expressed as

$\begin{array}{l}\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& m+a_{33}& 0& a_{35}\\ 0& 0& 1& 0\\ 0& a_{53}& 0& I+a_{55}\end{array}\right]\left[\begin{array}{l}{\dot{x}}_3\\ {\ddot{x}}_3\\ {\dot{x}}_5\\ {\ddot{x}}_5\end{array}\right]=\\ \left[\begin{array}{cccc}0& 1& 0& 0\\ -c_{33}& -b_{33}& -c_{35}& -b_{35}\\ 0& 0& 0& 1\\ -c_{53}& -b_{53}& -c_{55}& -b_{55}\end{array}\right]\left[\begin{array}{l}{x}_3\\ {\dot{x}}_3\\ x_5\\ {\dot{x}}_5\end{array}\right]+\left[\begin{array}{l}\text{0 0 }\\ \text{1 0}\\ \text{0 0}\\ \text{0 1}\end{array}\right]\left[\begin{array}{l}{F}_3\\ M_5\end{array}\right]\end{array}$

6

For convenience, the coefficient matrix is defined as

$\begin{array}{l}{\mathit{T}}_1=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& m+a_{33}& 0& a_{35}\\ 0& 0& 1& 0\\ 0& a_{53}& 0& I+a_{55}\end{array}\right],{\mathit{A}}_1=\left[\begin{array}{cccc}0& 1& 0& 0\\ -c_{33}& -b_{33}& -c_{35}& -b_{35}\\ 0& 0& 0& 1\\ -c_{53}& -b_{53}& -c_{55}& -b_{55}\end{array}\right]\\ {\mathit{g}}_{\text{1}}={\left[\begin{array}{cc}0& \begin{array}{ccc}1& 0& 0\end{array}\\ 0& \begin{array}{ccc}0& 0& 1\end{array}\end{array}\right]}^T\text{, }\mathit{A}=T_1{}^{-\text{1}}{\mathit{A}}_1,\mathit{g}=T_1{}^{-\text{1}}{\mathit{g}}_{\text{1}}\\ \mathit{C}=\left[\begin{array}{cc}1& \begin{array}{ccc}0& 0& 0\end{array}\\ 0& \begin{array}{ccc}0& 1& 0\end{array}\end{array}\right],\mathit{d}={\left[0F_3\text{ 0 }{M}_5\right]}^T\\ \mathit{x}={\left[\begin{array}{cccc}{x}_{\text{3}}& {\dot{x}}_{\text{3}}& x_{\text{5}}& {\dot{x}}_{\text{5}}\end{array}\right]}^T\end{array}$

Then, the trimaran motion equation without the control force is rewritten as

$\{\begin{array}{l}\dot{x}=\mathit{A}x+\mathit{g}d\\ \mathit{y}=\mathit{C}\mathit{x}\end{array}$

7

Appendage design and calculation

To reduce the heave and pitch motions of the trimaran, T-foil is selected as the appendage and installed on the trimaran. The airfoil type of the T-foil is NACA 0012; the size and 3-D view are shown in Figure 7.

Figure 7.

The parameters of T-foil.

Comparison of the bare trimaran and the trimaran with the fixed T-foil

The longitudinal distance of the installation is 180 mm from the pillar centerline to the bow, and the vertical distance is 26 mm from the horizontal foil profile chord line to the trimaran baseline. Figure 8(a) and (b) shows the T-foil installation position, while Figure 8(c) and (d) shows the comparison of the bare trimaran and the trimaran with T-foil in the CFD simulation environment. Through analyzing the heave and pitch of the trimaran with T-foil and without T-foil, the result shows that the stability effect of the trimaran with fixed T-foil is better than the bare trimaran. To improve the stability effect, a controller is designed to change the attack angle of the T-foil reasonably.

Figure 8.

(a)–(d) Trimaran with the fixed T-foil.

Lift force of the T-foil

To design the controller, the first step is to know the lift force provided by the T-foil. Also, the lift force of the T-foil is calculated by the CFD simulation, the step is the same as the CFD analysis in the second section. Figure 7 shows the parameters of the T-foil, and the T-foil grid is shown in Figure 9. Figure 9(a) shows the overall grid of the T-foil, Figure 9(b) shows the local grid, Figure 9(c) and (d) shows the streamline and vector graph, respectively. Then, the lift force of T-foil is obtained at different attack angles.

Figure 9.

(a)–(d) The grid, streamline, and vector figure of T-foil.

Figure 9 shows the attack angle of the T-foil is −5°. Using the same method, the lift force of different attack angles between −15° and 15° can be calculated, which is shown in Figure 10. Then, the lift force at different attack angles is obtained and presented in Table 5.

Figure 10.

The different attack angles of T-foil.

Table 5.

Lift force at different attack angles.

Attack angle (°)	Lift force (N)
−15°	−142.7520
−10°	−88.1742
−5°	−39.0657
5°	49.4885
10°	98.3871
15°	151.6214

The essence of antipitching is that the lift force and lift moment are provided by the T-foil counteract to the wave force and moment. For the controller design, the T-foil attack angle is controlled to change the lift force and moment. Therefore, the slope of the lift force requires to calculate to directly obtain the relationship between the lift force and the attack angles, which can be rewritten as ²⁸

$\{\begin{array}{l}{F}_{\text{T-foil}}=\frac{1}{2}\rho \mathit{A}{v}^2{C}_L\theta \hspace{0.17em}\text{=}\hspace{0.17em}{k}_T\theta \\ M_{\text{T-foil}}=F_{\text{T-foil}}{x}_T=k_T{l}_T\theta \end{array}$

8

where ρ is water density, A is the T-foil area, v is the trimaran velocity, C_L is the lift force coefficient, θ is the T-foil attack angle, x_T is the distance of the T-foil pillar center to the trimaran gravity center, and k_T is the slope of the lift force curve. k_T is obtained by fitting the lift force data at different attack angles. Figure 11 shows the data point and the fitting curve. Then, the slope k_T is solved and the value is 9.65.

Figure 11.

The lift force-fitting curve.

Disturbance observer and controller design

From equations (1) and (2), the heave and pitch motions of the trimaran are mainly caused by disturbance force and moment of waves. To improve the dynamic performance and robustness of the control system, the interference suppression is necessary. As a result, it is necessary to use a DOB to estimate the uncertainty disturbance in real time. The model of the disturbance control system can be approximately equivalent to the nominal model. Then, a model predictive controller based on DOB is designed to improve the antidisturbance performance of the control system.

Disturbance observer design

With the control force, equation (7) can be rewritten as the state space and shown as

$\{\begin{array}{l}\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}u+\mathit{g}\mathit{d}\\ \mathit{y}=\mathit{C}\mathit{x}\end{array}$

9

where

$\begin{array}{l}\mathit{B}={\mathit{T}}_1{}^{-1}{\mathit{B}}_2\\ B_2={\left[\begin{array}{l}0k_T\text{ 0 0}\\ \text{0 0 0 }{k}_T{l}_T\end{array}\right]}^T\\ k_T=\frac{1}{2}\rho \mathit{A}{v}^2{C}_L\end{array}$

After performing several simple mathematical manipulations, equation (10) is obtained as

$\mathit{g}\mathit{d}=\dot{\mathit{x}}-\mathit{A}(x)-\mathit{B}u$

10

The initial DOB is designed as

$\dot{\widehat{\mathit{d}}}=-\mathit{l}(x)\mathit{g}\widehat{\mathit{d}}+\mathit{l}(x)(\dot{x}-\mathit{A}(x)-Bu)=\mathit{l}(x)\mathit{g}(\mathit{d}-\widehat{\mathit{d}})$

11

where $l(x)=\partial p(x)/\partial x$ is the gain matrix of the DOB and $\widehat{d}$ is the estimator of $d$ .

However, the designed observer is not implemented due to that the disturbance is unknown. The auxiliary variable z is induced, and the DOB is designed as

$\{\begin{array}{l}\widehat{\mathit{d}}=z+\mathit{p}(x)\\ \dot{z}=-\mathit{l}(x)z-l(x)(\mathit{p}(x)+\mathit{A}(x)+\mathit{B}u)\end{array}$

12

where $\mathit{p}(\mathit{x})=[l_1{l}_2{l}_2]{[x_1{x}_2{x}_3]}^{\text{T}}$ .

The MPC with DOB is shown in Figure 12. The DOB is used to eliminate external disturbances $d(t)$ and the adverse effects of the model mismatch system. Not only it is easy to be implemented in physical form but also it needs to suppress external disturbance for ensuring robust stability of the control system. Therefore, the primary goal of the designed DOB is used to guarantee internal stability. The Lyapunov second method is used to prove the convergence stability of the designed DOB in equation (12).

Figure 12.

MPC with disturbance observer. MPC: model predictive control.

The estimation error of the DOB is defined as

$e_1=d-\widehat{d}$

13

${\dot{e}}_1(t)=\dot{d}-\dot{\widehat{d}}=l(x)\mathit{g}(\widehat{d}-d)=-l(x)\mathit{g}{e}_1(t)$

14

Then, the dynamic equation of the designed DOB is

${\dot{e}}_1(t)+l(x)\mathit{g}(x)e_1(t)=0$

15

The Lyapunov function is defined as

$V(e_1)=\frac{1}{2}{e}_1^T{e}_1$

16

where $V(e_1)$ is positive, and $V(0)=0$ .

Derivative $V(e_1)$ , the result is shown as

$\dot{V}(e_1)=e_1^T(t){\dot{e}}_1(t)=-e_1^T(t)l(x)\mathit{g}{e}_1(t)=-l(x)\mathit{g}{e}_1{}^2(t)$

17

There exists a positive number $\alpha >0$ to guarantee $\dot{V}(e_1)\le -\alpha e_1{}^2(t)<0$ , $\dot{V}(e_1)$ is negative, and $\dot{V}(0)=0$ . Based on the Lyapunov stability theory, the estimation error of the DOB is asymptotically stable over a wide range $e_1=0$ and $\underset{t\to \infty }{\text{lim}}{e}_1(t)=0$ . Therefore, a suitable gain matrix $l(x)$ can be selected for the DOB to guarantee $e_1(t)$ stable convergence and the disturbance estimation $\widehat{d}$ approach to d. Then, the designed DOB in equation (12) can be achieved for the purpose of estimating disturbance online.

Controller design

The control system of MPC with DOB is shown in Figure 13. Using a zero-order holder to discretize equation (9), the result is given by

$\{\begin{array}{l}\mathit{x}(k+1)={\mathit{A}}_k\mathit{x}(k)+{\mathit{B}}_{\mathit{k}}u(k)+{\mathit{g}}_{\mathit{k}}\mathit{d}(k)\\ \mathit{y}(k)={\mathit{C}}_{\mathit{k}}\mathit{x}(k)\end{array}$

18

Figure 13.

Control system of the trimaran.

Equation (18) is regarded as a predictive model, and the following prediction states are obtained by the iteration

$\{\begin{array}{l}\mathit{x}(k+1|k)={\mathit{A}}_{\mathit{k}}\mathit{x}(k)+{\mathit{B}}_{\mathit{k}}u(k)+{\mathit{g}}_{k}d(k)\\ \mathit{x}(k+2|k)={\mathit{A}}_{\mathit{k}}{}^{2}x(k)+{\mathit{A}}_{\mathit{k}}{\mathit{B}}_{\mathit{k}}u(k)+{\mathit{B}}_{\mathit{k}}u(k+1)\\ +{\mathit{A}}_{\mathit{k}}{\mathit{g}}_{k}d(k)+{\mathit{g}}_{k}d(k+1)\\ ⋮\\ x(k+N_p|k)={\mathit{A}}_k{}^{N_P}x(k)+{\mathit{A}}_k{}^{N_P-1}{\mathit{B}}_{\mathit{k}}u(k)+\cdots +\\ {\mathit{A}}_{\mathit{k}}{}^{N_P-N_C}{B}_{k}u(k+N_c-1)+\\ {\mathit{A}}_{\mathit{k}}{}^{N_P-1}{\mathit{g}}_{k}d(k)+\cdots +\\ {\mathit{A}}_{\mathit{k}}{}^{N_P-N_C}{\mathit{g}}_{k}d(k+N_c-1)\end{array}$

For convenience, some matrices are defined as

$\begin{array}{l}\mathit{Y}={[\begin{array}{cccc}{\mathit{y}}^T(k+1)& {\mathit{y}}^T(k+2)& \cdots & {\mathit{y}}^T(k+N_P)\end{array}]}^{\text{T}}\\ \mathit{U}={[\begin{array}{cccc}{u}^T(k)& u^T(k+1)& \cdots & u^T(k+N_P-1)\end{array}]}^{\text{T}}\\ \mathit{F}={[\begin{array}{cccc}{({\mathit{C}}_k{\mathit{A}}_k)}^T& {({\mathit{C}}_k{\mathit{A}}_k{}^2)}^T& \cdots & {({\mathit{C}}_k{\mathit{A}}_k{}^{N_P})}^T\end{array}]}^{\text{T}}\\ \mathit{D}={[\begin{array}{cccc}{d}^T(k)& d^T(k+1)& \cdots & d^T(k+N_P-1)\end{array}]}^{\text{T}}\end{array}$

$\mathit{\Phi }=\left[\begin{array}{cccc}{\mathit{C}}_{\mathit{k}}{\mathit{B}}_{\mathit{k}}& 0& \cdots & 0\\ {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{\mathit{B}}_{\mathit{k}}& {\mathit{C}}_{\mathit{k}}{B}_k& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C_k{\mathit{A}}_{\mathit{k}}{}^{N_P-1}{\mathit{B}}_{\mathit{k}}& {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{}^{N_P-2}{\mathit{B}}_{\mathit{k}}& \cdots & {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{}^{N_P-N_C}{\mathit{B}}_{\mathit{k}}\end{array}\right]$

$\mathit{T}=\left[\begin{array}{cccc}{\mathit{C}}_{\mathit{k}}{\mathit{g}}_k& 0& \cdots & 0\\ {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{\mathit{g}}_k& {\mathit{C}}_{\mathit{k}}{\mathit{g}}_k& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{}^{N_P-1}{\mathit{g}}_k& {\mathit{C}}_{\mathit{k}}{\mathit{A}}_{\mathit{k}}{}^{N_P-2}{\mathit{g}}_k& \cdots & {\mathit{C}}_{\mathit{k}}{\mathit{g}}_k\end{array}\right]$

where N_P is the predictive step and N_C is the control step. Combining equation (18) and substituting $\mathit{Y},\mathit{U},\mathit{F},\mathit{\Phi }$ into discrete state equations as

$\mathit{Y}=\mathit{F}\mathit{x}(k)+\mathit{\Phi }\mathit{U}+\mathit{T}\mathit{D}$

19

The disturbance matrix items $\mathit{D}{(k)}^T\mathit{H}\mathit{D}(k)$ add to the performance index to obtain a new time-domain optimization indicator function of the trimaran control system

$\overline{J}={(Y_d-Y)}^T\overline{\mathit{Q}}(Y_d-Y)+{\mathit{U}}^T\overline{\mathit{R}}\mathit{U}+{\mathit{D}}^T\mathit{H}\mathit{D}$

20

where $\overline{\mathit{Q}}>0$ , $\overline{\mathit{R}}>0$ , and $\mathit{H}>0$ are positive symmetric weighted matrices.

The necessary condition for minimizing the objective function is $\partial \overline{J}/\partial \mathit{U}=0$ . Based on the unconstrained disturbance optimization problem, $\mathit{U}$ can be obtained as

$\mathit{U}={({\mathit{\Phi }}^T\overline{\mathit{Q}}\mathit{\Phi }+\overline{\mathit{R}})}^{-1}[{\Phi }^T\overline{Q}(Y_d-\mathit{F}\mathit{x}(k)-\mathit{T}\mathit{D}]$

21

The error of control system $e(t)$ is defined as

$e(t)=y_d(t)-y(t)$

22

The Lyapunov function is selected as

$V=\frac{1}{2}{e}^T\overline{\mathit{Q}}e+\frac{1}{2}{e}_1{}^2$

23

where e ₁ is the estimated error of the designed DOB. Due to the matrix, $\overline{\mathit{Q}}$ is a positive definite symmetric constant matrix, and the Lyapunov function is satisfied with $V\ge 0$ . According to the principle of predictive control, in the limited prediction time domain, the output of the system $y(t)$ is to track the given composite output $y_d(t)\text{=}0$ . Through the predictive control laws in equation (21), the rolling optimization indicator in equation (20) can be approximated as

$J\approx e^T\overline{Q}e$

24

Derivative (23) and combining (24), equation (25) can be obtained as

$\dot{V}=\frac{1}{2}\frac{\partial J}{\partial U}+e_1^T{\dot{e}}_1$

25

Then, equation (17) is shown as

$\dot{V}(e_1)=-l(x){\mathit{g}}_2(x)e_1{}^2\le 0$

26

It can be seen from the prediction objective function in equation (20) that to achieve the minimum performance index, it must be satisfied with

$\frac{\partial J}{\partial U}=0$

27

From equations (25) to (27), the result is obtained as

$\dot{V}\le 0$

28

Then, the closed-loop system error asymptotically converges and the system is consistently stability.

Constraint handling

The range of T-foil attack angle is limited in actual manipulation. There exists a significant problem in the optimal control law obtained from the previous section. When the control quantity is higher than the constraints of actuator in long term, it damages the physical actuator. Consequently, it is necessary to limit the range of the attack angle of the T-foil during the control system design. When the predictive control step is N_p , at the time k, the constraint of the attack angle control amount in the future step j can be expressed as

$u_{\text{min}}\le u(k+j-1)\le u_{\text{max}},j=\begin{array}{ccc}1,& \cdots ,& N_P\end{array},\cdots$

And, the vector forms are shown as

$WU\le r$

29

where $W={\left[\begin{array}{cc}I& -I\end{array}\right]}^T,r={\left[\begin{array}{cc}{U}_{\text{max}}{}^T& -U_{\text{min}}{}^T\end{array}\right]}^T$ , and I is the $N_P\times N_P$ dimension matrix

$\begin{array}{l}U={\left[u\left(k\right),u\left(k+1\right),\cdots ,u\left(k+N_p-1\right)\right]}^T\\ U_{\text{max}}={\left[\begin{array}{ccc}{u}_{\text{max}}& u_{\text{max}}& \begin{array}{cc}\cdots & u_{\text{max}}\end{array}\end{array}\right]}^T{}_{{}_{N_P\times 1}}\\ U_{\text{min}}={\left[\begin{array}{ccc}{u}_{\text{min}}& u_{\text{min}}& \begin{array}{cc}\cdots & u_{\text{min}}\end{array}\end{array}\right]}^T{}_{N_P\times 1}\end{array}$

Equations (20) and (29) are organized to obtain the control problem with constraint conditions and written as

$\begin{array}{l}\text{min} J\\ \text{s.t.}\{\begin{array}{l}J={(Y_d-Y)}^T\overline{Q}(Y_d-Y)+U^T\overline{R}U\text{+}{D}^{T}HD\\ WU\le r\end{array}\end{array}$

30

With $Y_{\text{d}}\text{=}0,Y=Fx(k)+\Phi U+TD$ , the standard quadratic form can be obtained as

$\begin{array}{l}\text{min} J\\ \text{s.t.}\{\begin{array}{l}J=\frac{1}{2}{U}^{T}MU+c^{T}U+D^{T}HD\\ WU\le r\end{array}\end{array}$

31

where $M=2{\Phi }^T\overline{Q}\Phi +2\overline{R}$ and $c=2{\Phi }^T\overline{Q}Fx\left(k\right)$ .

From equation (31), the original problem is converted to a quadratic programming problem. The Lagrange multiplier method and the effective set method in predictive control are used to deal with the constraint conditions. In this way, the inequality problem is converted to an equality constraint by adding a Lagrange multiplier vector, and the Lagrange function is given by

$L_x={(Y_d-Y)}^T\overline{Q}(Y_d-Y)+U^T\overline{R}U+D^{T}HD+{\gamma }^T(Wu-r)$

32

where $\gamma$ is the multiplier vector.

Then, L_x is partial derivative to U and γ, the matrix form is obtained as

$\left[\begin{array}{cc}\mathit{M}& {\mathit{W}}^T\\ \mathit{W}& 0\end{array}\right]\left[\begin{array}{c}U\\ \gamma \end{array}\right]=\left[\begin{array}{c}-c^T\\ r\end{array}\right]$

33

The effective set method is used in the predictive control, the idea is: For any feasible solution $U^{*}$ , it must be satisfied with all constraints. The constraints that satisfied the equation are called effective constraints, the set of all effective constraints is called an active set or positive set with the remark $\Omega \left(U\right)$ . Then, the problem of solving the specific quadratic programming corresponding to the effective set is expressed as

$\begin{array}{l}\text{min} J\\ \text{s.t.}\{\begin{array}{l}J=\frac{1}{2}{U}^{T}MU+c^{T}U+D^{T}HD\\ w_{i}U=r_i\hspace{0.17em}\hspace{0.17em},i\in \Omega (U)\end{array}\end{array}$

34

where w_i is the i’th row vector of the constraint set $W$ and r_i is the i’th element of r.

The feasible solutions $U^{*}$ and ${\gamma }^{*}$ are obtained by equation (29), and the nonsingular matrix $W$ is the set of all valid constraints at this time. The constraint set $W$ must be a nonsingular matrix. The feasible solution can be expressed as

$\left[\begin{array}{c}{U}^{*}\\ {\gamma }^{*}\end{array}\right]={\left[\begin{array}{cc}\mathit{M}& {\mathit{W}}^T\\ \mathit{W}& 0\end{array}\right]}^{-1}\left[\begin{array}{c}-c^T\\ r\end{array}\right]$

35

To facilitate the iterations for predictive control, the feasible solution $U^{*}$ is changed to $U^{*}+\delta$ and substituted it with equation (34). Then, equation (34) is rewritten as

$\begin{array}{l}\text{min} J\\ \text{s.t.}\{\begin{array}{l}J=\frac{1}{2}{(U^{*}+\delta )}^{T}M(U^{*}+\delta )+c^T(U^{*}+\delta )+D^{T}HD\\ w_i\delta =0\hspace{0.17em},i\in {\Omega }^{*}(U^{*})\end{array}\end{array}$

36

At each iteration, $U^{*}$ is given and equation (31) is rewritten as

$\begin{array}{l}\text{min}J\\ \text{s.t.}\{\begin{array}{l}{J}^{\prime }=\frac{1}{2}{\delta }^{T}M\delta +(U^{*}M+c^T)\delta +D^{T}HD\\ w_i\delta =0\hspace{0.17em},i\in {\Omega }^{*}(U^{*})\end{array}\end{array}$

37

The Lagrange multiplier method is also derivative to δ and γ and the feasible solutions ${\delta }^{*}$ and ${\gamma }^{*}$ are obtained by matrix division. For different value assignments of ${\delta }^{*}$ and ${\gamma }^{*}$ , the following cases are discussed.

Remark 1

When ${\delta }^{*}=0$ and the whole elements ${\gamma }_i\ge 0$ in the vector ${\gamma }^{*}$ , $U^{*}$ is the final solution to the quadratic problem. ${\gamma }_k{}^{*}={\gamma }_k$ is taken, the control sequence under constraint conditions is formulated as

$U={({\Phi }^T\overline{Q}\Phi +\overline{R})}^{-1}[{\Phi }^T\overline{Q}(Y_d-Fx(k))-\frac{1}{2}{W}^T{\gamma }^{*}]$

38

To improve the robustness of predictive control, a rolling horizon optimization strategy is adopted. The first item of the optimal control sequence is applied to the controlled object, such that the control law with the constraint condition is shown as

$\begin{array}{l}{\delta }_{mpc}(k)=[I\text{ 0 0 }\dots \text{0]}U\\ \text{=}[I\text{ 0 0 }\dots \text{0]}{({\Phi }^T\overline{Q}\Phi +\overline{R})}^{-1}\times \\ [{\Phi }^T\overline{Q}(Y_d-Fx(k)-TD)-\frac{1}{2}{W}^T{\gamma }^{*}]\end{array}$

39

where ${\gamma }_k{}^{*}={\left[\begin{array}{ccc}{\gamma }_1{}^{*}& {\gamma }_2{}^{*}& \cdots \end{array}\right]}^T$ is the Lagrange multiplier vector.

Remark 2

When ${\delta }^{*}=0$ and some element ${\gamma }_j<0$ in the vector ${\gamma }^{*}$ , it represents that the corresponding constraints ${\gamma }_j<0$ are not satisfied and $U^{*}$ is not an optimal solution. Then, the corresponding constraint ${\gamma }_j<0$ is removed from the effective set. If multiple constraints are not satisfied, then, the smallest one ${\gamma }_j$ is removed and resolved equation (37).

Remark 3

When ${\delta }^{*}\ne 0$ and $U^{*}+\delta$ are a feasible solution, $U^{*}=U^{*}+\delta$ is taken and resolved equation (37).

Remark 4

When ${\delta }^{*}\ne 0$ and $U^{*}+{\delta }^{*}$ are not a feasible solution, it represents that the step size ${\delta }^{*}$ exceeds the specified constraint limitation. It requires to consider whether the step size ${\alpha }^{*}=\text{max}(\alpha >0|U^{*}+\alpha {\delta }^{*})$ is feasible, which starts from $U^{*}$ and alone ${\delta }^{*}$ to find the optimal value that satisfied the constraint conditions. During this process, there must be at least one constraint $l\notin {\Omega }^{*}(U^{*})$ satisfies $w_l(U^{*}+\alpha \delta )=r_l$ , then the effective set is modified to ${\Omega }^{*}(U^{*})\cup \left\{l\right\}$ , and resolved equation (37).

Summarizing the conditions discussed above: If a step size is ${\delta }^{*}\ne 0$ , it indicates that the solution $U^{*}$ is not an optimal solution, which can be further improved. If ${\delta }^{*}=0$ , it is impossible to improve under the current constraints. In addition, for the parameter ${\gamma }^{*}$ , when some of these elements are ${\gamma }_j<0$ , it indicates that the corresponding constraints limit the further improvement of the feasible solution. If all ${\gamma }_j>0$ are satisfied, it shows that the obtained solution is the optimal solution with current constraints.

Simulation analysis

In the simulation, the random wave spectrum is used to simulate the irregular wave. ²⁹ The sea state is 4 with the wave height range 1.25–2.5 m. As the trimaran model is a 1:23 model, the wave range is changed to 0.054–0.110 m and the heading wave is selected. Considering the trimaran under the condition of 40 knots, the above matrix parameters in “Trimaran motion equation” section, which can be calculated are

$\begin{array}{l}{a}_{33}=165.1808, a_{35}=33.34065, a_{53}=\text{35.55758}, \\ a_{55}=\text{963.1973,}{b}_{33}=\text{0.6493438}, b_{35}=\text{883.4414}, \\ b_{53}=-\text{257.0464}, b_{33}=\text{113.0062, }{l}_T=0.72, k_T=9.65\end{array}$

Remark: For the control weigh matrix, the larger the error weight matrix is, the shorter the response time of each state and the greater the overshoot are. For the error weight matrix, the larger the control weight matrix is, the slower the response of the system and the smoother transition between every state are. For a trimaran with a relatively high-speed requirement, a larger error weight matrix should be selected so that each state of the trimaran can respond quickly, and in this article, the value $Q\text{ = }200$ and $R\text{ = }1$ is chosen. After determining the two-weight matrices, the larger the prediction step, the faster the system’s state response, so that the faster adjustments can be made, but larger oscillations may occur. The small prediction step can also lead to system instability. With the compromise, the predictive control step is chosen $N_p=3$ .

Using the determined weight matrix and prediction step size in the controller, a PID (Proportional, Integral, Differential) controller is used as a comparison. The simulation result is shown in Figure 14. Figure 14(a) and (b) shows the heave and pitch motion comparison, while Figure 14(c) and (d) shows the attack angle of the PID controller and the controller of MPC with DOB.

Figure 14.

The control method comparison. (a) and (b) The heave comparison, (c) the attack angle changed in PID controller, (d) the attack angle changed in MPC + DOB controller. MPC: model predictive control; DOB: disturbance observer.

Conclusion

To reduce the heave and pitch motions of the trimaran, a T-foil is selected as the appendage and installed on the bottom of the trimaran. Then, the controller of MPC with DOB is designed to change the T-foil attack angle and perform the algorithm through simulation. As a comparison, the PID controller is performed at the same time. Through analyzing the result, the heave motion decreased by 56.73% and the pitch motion decreased by 42.15% using the controller of MPC with DOB, while the heave motion decreased by 32.75% and the pitch motion decreased by 34.26% using the PID controller. It is seen from Figure 14 that the MPC with DOB can greatly suppress the trimaran longitudinal motion.