Artificial Neural Networks in Coordinated Control of Multiple Hovercrafts with Unmodeled Terms

Abstract

In this paper, the problem of coordinated control of multiple hovercrafts is addressed. For a single hovercraft, by using the backstepping technique, a nonlinear controller is proposed, where Radial Basis Function Neural Networks (RBFNNs) are adopted to approximate unmodeled terms. Despite the application of RBFNNs, integral terms are introduced, improving the robustness of controller. As a result, global uniformly ultimate boundedness is achieved. Regarding the communication topology, two different directed graphs are chosen under the assumption that there are no delays when they communicate with each other. In order to testify the performance of the proposed strategy, simulation results are presented, showing that vehicles can move forward in a specific formation pattern and RBFNNs are able to approximate unmodeled terms.

1. Introduction

In recent years, Wireless Sensor Networks (WSNs) have attracted growing interests from researchers, because they have merits, compared with traditional networking solutions, such as reliability, flexibility, and an ease of deployment, that enable their use in a wide range of varied application scenarios [1]. They can be applied to track moving objects, to monitor special areas so as to trigger alarm systems when some dangerous signals are detected, etc. As the eyes and ears of the IoT, WSNs can work as bridges to build connections between the real-world and the digital-world. In light of this promising application scenario, this paper mainly focuses on a case study of mobile WSNs, where a group of hovercrafts equipped with specific sensors are chosen as test platforms. The objective is to enable them to move around and interact with the physical environment [2] and thus execute a mission of mapping, searching, and monitoring in a specific area.

Coordinated control of a fleet of hovercrafts is challenging, especially when we take into account their complex dynamic models. Until now, for a single surface vehicle, many research results have been reported. For example, a linear fuzzy-PID controller was proposed in [3]. Compared with the ordinary PID controller, the proposed controller therein performs better in term of improving settling time and reducing overshoot of the control signal. However, their works just consider the kinematic models of the vehicle without considering the dynamic models, which is not realistic in real operation scenarios. Another weakness of the linear controller is that it usually achieves local stability, e.g., [4], where velocity and position controllers were developed based on a linearized system, which is controllable only when the angular velocity is nonzero. Considering the limitations of linear controllers, in [5,6], nonlinear controllers for underactuated ships were designed, and global asymptotic stability is achieved. In [7], a nonlinear Lyapunov-based tracking controller was presented, and it was able to exponentially stabilize the position tracking error to a neighborhood of the origin that can be made arbitrarily small. A method of incorporating multiobjective controller selection into a closed-loop control system was presented in [8], where the authors designed three controllers so as to capture three “behaviors” representative of typical maneuvers that would be performed in a port environment. However, none of the works mentioned above consider disturbances and unmodeled terms of the vehicle. In order to ensure that the vehicle is robust to external disturbances, in [9], two controllers with application to a surface vehicle (named Qboat) and a hovercraft were proposed, and the authors designed disturbance estimators to estimate external constant disturbances. The disadvantage of this control strategy is that they did not consider unmodeled terms involved in the dynamic model of the vehicle. Considering this constraint, an estimator was developed in [10], where a fuzzy system was used to approximate unknown kinetics. A fault tolerant tracking controller was designed in [11] for a surface vessel. In addition, a self-constructing adaptive robust fuzzy neural control scheme for tracking surface vessels was proposed in [12], where simulation results were shown to testify the efficiency of the proposed method therein.

With respect to coordinated control strategy for multiple vehicles, many authors have presented their own approaches. In [13], a cooperative path following methodology was proposed under the assumption that the communication among a group of fully-actuated surface vehicles is undirected and continuous. A coordinated path following with a switching communication topology was designed in [14], while a null-space-based behavioral control technique was proposed in [15,16]. In [17,18,19], a leader–follower control strategy was presented. In [20], an adaptive coordinate tracking control problem for a fleet of nonholonomic chained systems was discussed under the assumption that the desired trajectory is available only to part of the neighbors. The reader is also referred to [21] for more results about multi-vehicle control approaches.

Inspired and motived by those works mentioned above, in this paper, we first develop a controller that is able to drive a single hovercraft to the neighborhood of a desired smooth path, where a Radial Basis Function Neural Network (RBFNN) is applied to approximate unmodeled dynamic terms of the vehicle while integral error terms are introduced, thus improving the robust performance of the controller. It is relevant to point out that all elements of the estimation weight matrix are always bounded through the use of a smooth projection function. We also derive a consensus strategy to make sure the desire paths progress in a specific formation. In order to validate the effectiveness of the proposed strategy, simulation results are presented.

The rest of the paper is organized as follows: Section 2 presents robot modeling, graph theory, RBFNNs, and coordinated control problem. A single controller is proposed in Section 3, while Section 4 devises a consensus strategy. Simulation results are given in Section 5 to validate the performance of the proposed approach herein. At last, Section 6 summarizes our work and describes the future work.

2. Problem Formulation

2.1. Vehicle Modeling

We first define a global coordinate frame $\left\{U\right\}$ and a body frame $\left\{B\right\}$ as shown in Figure 1. The kinematic equations of the vehicle are written as

$$\begin{array}{cc}\hfill \dot{\mathbf{p}}& =\mathbf{R}\mathbf{v}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \dot{\psi }& =\omega \hfill \end{array}$$

(2)

where $\dot{\mathbf{p}}={[x,y]}^T$ denotes the coordinates of its center of mass, $\mathbf{v}={[u,v]}^T$ represents linear velocity, $\psi$ is the orientation of the vehicle, and its angular velocity is represented by $\omega$. Moreover, the rotation matrix $\mathbf{R}\left(\psi \right)$ is given by

$$\mathbf{R}\left(\psi \right)=\left[\begin{array}{cc}\cos \left(\psi \right)& -\sin \left(\psi \right)\\ \sin \left(\psi \right)& \cos \left(\psi \right)\end{array}\right].$$

(3)

Its dynamic equations are

$$\begin{array}{cc}\hfill \dot{\mathbf{v}}& =-\mathbf{S}\left(\omega \right)\mathbf{v}+m^{-1}T{\mathbf{n}}_1+{\Delta }_v\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \dot{\omega }& =-J^{-1}\omega +J^{-1}\tau +{\Delta }_{\omega }\hfill \end{array}$$

(5)

where $\mathbf{S}\left(\omega \right)$ is a skew symmetric matrix, given by

$$\begin{array}{c}\hfill \mathbf{S}\left(\omega \right)=\left[\begin{array}{cc}0& -\omega \\ \omega & 0\end{array}\right].\end{array}$$

(6)

${\mathbf{n}}_1={[1,0]}^T$, m and J denote the car’s mass and rotational inertial, respectively. The force used to make the car move forward is denoted by T, and $\tau$ represents the torque that can steer the vehicle. Unmodeled dynamic terms are represented by ${\Delta }_v$ and ${\Delta }_{\omega }$. For more details about modeling surface vehicles, the reader is referred to [22].

2.2. Graph Theory

In this paper, $\mathcal{G}=\mathcal{G}(\mathcal{V},\mathcal{E})$ denotes a directed graph that can be used to model the interaction communication topology among mobile robots. The graph $\mathcal{G}$ consists of a finite set $\mathcal{V}=\{1,2,\dots ,n\}$ of n vehicles and a finite set $\mathcal{E}$ of m pairs of vertices ${\mathcal{V}}_{ij}=\{i,j\}\in \mathcal{E}$. If ${\mathcal{V}}_{ij}$ belongs to $\mathcal{E}$, then i and j are said to be adjacent. A graph from i to j is a sequence of distinct vertices starting with i and ending with j such that consecutive vertices are adjacent. In this case, ${\mathcal{V}}_{ij}$ also represents a directional communication link from agent i to agent j. The adjacency matrix of the graph $\mathcal{G}$ is denoted by $\mathcal{A}=\left[a_{ij}\right]\in {\mathbb{R}}^{n\times n}$, which is a square matrix where $a_{ij}$ equals to one if $\{j,i\}\in \mathcal{E}$ and zero otherwise. Moreover, the Laplacian matrix $\mathcal{L}$ is defined as $\mathcal{L}=\mathcal{D}-\mathcal{A}$, where the degree matrix $\mathcal{D}=\left[d_{ij}\right]\in {\mathbb{R}}^{n\times n}$ of the graph $\mathcal{G}$ is a diagonal matrix and $d_{ij}$ equals the number of adjacent vertices of vertex i.

2.3. Radial Basis Function Neural Networks

Radial Basis Function Neural Networks (RBFNNs) can be used to approximate the unmodeled nonlinear dynamic terms due to their universal approximation capability [23]. For any unknown smooth function $f\left(\mathbf{x}\right):{\mathbb{R}}^n\to {\mathbb{R}}^m$ can be approximated by RBFNNs in the following form, given by

$$\begin{array}{c}\hfill \widehat{f}\left(\mathbf{x}\right)={\mathbf{W}}^T\mathit{\sigma }\left(\mathbf{x}\right)\end{array}$$

(7)

where $\mathbf{x}\in \Omega \subset {\mathbb{R}}^n$, $\Omega$ is a compact set. The adjustable weight matrix with n neurons is denoted by $\mathbf{W}\in {\mathbb{R}}^{n\times m}$ under the assumption that it is a bounded matrix, that is

$$\begin{array}{c}\hfill \mathbf{W}\le {\mathbf{W}}_{\max }.\end{array}$$

(8)

It is important to point out that here when we say matrix $\mathbf{x}\in {\mathbb{R}}^{m\times n}$ is smaller than or equal to ${\mathbf{x}}_{\max }\in {\mathbb{R}}^{m\times n}$, we mean all elements of $\mathbf{x}$ are smaller than or equal to their corresponding elements of ${\mathbf{x}}_{\max }\in {\mathbb{R}}^{m\times n}$. Moreover, $\mathit{\sigma }\left(\mathbf{x}\right)$ is the basis function vector and ${\sigma }_i\left(x_i\right)=\exp (-{(x_i-{\mu }_i)}^T(x_i-{\mu }_i)/c_i^2),i=1,2,\dots ,n$ denotes its component, ${\mu }_i$ is the center of the receptive field, and $c_i$ represents the width of the Gaussian function. Moreover, it is relevant to point out that, in order to achieve better approximate results, we should make the neuron number n large enough and choose the parameters properly. Going back to the smooth function $f\left(\mathbf{x}\right)$ mentioned above, there is an ideal weight ${\mathbf{W}}_d$ such that

$$\begin{array}{c}\hfill f\left(\mathbf{x}\right)={\mathbf{W}}_d^T\mathit{\sigma }\left(\mathbf{x}\right)+\mathit{ϵ}\left(\mathbf{x}\right)\end{array}$$

(9)

where $\mathit{ϵ}\left(\mathbf{x}\right)$ denotes the approximation error and satisfies $\left|\right|\mathit{ϵ}\left(\mathbf{x}\right)\left|\right|\le {ϵ}_{\max }$, where ${ϵ}_{\max }$ is a positive number. It is noted that ${\mathbf{W}}_d$ is an “artificial” quantity for the purpose of mathematic analysis, in the process of controller design, we need to estimate it [24]. A simple RBFNNs is given by Figure 2.

2.4. Problem Statement

Now we can state our problem: Through designing a controller for each robot and proposing consensus strategies for their corresponding desired paths, we want a group of mobile robots to move forward in a specific formation pattern; that is,

(1)

for an individual vehicle, $\left|\right|\mathbf{p}-{\mathbf{p}}_d\left|\right|\to \delta$, where $\delta$ is an arbitrarily small constant value;

(2)

for a group of n desired paths, ${\gamma }_i-{\gamma }_j\to 0$ and ${\dot{\gamma }}_i-{\dot{\gamma }}_d\to 0$, where agent j is the neighbor of agent i, ${\dot{\gamma }}_d$ represents the desired value of ${\dot{\gamma }}_i$, which is a known value.

3. Controller Design

Following the works of [7] and [9], we define the position error in the body frame as

$$\begin{array}{c}\hfill {\mathbf{e}}_1={\mathbf{R}}^T(\mathbf{p}-{\mathbf{p}}_d)\end{array}$$

(10)

where $\mathbf{R}=\mathbf{R}\left(\psi \right)$ denotes the rotation matrix. The time derivative of ${\mathbf{e}}_1$ yields

$$\begin{array}{c}\hfill {\dot{\mathbf{e}}}_1=-\mathbf{S}\left(\omega \right){\mathbf{e}}_1+\mathbf{v}-{\mathbf{R}}^T{\dot{\mathbf{p}}}_d.\end{array}$$

(11)

Define our first Lyapunov function as

$$\begin{array}{c}\hfill V_1={\displaystyle \frac{1}{2}}{\mathbf{e}}_1^T{\mathbf{e}}_1\end{array}$$

(12)

and compute its time derivative, we have

$$\begin{array}{c}\hfill {\dot{V}}_1=-{\alpha }_1+{\mathbf{e}}_1^T(\mathbf{v}-{\mathbf{R}}^T{\dot{\mathbf{p}}}_d+k_1{\mathbf{e}}_1)\end{array}$$

(13)

where ${\alpha }_1=k_1{\mathbf{e}}_1^T{\mathbf{e}}_1$ is positive definite, and $k_1$ denotes gain, which is a positive value.

In order to continue to use the backstepping technique, a second error term is defined as

$$\begin{array}{c}\hfill {\mathbf{e}}_2=\mathbf{v}-{\mathbf{R}}^T{\dot{\mathbf{p}}}_d+k_1{\mathbf{e}}_1-\mathit{\eta },\end{array}$$

(14)

and its integral term

$$\begin{array}{c}\hfill {\mathbf{e}}_{2n}={\int }_0^t{\mathbf{e}}_{2}dt\end{array}$$

(15)

where $\mathit{\eta }={[{\eta }_1,{\eta }_2]}^T,{\eta }_1\ne 0$ is a constant vector. Define our second Lyapunov function as

$$\begin{array}{c}\hfill V_2=V_1+{\displaystyle \frac{1}{2}}{\mathbf{e}}_2^T{\mathbf{e}}_2+{\displaystyle \frac{1}{2}}{\mathbf{e}}_{2n}^T{\mathbf{e}}_{2n},\end{array}$$

(16)

and its time derivative yields

$$\begin{array}{cc}{\dot{V}}_2\hfill & =-{\alpha }_2+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T(-\mathbf{S}\left(\omega \right)\mathit{\eta }+m^{-1}T{\mathbf{n}}_1-{\mathbf{R}}^T{\ddot{\mathbf{p}}}_d+k_1({\mathbf{e}}_2-k_1{\mathbf{e}}_1+\mathit{\eta })\hfill \\ & +k_2{\mathbf{e}}_2+{\Delta }_v+{\mathbf{e}}_{2n})\hfill \end{array}$$

(17)

where ${\alpha }_2=k_1{\mathbf{e}}_1^T{\mathbf{e}}_1+k_2{\mathbf{e}}^T{\mathbf{e}}_2$, which is positive definite, and $k_2$ is a positive number. It is relevant to point out that ${\mathbf{e}}_{2n}$ is introduced to eliminate the external slow-varying disturbances that act on the dynamic of the linear velocity $\mathbf{v}$. Moreover, notice that we do not know ${\Delta }_v$, thereby we use RBFNNs mentioned before to approximate it, given by

$$\begin{array}{c}\hfill {\Delta }_v={\mathbf{W}}_{d1}^T\mathit{\sigma }\left({\mathbf{x}}_1\right)+{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)\end{array}$$

(18)

where ${\mathbf{x}}_1={[1,{\mathbf{v}}^T]}^T\in {\mathbb{R}}^3$. Moreover, notice that

$$\begin{array}{c}\hfill -\mathbf{S}\left(\omega \right)\mathit{\eta }+m^{-1}T{\mathbf{n}}_1=\mathbf{N}\mathbf{I}\end{array}$$

(19)

where

$$\mathbf{N}=\left[\begin{array}{cc}{m}^{-1}& {\eta }_2\\ 0& -{\eta }_1\end{array}\right].\mathbf{I}={[T,\omega ]}^T.$$

Rewrite Equation (17), and we have

$$\begin{array}{c}\hfill {\dot{V}}_2=-{\alpha }_2+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T\left(\mathbf{N}\mathbf{I}+\mathit{\beta }+{\mathbf{W}}_{d1}^T\mathit{\sigma }\left({\mathbf{x}}_1\right)\right)+{\mathbf{e}}_2^T{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)\end{array}$$

(20)

where $\mathit{\beta }=-{\mathbf{R}}^T{\ddot{\mathbf{p}}}_d+k_1({\mathbf{e}}_2-k_1{\mathbf{e}}_1+\mathit{\eta })+k_2{\mathbf{e}}_2+{\mathbf{e}}_{2n}$.

Now, we can define our third Lyapunov function as

$$\begin{array}{c}\hfill V_3=V_2+{\displaystyle \frac{1}{2}}\left({\tilde{\mathbf{W}}}_{d1}^T{\Gamma }_1^{-1}{\tilde{\mathbf{W}}}_{d1}\right)\end{array}$$

(21)

where ${\tilde{\mathbf{W}}}_{d1}={\mathbf{W}}_d-{\widehat{\mathbf{W}}}_{d1}$ denotes the estimation error, ${\Gamma }_1=({\lambda }_{11},and{\lambda }_{12})$ is a matrix, where ${\lambda }_{11}$ and ${\lambda }_{12}$ are positive values. Compute the time derivative of $V_3$, and one obtains

$$\begin{array}{cc}{\dot{V}}_3\hfill & =-{\alpha }_2+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)+{\mathbf{e}}_2^T\left(\mathbf{N}\mathbf{I}+\mathit{\beta }+{\widehat{\mathbf{W}}}_{d1}^T\mathit{\sigma }\left({\mathbf{x}}_1\right)\right)\hfill \\ & +{\mathbf{e}}_2^T{\tilde{\mathbf{W}}}_{d1}^T\sigma \left({\mathbf{x}}_1\right)-\left({\tilde{\mathbf{W}}}_{d1}^T{\Gamma }_1^{-1}{\dot{\widehat{\mathbf{W}}}}_{d1}\right).\hfill \end{array}$$

(22)

Therefore, we choose our desired input ${\mathbf{I}}_d$ as

$$\begin{array}{c}\hfill {\mathbf{I}}_d=-{\mathbf{N}}^{-1}\left(\mathit{\beta }+{\widehat{\mathbf{W}}}_{d1}^T\mathit{\sigma }\left({\mathbf{x}}_1\right)\right),\end{array}$$

(23)

and thereby our first controller T is chosen as

$$\begin{array}{c}\hfill T={\mathbf{n}}_1^T{\mathbf{I}}_d.\end{array}$$

(24)

Correspondingly, the desired angular velocity is

$$\begin{array}{c}\hfill {\omega }_d={\mathbf{n}}_2^T{\mathbf{I}}_d\end{array}$$

(25)

where ${\mathbf{n}}_2={[0,1]}^T$.

Notice that, if the updated law for ${\dot{\widehat{\mathbf{W}}}}_{d1}$ is set as

$$\begin{array}{c}\hfill {\dot{\widehat{\mathbf{W}}}}_{d1}={\Gamma }_1\sigma \left({\mathbf{x}}_1\right){\mathbf{e}}_2^T,\end{array}$$

(26)

we cannot ensure that it is bounded by ${\mathbf{W}}_{\max }$. To solve this problem, a projection operator, which is Lipschitz continuous [25], is applied in our case, which is given by

$$\begin{array}{cc}\mathbf{proj}(\rho ,\widehat{\ell })\hfill & =\left\{\begin{array}{cc}\rho \hfill & \mathrm{if} \Theta (\widehat{\ell })\le 0\hfill \\ \rho \hfill & \mathrm{if} \Theta (\widehat{\ell })\ge 0 \mathrm{and} {\Theta }_{\widehat{\ell }}(\widehat{\ell })\rho \le 0\hfill \\ (1-\Theta (\widehat{\ell }))\rho \hfill & \mathrm{if} \Theta (\widehat{\ell })>0 \mathrm{and} {\Theta }_{\widehat{\ell }}(\widehat{\ell })\rho >0\hfill \end{array}\right.\hfill \end{array}$$

(27)

where

$$\begin{array}{c}\hfill \Theta (\widehat{\ell })={\displaystyle \frac{{\widehat{\ell }}^2-{\ell }_{\max }^2}{{ϵ}^2+2ϵ{\ell }_{\max }^2}},{\Theta }_{\widehat{\ell }}(\widehat{\ell })={\displaystyle \frac{\partial \Theta (\widehat{\ell })}{\partial \widehat{\ell }}},\end{array}$$

(28)

with the following condition: if $\dot{\widehat{\ell }}=proj(\rho ,\widehat{\ell })$ and $\widehat{\ell }(t_0)\le {\ell }_{\max }$, then

(1)

$\widehat{\ell }\le {\ell }_{\max }+ϵ,\forall 0\le t<\infty$;

(2)

$\mathbf{proj}(\rho ,\widehat{\ell })$ is Lipschitz continuous;

(3)

$|\mathbf{proj}(\rho ,\widehat{\ell })|\le \left|\rho \right|$;

(4)

$\tilde{\ell }\mathbf{proj}(\rho ,\widehat{\ell })\ge \tilde{\ell }\rho$, where $\tilde{\ell }=\ell -\widehat{\ell }$.

Therefore, to make sure all elements of ${\widehat{\mathbf{W}}}_{d1}$ are upper-bounded, the update law for ${\dot{\widehat{\mathbf{W}}}}_{d1}$ is finally set as

$$\begin{array}{c}\hfill {\dot{\widehat{\mathbf{W}}}}_{d1}={\Gamma }_1\mathbf{proj}\left(\sigma \left({\mathbf{x}}_1\right){\mathbf{e}}_2,{\widehat{\mathbf{W}}}_{d1}\right).\end{array}$$

(29)

To keep using the backstepping technique, we define a new error term

$$\begin{array}{c}\hfill {\mathbf{e}}_3=\omega -{\omega }_d,\end{array}$$

(30)

and its corresponding integral term

$$\begin{array}{c}\hfill {\mathbf{e}}_{3n}={\int }_0^t{\mathbf{e}}_{3}dt.\end{array}$$

(31)

Then, we define a new Lyapunov function as

$$\begin{array}{c}\hfill V_4=V_3+{\displaystyle \frac{1}{2}}{\mathbf{e}}_3^T{\mathbf{e}}_3+{\displaystyle \frac{1}{2}}{\mathbf{e}}_{3n}^T{\mathbf{e}}_{3n}.\end{array}$$

(32)

Compute its time derivatives, substitute Equation (29) into Equation (24), and combine the 4th property of projector $(\tilde{\ell }\mathbf{proj}(\rho ,\widehat{\ell })\ge \tilde{\ell }\rho )$, and one obtains

$$\begin{array}{cc}{\dot{V}}_4\hfill & \le -{\alpha }_3+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)+{\mathbf{e}}_3^T({\mathbf{n}}_2^T{\mathbf{N}}^T{\mathbf{e}}_2-{\displaystyle \frac{1}{J}}\tau +{\mathbf{n}}_2^T{\mathbf{N}}^{-1}(\dot{\mathit{\beta }}+{\dot{\widehat{\mathbf{W}}}}_{d1}\mathit{\sigma }\left({\mathbf{x}}_1\right)\hfill \\ & +{\widehat{\mathbf{W}}}_{d1}\dot{\mathit{\sigma }}\left({\mathbf{x}}_1\right))+k_3{\mathbf{e}}_3+{\mathbf{e}}_{3n}+{\Delta }_{\omega })\hfill \end{array}$$

(33)

where ${\alpha }_3={\alpha }_2+k_3{\mathbf{e}}_3^T{\mathbf{e}}_3\ge 0$, $k_3>0$. However, it is noted that both $\dot{\mathit{\beta }}$ and $\dot{\mathit{\sigma }}\left({\mathbf{x}}_1\right)$ contain unmodeled term ${\Delta }_v$, so we need to separate ${\Delta }_v$ out from $\dot{\mathit{\beta }}$ and $\dot{\mathit{\sigma }}\left({\mathbf{x}}_1\right)$. After that, we can use RBFNNs to estimate it. Similar to ${\Delta }_{\omega }$, we also need to approximate it. This is given by,

$$\begin{array}{cc}\hfill {\Delta }_v& ={\mathbf{W}}_{d2}^T\mathit{\sigma }\left({\mathbf{x}}_2\right)+{\mathit{ϵ}}_2\left({\mathbf{x}}_2\right)\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\Delta }_{\omega }& ={\mathbf{W}}_{d3}^T\mathit{\sigma }\left({\mathbf{x}}_3\right)+{\mathit{ϵ}}_3\left({\mathbf{x}}_3\right)\hfill \end{array}$$

(35)

where ${\mathbf{x}}_2={[1,{\mathbf{v}}^T]}^T\in {\mathbb{R}}^3$ and ${\mathbf{x}}_3={[1,{\mathbf{v}}^T,\omega ]}^T\in {\mathbb{R}}^4$.

Now we define our last Lyapunov function as

$$\begin{array}{c}\hfill V_5=V_4+{\displaystyle \frac{1}{2}\mathrm{tr}}\left({\tilde{\mathbf{W}}}_{d2}^T{\Gamma }_2^{-1}{\tilde{\mathbf{W}}}_{d2}\right)+{\displaystyle \frac{1}{2}}\mathrm{tr}\left({\tilde{\mathbf{W}}}_{d3}^T{\Gamma }_3^{-1}{\tilde{\mathbf{W}}}_{d3}\right)\end{array}$$

(36)

where ${\tilde{\mathbf{W}}}_{d2}={\mathbf{W}}_d-{\widehat{\mathbf{W}}}_{d2}$, ${\tilde{\mathbf{W}}}_{d3}={\mathbf{W}}_d-{\widehat{\mathbf{W}}}_{d3}$, ${\Gamma }_2=\mathrm{diag}({\lambda }_{21},{\lambda }_{22})$, and ${\Gamma }_3=\mathrm{diag}({\lambda }_{31},{\lambda }_{32})$ are positive definite gain matrices. Then, we compute the time derivative of $V_5$ as

$$\begin{array}{cc}{\dot{V}}_5\hfill & \le -{\alpha }_3+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)+{\mathbf{e}}_3^T\left({\mathbf{n}}_2^T\mathbf{M}\right){\mathit{ϵ}}_2\left({\mathbf{x}}_2\right)+{\mathbf{e}}_3^T{\mathit{ϵ}}_3\left({\mathbf{x}}_3\right)+{\mathbf{e}}_3^T({\mathbf{n}}_2^T{\mathbf{N}}^T{\mathbf{e}}_2+k_3{\mathbf{e}}_3\hfill \\ & -{\displaystyle \frac{\tau }{J}}-{\mathbf{n}}_2^T{\widehat{\dot{\mathbf{I}}}}_d-{\mathbf{n}}_2^T\mathbf{M}{\widehat{\mathbf{W}}}_{2d}^T\sigma \left({\mathbf{x}}_2\right)+{\widehat{\mathbf{W}}}_{3d}^T\sigma \left({\mathbf{x}}_3\right)+{\mathbf{e}}_{3n})+{\mathbf{e}}_3^T{\mathbf{n}}_2^T\mathbf{M}{\tilde{\mathbf{W}}}_{2d}^T\sigma \left({\mathbf{x}}_2\right)\hfill \\ & +{\mathbf{e}}_3^T{\tilde{\mathbf{W}}}_3^T\sigma \left({\mathbf{x}}_3\right)-\mathrm{tr}\left({\tilde{\mathbf{W}}}_{d2}^T{\Gamma }_2^{-1}{\dot{\widehat{\mathbf{W}}}}_{d2}\right)-\mathrm{tr}\left({\tilde{\mathbf{W}}}_{d3}^T{\Gamma }_3^{-1}{\dot{\widehat{\mathbf{W}}}}_{d3}\right)\hfill \end{array}$$

(37)

where

$$\begin{array}{c}\hfill {\widehat{\dot{\mathbf{I}}}}_d=\widehat{\dot{\mathit{\beta }}}+{\dot{\widehat{\mathbf{W}}}}_{d1}\mathit{\sigma }\left({\mathbf{x}}_1\right)+{\widehat{\mathbf{W}}}_{d1}\widehat{\dot{\mathit{\sigma }}}\left({\mathbf{x}}_1\right)\end{array}$$

(38)

and

$$\begin{array}{c}\hfill \mathbf{M}={\mathbf{N}}^{-1}{\widehat{\mathbf{W}}}_{d1}^T\mathbf{G}+(k_1+k_2){\mathbf{N}}^{-1}\end{array}$$

(39)

with $\mathbf{G}={[{\mathbf{0}}_{2\times 1}^T,{\mathbf{n}}_1{\mathit{\sigma }}^{\prime }\left({\mathbf{n}}_1^T\mathbf{v}\right),and{\mathbf{n}}_1{\mathit{\sigma }}^{\prime }\left({\mathbf{n}}_2^T\mathbf{v}\right)]}^T$. Then, we define our second control law, torque, as

$$\begin{array}{c}\hfill \tau =J(-{\mathbf{n}}_2^T{\widehat{\dot{\mathbf{I}}}}_d+{\mathbf{n}}_2^T{\mathbf{N}}^T{\mathbf{e}}_2+k_3{\mathbf{e}}_3+{\mathbf{n}}_2^T\mathbf{M}{\widehat{\mathbf{W}}}_{d2}^T\mathit{\sigma }\left({\mathbf{x}}_2\right)+{\widehat{\mathbf{W}}}_{d3}^T\mathit{\sigma }\left({\mathbf{x}}_3\right)+{\mathbf{e}}_{3n}),\end{array}$$

(40)

and estimate laws for ${\dot{\widehat{\mathbf{W}}}}_{d2}$ and ${\dot{\widehat{\mathbf{W}}}}_{d3}$

$$\begin{array}{cc}\hfill {\dot{\widehat{\mathbf{W}}}}_{d2}& ={\Gamma }_2\mathbf{proj}\left(\mathit{\sigma }\left({\mathbf{x}}_2\right)\left({\mathbf{n}}_2^T\mathbf{M}\right){\mathbf{e}}_3^T{\widehat{\mathbf{W}}}_{d2}\right),\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\dot{\widehat{\mathbf{W}}}}_{d3}& ={\Gamma }_3\mathbf{proj}\left(\mathit{\sigma }\left({\mathbf{x}}_3\right){\mathbf{e}}_3^T,{\widehat{\mathbf{W}}}_{d3}\right).\hfill \end{array}$$

(42)

Substitute Equation (40)–(42) into Equation (37), one obtains

$$\begin{array}{cc}\hfill {\dot{V}}_5\le & -{\alpha }_3+{\mathbf{e}}_1^T\mathit{\eta }+{\mathbf{e}}_2^T{\mathit{ϵ}}_1\left({\mathbf{x}}_1\right)+{\mathbf{e}}_3^T\left({\mathbf{n}}_2^T\mathbf{M}\right){\mathit{ϵ}}_2\left({\mathbf{x}}_2\right)+{\mathbf{e}}_3^T{\mathit{ϵ}}_3\left({\mathbf{x}}_3\right).\hfill \end{array}$$

(43)

where ${\alpha }_3=k_1{\mathbf{e}}_1^T{\mathbf{e}}_1+k_2{\mathbf{e}}_2^T{\mathbf{e}}_2+k_3{\mathbf{e}}_3^T{\mathbf{e}}_3$.

3.1. Stability Analysis

Theorem 1.

For a single mobile robot, by applying control laws, Equations (24) and (40), and updated laws, Equation (29), Equation (41), and Equation (42), for any large initial position, the robot will converge to the neighborhood of its corresponding desired path ${\mathit{p}}_d\left(\gamma \right)$, whose partial derivatives with respect to γ are all bounded. As a consequence, global uniformly ultimately boundedness is achieved.

Proof. 

Let’s go back to Equation (43). Rewrite it, and we obtain

$$\begin{array}{cc}{\dot{V}}_5\hfill & \le -{\mathbf{X}}^T\mathbf{K}\mathbf{X}+{\mathbf{X}}^T\rho \hfill \\ & \le {\left|\mathbf{X}\right|}^T\left(-k_{\min }\left|\mathbf{X}\right|+\left|\rho \right|\right)\hfill \end{array}$$

(44)

where $\mathbf{X}=[{\mathbf{e}}_1^T,{\mathbf{e}}_2^T,{\mathbf{e}}_3^T]$, $k_{\min }$ is the smallest eigenvalue of $\mathbf{K}=(k_1,k_2,k_3)$, and $\rho ={[{\mathit{\eta }}^T,{\mathit{ϵ}}_1{\left({\mathbf{x}}_1\right)}^T,\left({\mathbf{n}}_2^T\mathbf{M}\right){\mathit{ϵ}}_2\left({\mathbf{x}}_2\right),{\mathit{ϵ}}_3\left({\mathbf{x}}_3\right)]}^T$, which is bounded due to the fact that $\left|\right|{\mathit{ϵ}}_i\left(\mathbf{x}\right)\left|\right|\le {ϵ}_{\max }$, and the upper bound of $\mathit{\rho }$ is

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\max }={[{\mathit{\eta }}^T,{\mathit{ϵ}}_{1\max }{\left({\mathbf{x}}_1\right)}^T,\left({\mathbf{n}}_2^T\mathbf{M}\right){\mathit{ϵ}}_{2\max }\left({\mathbf{x}}_2\right),{\mathit{ϵ}}_{3\max }\left({\mathbf{x}}_3\right)]}^T.\end{array}$$

(45)

Thereby, we can obtain that ${\dot{V}}_5$ is negative for $\left|\right|\mathbf{X}\left|\right|\ge \left|\right|{\mathit{\rho }}_{\max }/k_{\min }\left|\right|$, which can be made as small as possible by tuning the value of $k_{\min }$. As a result, the system is uniformly ultimate bounded, global uniformly ultimate boundedness is achieved. ☐

4. Consensus Strategy

Building upon the work of [14], the proposed solution is given by

$$\begin{array}{cc}\hfill {\dot{\gamma }}_i& =v_d\mathbf{1}-a_1\sum _{j\in N_i}({\gamma }_i-{\gamma }_j)+z_i\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\dot{z}}_i& =-a_2{z}_i+a_3\sum _{j\in N_i}({\gamma }_i-{\gamma }_j)\hfill \end{array}$$

(47)

where $a_1,a_2,and{a}_3$ are positive numbers, and $v_d$ denotes the desired value of ${\dot{\gamma }}_i$. It is relevant to point out that $z_i$ can be viewed as an auxiliary state that helps n paths to reach consensus.

4.1. Stability Analysis

Theorem 2.

For a group of n desired paths ${\mathit{p}}_d\left({\gamma }_i\right)(i=1,2,\dots ,n)$, by applying Equations (46) and (47), ${\gamma }_i-{\gamma }_j$ and ${\dot{\gamma }}_i-{\dot{\gamma }}_d$ converge to zero.

Proof. 

We first choose Laplacian matrix $\mathcal{L}$, and define the coordinate error as

$$\begin{array}{c}\hfill {\Gamma }_e=\mathcal{L}\Lambda \end{array}$$

(48)

where $\Lambda ={\left[{\gamma }_i\right]}_{n\times 1}$. Rewrite Equations (46) and (47), and one obtains

$$\begin{array}{cc}\hfill {\dot{\Gamma }}_e& =-{\mathbf{A}}_1\mathcal{L}{\Gamma }_e+\mathcal{L}\mathbf{Z},\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill \dot{\mathbf{Z}}& =-{\mathbf{A}}_2\mathbf{Z}+{\mathbf{A}}_3{\Gamma }_e,\hfill \end{array}$$

(50)

where ${\mathbf{A}}_i=\left(a_i\right),(i=1,2,3)$ is a positive definite matrix. Define $\mathbf{x}={[{\Gamma }_e,\mathbf{Z}]}^T$, and rewrite Equations (49) and (50), and we have

$$\begin{array}{c}\hfill \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}\end{array}$$

(51)

where

$$\begin{array}{c}\hfill \mathbf{A}=\left[\begin{array}{cc}-{\mathbf{A}}_1\mathcal{L}& \mathcal{L}\\ {\mathbf{A}}_3& -{\mathbf{A}}_2\end{array}\right].\end{array}$$

(52)

In order to ensure that Equation (51) is stable, we need to guarantee that all the eigenvalues of $\mathbf{A}$ have negative or zero real parts and all the Jordan blocks corresponding to eigenvalues with zero real parts are $1\times 1$. We consider n agents, where

$$\begin{array}{c}\hfill \mathcal{L}={\left|\begin{array}{cccccc}0& 0& 0& \cdots & 0& 0\\ -1& 1& 0& \cdots & 0& 0\\ 0& -1& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & -1& 1\end{array}\right|}_{n\times n}.\end{array}$$

(53)

For the sake of saving space, here we just present the eigenvalues of $\mathbf{A}$ directly, they are ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$, and ${\lambda }_4$, and their corresponding multiplicities are 1, 1, $n-1$, and $n-1$, with

$$\begin{array}{cc}\hfill {\lambda }_1& =0\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill {\lambda }_2& =-a_2\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\lambda }_3& =-(a_1+a_2)/2+\sqrt{a_1^2+a_2^2-2a_1{a}_2+4a_3}/2\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill {\lambda }_4& =-(a_1+a_2)/2-\sqrt{a_1^2+a_2^2-2a_1{a}_2+4a_3}/2.\hfill \end{array}$$

(57)

Choosing $a_1,a_2,and{a}_3$ properly, we can guarantee that ${\lambda }_2,{\lambda }_3$, and ${\lambda }_4$ are negative-definite. Moreover, we also have that the Jordan block corresponding to ${\lambda }_1$ is $1\times 1$. As a consequence, Equation (51) is stable [26]. ☐

To summarize, a fleet of n desired paths can progress in a specific formation, while each individual mobile robot converges to the neighborhood of its corresponding desired path. As a result, all those robots can move forward in a specific formation.

5. Simulation Results

In this section, we present simulation results with two different communication graphs including a cascade-directed communication graph (CDCG) and a parallel-directed communication graph (PDCG). Figure 3 and Figure 4 show the sketches of the control blocks in Simulink/Matlab.

5.1. The Cascade-Directed Communication Graph

The CDCG used in this study is shown in Figure 5, where agent 1 can be viewed as the leader. Its corresponding Laplacian matrix ${\mathcal{L}}_1$ is

$$\begin{array}{c}\hfill {\mathcal{L}}_1=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ -1& 1& 0& 0& 0\\ 0& -1& 1& 0& 0\\ 0& 0& -1& 1& 0\\ 0& 0& 0& -1& 1\end{array}\right].\end{array}$$

(58)

The desired paths are defined as

$$\begin{array}{c}\hfill {\mathbf{p}}_{di}\left({\gamma }_i\right)=R_i\left[\begin{array}{c}\cos \left({\gamma }_i\right)\\ \sin \left({\gamma }_i\right)\end{array}\right]\left(m\right)\end{array}$$

(59)

where $R_i=6-i,(i=1,2,\dots ,5)$ denotes the radius of the circle, and unmodeled terms ${\Delta }_V$ and ${\Delta }_{\omega }$ are denoted by ${[0.1u^2+0.01uv,0.01uv+0.1v^2]}^T$ and $0.01uv+0.05u\omega +0.01v\omega$, respectively. The parameters used herein are as follows: $m=0.6$, $J=0.1$, $c_i=2$, $a_1=1$, $a_2=2$, $a_3=1$, $k_1=6,$$k_2=3,k_3=2,{\Gamma }_1=(0.6,0.6),{\Gamma }_2=(0.01,0.01),{\Gamma }_3=(0.022,0.022)$, and $\mathit{\eta }={[0.2,0]}^T$.

Figure 6a shows the actual trajectory of the mobile robots, and we can see that they move forward in a line formation. Moreover, Figure 7a,b display the convergence of $\left|\right|{\mathbf{e}}_{1i}\left|\right|$ and $\left|\right|{\mathbf{e}}_{2i}\left|\right|$, respectively, showing that all of them converge to the neighborhood of zero. The consensus performances are shown in Figure 8a,b, showing that ${\gamma }_{ij}={\gamma }_i-{\gamma }_j$ converges to zero, where agent j is the neighbor of agent i. Moreover, we can also see that ${\dot{\gamma }}_i$ converges to the desired value ${\dot{\gamma }}_d=0.5$. It is important to point out that, in this work, we chosen agent 1 as our leader, and this satisfies ${\dot{\gamma }}_1={\dot{\gamma }}_d$. The approximate performance of RBFNNs can be found in Figure 6b, where the blue lines denote the real values of ${\Delta }_v$ and ${\Delta }_{\omega }$, while the red lines represent their estimates ${\widehat{\Delta }}_v$ and ${\widehat{\Delta }}_{\omega }$. Thus, both estimates converge to their corresponding real values.

5.2. Parallel Communication Graph

The parallel communication graph is depicted in Figure 9.

In this case, the Laplacian matrix ${\mathcal{L}}_2$ is

$$\begin{array}{c}\hfill {\mathcal{L}}_2=\left[\begin{array}{ccccc}-1& 1& 0& 0& 0\\ -1& 0& 1& 0& 0\\ 0& -1& 0& 1& 0\\ 0& 0& -1& 0& 1\\ 0& 0& 0& 0& 0\end{array}\right].\end{array}$$

(60)

Moreover, it is noted that the graph presented in Figure 9 can be viewed as the combination of two cascade-directed graphs, where $\mathbf{1}\to \mathbf{2}\to \mathbf{4}$ and $\mathbf{1}\to \mathbf{3}\to \mathbf{5}$ and where agent $\mathbf{1}$ is the leader whose state is known. The desired paths are as follows:

(1)

1st vehicle:

$${\mathbf{p}}_{d1}\left({\gamma }_1\right)=R_1\left[\begin{array}{c}\cos \left({\gamma }_1\right)\\ \sin \left({\gamma }_1\right)\end{array}\right]\left(\mathrm{m}\right)$$

(61)

(2)

2nd vehicle:

$${\mathbf{p}}_{d2}\left({\gamma }_2\right)=R_2\left[\begin{array}{c}\cos ({\gamma }_2-\pi /24)\\ \sin \left({\gamma }_2\right)-\pi /24\end{array}\right]\left(\mathrm{m}\right)$$

(62)

(3)

3rd vehicle:

$${\mathbf{p}}_{d3}\left({\gamma }_3\right)=R_3\left[\begin{array}{c}\cos ({\gamma }_3-\pi /12)\\ \sin ({\gamma }_3-\pi /12)\end{array}\right]\left(\mathrm{m}\right)$$

(63)

(4)

4th vehicle:

$${\mathbf{p}}_{d4}\left({\gamma }_4\right)=R_4\left[\begin{array}{c}\cos ({\gamma }_4-\pi /24)\\ \sin ({\gamma }_4-\pi /24)\end{array}\right]\left(\mathrm{m}\right)$$

(64)

(5)

5th vehicle:

$${\mathbf{p}}_{d5}\left({\gamma }_5\right)=R_5\left[\begin{array}{c}\cos ({\gamma }_5-\pi /12)\\ \sin ({\gamma }_5-\pi /12)\end{array}\right]\left(\mathrm{m}\right)$$

(65)

where $R_1=3,R_2=4,R_3=2,R_4=5,$ and $R_5=1.$ Figure 10a displays the actual paths of the robots. From Figure 11a,b, we can obtain that the norm of position errors and linear velocity errors converge to a ball centered at the origin. Moreover, Figure 12a,b show the performance of the consensus strategy introduced herein.

It is interesting to remark that, by using the proposed consensus strategy—Equations (46) and (47)—a consensus is researched if and only if the directed graph has a directed spanning tree [21]. However, in our case, the root must be the leader whose states are known beforehand.

6. Conclusions

In this paper, we mainly focus on designing coordinated control algorithms for multiple agents, where a group of underactuated hovercrafts were chosen as test platforms. In order to testify the efficiency of the devised control strategy, we implemented it by using Simulink/Matlab. Moreover, it is necessary to point out that agents can also be mobile robots, unmanned air vehicles, etc. For a single vehicle, we used RBFNNs to approximate unmodeled terms and introduce integral terms, which can improve the robustness of the controller. For multiple vehicles, we consider directed topology under the assumption that the communication among vehicles are continuous.

With respect to our future works, we plan to (i) use deep neural networks to estimate unmodeled terms so as to enhance the performance of approximation, (ii) build a mathematical model for external disturbance, such as winds, waves, or currents, (iii) take into account time-delays when we develop communication strategy for vehicles, and (iv) propose collision–avoidance algorithms so that we can ensure the operation is safe.