**9. Conclusions and Future Work**

This paper addressed the problem of encircling an environmental boundary caused by a chemical spill using a team of robots composed of an aerial quadrotor and Medusa marine vehicles. The path following problem was introduced, and a non-linear control law derived for the ASV, exploiting the technique described in P. Aguiar and F. Vanni [10–12]. Inspired by this control law, a new one was derived for a quadrotor following the same methodology, with some key differences due to the nature of the aircraft. For the section that followed, the CPF problem was formulated and a proposal to solve the problem was presented, such that the synchronisation controller was distributed and the same for all vehicles (aerial and marine) using event-triggered communications based on previous work by N. Hung and F. Rego [13]. In addition, a new real-time path planning algorithm was developed that made use of the camera sensor onboard of the quadrotor to have a local view of the boundary and generate a point cloud expressed in the inertial frame. This data was then used to solve an optimisation problem which generates a B-spline-based path that grows dynamically as the drone moves along the boundary and acquires more data. The path is then shared with all ASV vehicles in the network in real time. The proposed algorithms were implemented in ROS, and a 3-D virtual scenario was generated, allowing for a mixture of real and simulated results. Future work includes making the height at which the quadrotor operates dynamic and introducing curvature limits as inequality constraints to the path planning problem, as well as obstacle avoidance before carrying out integrated experiments with real vehicles.

**Author Contributions:** The individual contributions of the authors are as follows: Conceptualisation, software, formal analysis, investigation, writing, review, and editing by M.J.; Conceptualisation, validation, review, editing, project administration, and funding acquisition by A.P. and R.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially funded by the H2020-INFRAIA-2017-1-two-stage EUMarineRobots-Marine Robotics Research Infrastructure Network (Grant agreement ID: 731103), the H2020-FETPROACT-2020-2 RAMONES-Radioactivity Monitoring in Ocean Ecosystems (Grant agreement ID: 101017808), H2020-MSCA-RISE-2018 ECOBOTICS.SEA-Bio-inspired Technologies for a Sustainable Marine Ecosystem (Grant agreement ID: 824043).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The path planning library that implements the algorithms proposed in Section 6 is available at https://github.com/MarceloJacinto/BSplineFit (accessed on 8 March 2022).

**Acknowledgments:** The authors of this work would like to express a deep gratitude to Nguyen Hung, Francisco Rego, João Quintas, and João Cruz for all the support provided during the preparation of the results presented in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Proof of Preposition 1**

**Proof.** Consider the candidate Lyapunov Function given by:

$$V\_1(\mathbf{e}\_p) = \frac{1}{2}(\mathbf{e}\_p - \delta)^T(\mathbf{e}\_p - \delta). \tag{A1}$$

Taking the first time derivative of (A1) and replacing in (17) and (14) leads to:

$$\begin{split} \dot{V}\_{1}(\mathbf{e}\_{p}) &= (\mathbf{e}\_{p} - \boldsymbol{\delta})^{T} \Big( -S(r)(\mathbf{e}\_{p} - \boldsymbol{\delta}) + \Delta \mathbf{u} + \begin{bmatrix} 0 \\ v \end{bmatrix} \\ &+ \boldsymbol{\sigma}\_{c} - \, ^{B}\_{ll} \mathcal{R}(\boldsymbol{\psi}) \frac{\partial \mathbf{p}\_{d}(\boldsymbol{\gamma})}{\partial \boldsymbol{\gamma}} \Big( \boldsymbol{e}\_{\boldsymbol{\gamma}} + \boldsymbol{v}\_{d}(\boldsymbol{\gamma}, t) \Big) ). \end{split} \tag{A2}$$

Taking into account the properties of the skew-symmetric matrix S:

$$(\mathbf{e}\_p - \boldsymbol{\delta})^T \mathbf{S}(r) (\mathbf{e}\_p - \boldsymbol{\delta}) = 0. \tag{A3}$$

Replacing (20) and (21) in *V*˙ <sup>1</sup> yields:

$$\begin{split} \dot{V}\_{1}(\mathbf{e}\_{p}) &= (\mathbf{e}\_{p} - \boldsymbol{\delta})^{T} \Big( \Delta(\hat{\mathbf{u}} + \mathbf{u}\_{d}) + \begin{bmatrix} 0 \\ \upsilon \end{bmatrix} + \vec{\sigma}\_{\boldsymbol{\varepsilon}} + \mathfrak{e}\_{\boldsymbol{\varepsilon}} - \,^{\mathrm{B}}\_{\mathrm{I}} \mathcal{R}(\boldsymbol{\psi}) \frac{\partial \mathbf{p}\_{d}(\boldsymbol{\gamma})}{\partial \boldsymbol{\gamma}} \big( \boldsymbol{e}\_{\boldsymbol{\gamma}} + \boldsymbol{v}\_{d}(\boldsymbol{\gamma}, t) \big) \Big) \\ &= -(\mathbf{e}\_{p} - \boldsymbol{\delta})^{T} \mathcal{K}\_{p} \boldsymbol{\sigma} (\mathbf{e}\_{p} - \boldsymbol{\delta}) - (\mathbf{e}\_{p} - \boldsymbol{\delta})^{T} \prescript{B}{}{\mathrm{I}} \mathcal{R}(\boldsymbol{\psi}) \frac{\partial \mathbf{p}\_{d}(\boldsymbol{\gamma})}{\partial \boldsymbol{\gamma}} \boldsymbol{e}\_{\boldsymbol{\gamma}} + \Delta \bar{\mathbf{u}} + \vec{\sigma}\_{\boldsymbol{\varepsilon}}. \end{split} \tag{A4}$$

By taking a backstepping approach, consider a second candidate Lyapunov function:

$$V\_2(\mathbf{e}\_p, \mathbf{e}\_\gamma) = V\_1(\mathbf{e}\_p) + \frac{1}{2}e\_\gamma^2. \tag{A5}$$

Taking the first derivative and replacing in the control law (22) for the virtual target, we obtain:

$$\begin{split} \mathcal{V}\_{2}(\mathbf{e}\_{p},\mathbf{e}\_{\gamma}) &= \dot{\mathcal{V}}\_{1}(\mathbf{e}\_{p}) + e\_{\gamma}(\dot{\gamma} - \dot{\upsilon}\_{d}(\gamma,t)) \\ &= -\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right)^{T} \mathcal{K}\_{p}\boldsymbol{\sigma}\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right) - k\_{\gamma}e\_{\gamma}^{2} + \Delta \dot{\mathbf{u}} + \tilde{\boldsymbol{\sigma}}\_{\mathcal{L}} \\ &\leq -\left(1 - \theta\right)\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right)^{T} \mathcal{K}\_{p}\boldsymbol{\sigma}\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right) - \theta\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right)^{T} \mathcal{K}\_{p}\boldsymbol{\sigma}\left(\mathbf{e}\_{p} - \boldsymbol{\delta}\right) \\ &\quad - k\_{\gamma}|\mathbf{e}\_{\gamma}|^{2} + ||\mathbf{e}\_{p} - \boldsymbol{\delta}|| \left|\Delta \tilde{\mathbf{u}} + \tilde{\mathbf{v}}\_{\mathcal{L}}\right|. \end{split} \tag{A6}$$

where 0 < *θ* < 1. The term:

$$\begin{split} & -\theta \left( \mathbf{e}\_p - \boldsymbol{\delta} \right)^T K\_p \sigma \left( \mathbf{e}\_p - \boldsymbol{\delta} \right) + \left\| \mathbf{e}\_p - \boldsymbol{\delta} \right\| \left\| \Delta \tilde{\mathbf{u}} + \vec{\sigma}\_{\varepsilon} \right\| \\ & = -\theta (\mathbf{e}\_p - \boldsymbol{\delta})^T K\_p \frac{\mathbf{e}\_p - \boldsymbol{\delta}}{\left\| \mathbf{e}\_p - \boldsymbol{\delta} \right\|} \sigma \left( \left\| \mathbf{e}\_p - \boldsymbol{\delta} \right\| \right) + \left\| \mathbf{e}\_p - \boldsymbol{\delta} \right\| \left\| \Delta \tilde{\mathbf{u}} + \vec{\sigma}\_{\varepsilon} \right\|, \end{split} \tag{A7}$$

will be ≤ 0 if:

$$\theta \lambda\_{\min}(\mathsf{K}\_{\mathcal{P}}) \sigma(\left||\mathsf{e}\_{\mathcal{P}} - \mathcal{S}\right||) \geq \left||\Delta \mathbf{\tilde{u}} + \vec{\sigma}\_{\mathcal{E}}\right||,\tag{A8}$$

which, in turn, implies that:

$$\left\|\left\|\mathbf{e}\_{p}-\mathcal{S}\right\|\right\|\geq\sigma^{-1}\left(\frac{1}{\theta\lambda\_{\min}(K\_{p})}\left\|\Delta\tilde{\mathbf{u}}+\tilde{\sigma}\_{\mathfrak{c}}\right\|\right),\tag{A9}$$

and:

$$\dot{V}\_2 \le -(1-\theta)(\mathbf{e}\_{\mathcal{P}} - \boldsymbol{\delta})^T K\_{\mathcal{P}} \sigma (\mathbf{e}\_{\mathcal{P}} - \boldsymbol{\delta}) - k\_{\mathcal{I}} |\boldsymbol{e}\_{\mathcal{I}}|^2 \tag{A10}$$

as the right side of inequality (A9) can be made arbitrarily small through the choice of the gain matrix *Kp*. It follows directly from H. Khalil ([34], Theorem 4.19) that the controlled system is ISS.

#### **Appendix B. Proof of Preposition 2**

**Proof.** Consider the candidate Lyapunov function given by:

$$V\_1(\mathbf{e}\_p) := \frac{1}{2} \mathbf{e}\_p^T \mathbf{e}\_p. \tag{A11}$$

Taking the first derivative of (A11) and replacing in (24)–(26), yields:

$$\dot{V}\_1(\mathbf{e}\_p) = -\mathbf{e}\_p^T K\_1 \mathbf{e}\_p + \mathbf{e}\_p^T \left(\mathbf{z} - \frac{\partial \mathbf{p}\_d}{\partial \gamma} \boldsymbol{e}\_\gamma\right). \tag{A12}$$

With a view to applying backstepping techniques, consider:

$$V\_2(\mathbf{e}\_{\mathcal{V}}, \mathbf{e}\_{\mathcal{V}}) := V\_1(\mathbf{e}\_{\mathcal{V}}) + \frac{1}{2} \mathbf{z}^T \mathbf{z}. \tag{A13}$$

Replacing (26), (27), and (29) in *V*˙ <sup>2</sup> yields:

$$\begin{split} \dot{V}\_{2} &= -\mathbf{e}\_{p}^{T} K\_{1} \mathbf{e}\_{p} + \mathbf{e}\_{p}^{T} \mathbf{z} - \mathbf{e}\_{p}^{T} \frac{\partial \mathbf{p}\_{d}}{\partial \gamma} \boldsymbol{e}\_{\gamma} \\ &+ \mathbf{z}^{T} \Big( \mathbf{u} + \mathbf{d} - \mathbf{h}(\gamma) (\boldsymbol{e}\_{\gamma} + \boldsymbol{v}\_{d}(\gamma\_{\prime}, t)) - \frac{\partial \mathbf{p}\_{d}}{\partial \gamma} \boldsymbol{\dot{v}}^{\rm cond} (t) + K\_{1} \mathbf{e}\_{v} - K\_{1} \frac{\partial \mathbf{p}\_{d}}{\partial \gamma} \boldsymbol{e}\_{\gamma} \Big). \end{split} \tag{A14}$$

By replacing (35)–(37) in the Lyapunov time derivative, it follows that:

$$\dot{V}\_2 = -\mathbf{e}\_p^T K\_1 \mathbf{e}\_p - \mathbf{z}^T K\_2 \mathbf{z} - \mathbf{e}\_p^T \frac{\partial \mathbf{p}\_d}{\partial \gamma} \mathbf{e}\_\gamma - \mathbf{z}^T \left(\mathbf{h}(\gamma) + K\_1 \frac{\partial \mathbf{p}\_d}{\partial \gamma}\right) \mathbf{e}\_\gamma + \mathbf{z}^T (\tilde{\mathbf{u}} + \tilde{\mathbf{d}}).\tag{A15}$$

Consider a third candidate Lyapunov, obtained by backstepping, defined as:

$$V\_3 := V\_2 + \frac{1}{2}e\_\gamma^2. \tag{A16}$$

Taking its derivative, with respect to time, and replacing in (38), we obtain:

$$\dot{\mathcal{V}}\_3 = -\mathbf{e}\_p^T K\_1 \mathbf{e}\_p - \mathbf{z}^T K\_2 \mathbf{z} - k\_\gamma e\_\gamma^2 + \mathbf{z}^T (\tilde{\mathbf{u}} + \tilde{\mathbf{d}}).\tag{A17}$$

Consider one last backstepping that involves the construction of:

$$V\_4 = V\_3 + \frac{1}{2} \mathbf{\dot{d}}^T \mathbf{K}\_d^{-1} \mathbf{\dot{d}}.\tag{A18}$$

Taking its time derivative, and taking into consideration (33):

$$\begin{split} \dot{W}\_{4} &= -\mathbf{e}\_{p}^{T} K\_{1} \mathbf{e}\_{p} - \mathbf{z}^{T} K\_{2} \mathbf{z} - k\_{\gamma} \mathbf{e}\_{\gamma}^{2} + \underbrace{\dot{\mathbf{d}}^{T} (\mathbf{z} - \text{Proj}(\mathbf{z}, \dot{\mathbf{d}}))}\_{\leq 0} + \mathbf{z}^{T} \tilde{\mathbf{u}} \\ &\leq -\mathcal{W}(\mathbf{e}\_{p}, \mathbf{e}\_{\upsilon}, \mathbf{e}\_{\gamma}) + \mathbf{z}^{T} \tilde{\mathbf{u}}. \end{split} \tag{A19}$$

Assuming that the quadrotor is equipped with a generic inner loop that is capable of keeping the tracking error **u˜** small and bounded, the right side of inequality (A19) can be made small enough such that the controlled system is stable. A more in-depth stability analysis can be conducted for the inner–outer loop control system, but this will be dependent directly on the type of inner loop adopted. This results from the fact that the desired accelerations **ud** must be decoupled in a set of desired thrusts and attitudes for the quadrotor to track.

In order to simplify the designed control law **u***d*, consider the final algebraic manipulation:

$$\begin{split} \mathbf{u}\_{d}^{\diamond} &= -\mathbf{\hat{d}} + \mathbf{h}(\gamma)v\_{d}(\gamma, t) + \frac{\partial \mathbf{p}\_{d}}{\partial \gamma} \dot{v}^{\text{corord}}(t) - K\_{1} \mathbf{e}\_{v} - \mathbf{e}\_{p} - K\_{2} \mathbf{z} \\ &= -\mathbf{\hat{d}} + \mathbf{h}(\gamma)v\_{d}(\gamma, t) + \frac{\partial \mathbf{p}\_{d}}{\partial \gamma} \dot{v}^{\text{corord}}(t) - \mathbf{e}\_{v} \underbrace{(K\_{1} + K\_{2})}\_{K\_{v}} - \mathbf{e}\_{p} \underbrace{(I + K\_{1}K\_{2})}\_{K\_{p}}. \end{split} \tag{A20}$$

#### **Appendix C. Proof of Preposition 3**

**Proof.** Consider that:

$$\mathbf{v}\_{\mathbf{L}}(\gamma) = v\_L \mathbf{1} + \tilde{\mathbf{v}}\_{\mathbf{L}\prime} \tag{A21}$$

where **v˜ <sup>L</sup>** is a bounded and arbitrarily small term that accounts for a transient period in which the vehicles are on different sections of the path, with slightly different desired speed profiles. Replacing (39), the speed correction term proposed in (46), and (A21) in (42) yields:

$$\begin{split} \dot{\varepsilon} &= L(\mathbf{v}\_{\mathcal{L}}(\gamma) - k\_{\varepsilon}(\varepsilon + \mathcal{A}\tilde{\gamma})) \\ &= v\_{\mathcal{L}}L\mathbf{\mathcal{A}} + L\tilde{\mathbf{v}}\_{\mathcal{L}} - k\_{\varepsilon}L(\varepsilon + \mathcal{A}\tilde{\gamma}) \\ &= -k\_{\varepsilon}L(\varepsilon + \mathbf{d}) \text{ with } \mathbf{d} = \frac{\tilde{\mathbf{v}}\_{\mathcal{L}}}{k\_{\varepsilon}} + \mathcal{A}\tilde{\gamma} \end{split} \tag{A22}$$

where **d** is a disturbance that results from combining the terms dependent on **v˜ <sup>L</sup>** and *γ***˜**. Consider the Laplacian matrix *L*, expressed in canonical Jordan form as:

$$L = V \Lambda V^{-1} \, , \tag{A23}$$

and the change of variables:

$$
\overline{\mathfrak{e}} = V^{-1} \mathfrak{e}.\tag{A24}
$$

Applying (A24) to (A22) yields:

$$\mathfrak{E} = -k\_{\mathfrak{e}} \Lambda(\mathfrak{E} + \mathfrak{A}) \text{, with } \mathfrak{A} = V^{-1} \mathfrak{A}. \tag{A25}$$

It is possible to decompose the above equality according to the notation:

$$
\begin{bmatrix}
\frac{\dot{\varepsilon}\_1}{\dot{\varepsilon}\_2}
\end{bmatrix} = \begin{bmatrix}
0\\-k\_\ell \Lambda\_2 (\mathfrak{E}\_2 + \overleftarrow{\mathbf{d}}\_2)
\end{bmatrix},
\tag{A26}
$$

where the first half of the vector denotes the term that depends on the null eigenvalue of the Laplacian, while the second term is a vector that depends only on the positive eigenvalues of the Laplacian. Consider now the candidate Lyapunov function:

$$V\_{\mathbb{P}\_2} = \frac{1}{2} \mathbb{E}\_2^T \mathbb{E}\_{2\prime} \tag{A27}$$

and its time derivative, given by:

$$\begin{split} \dot{V}\_{\mathbb{E}\_2} &= -k\_{\mathbb{E}} \mathbb{E}\_2^T \Lambda\_2 (\mathbb{E}\_2 + \mathbf{\tilde{d}}\_2) \\ &= -(1 - \theta) k\_{\mathbb{E}} \mathbb{E}\_2^T \Lambda\_2 \mathbb{E}\_2 - \theta k\_{\mathbb{E}} \mathbb{E}\_2^T \Lambda\_2 \mathbb{E}\_2 - k\_{\mathbb{E}} \mathbb{E}\_2^T \Lambda\_2 \mathbf{\tilde{d}}\_2, \end{split} \tag{A28}$$

where 0 < *θ* < 1. The term:

$$-\theta k\_{\ell} \mathbb{E}\_2^T \Lambda\_2 \mathbb{E}\_2 - k\_{\ell} \mathbb{E}\_2^T \Lambda\_2 \overline{\mathbf{d}}\_{2,\prime} \tag{A29}$$

will be ≤ 0 if:

$$\|\|\mathbb{E}\_2\|\| \ge \frac{1}{\theta} \|\|\overline{\mathbf{d}}\_2\|\|\_{\prime} \tag{A30}$$

and, therefore:

$$
\dot{V}\_{\mathbb{Z}\_2} \le -(1-\theta)k\_{\varepsilon}\mathbb{E}\_2^T \Lambda\_2 \mathbb{1}\_2. \tag{A31}
$$

The term *γ***˜** can be made arbitrarily small by controlling the gains that dictate the broadcasting scheme. Moreover, the term **v˜ <sup>L</sup>** can be dominated by a proper choice *kε*. Hence, **<sup>d</sup>** can be made arbitrarily small and, thus, # #**d¯** 2 # # can be as well. It follows directly from H. Khalil ([34], Theorem 4.19) that the controlled system is ISS with respect to the error vector *ε* and the inputs *γ***˜** and **v˜ <sup>L</sup>**.

#### **Appendix D. Computing the Regularisation Term Using Vectorial Notation**

Consider the simplest unclamped uniform cubic B-spline with only one segment, such that *γ* ∈ [0, 1] and is described by (6). Then, its first derivative **C** (*γ*) is given by:

$$
\underbrace{\frac{\partial \mathbb{C}^{x/y}}{\partial \boldsymbol{\gamma}}(\boldsymbol{\gamma})}\_{\mathbf{T}(\boldsymbol{\gamma})}(\boldsymbol{\gamma}) = \underbrace{\underbrace{\begin{bmatrix} \boldsymbol{\gamma}^{2} & \boldsymbol{\gamma} & 1 & 0 \end{bmatrix}}\_{\mathbf{T}(\boldsymbol{\gamma})} \underbrace{\frac{1}{6} \begin{bmatrix} 3 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} -1 & 3 & -3 & 1 \\ 3 & -6 & 3 & 0 \\ -3 & 0 & 3 & 0 \\ 1 & 4 & 1 & 0 \end{bmatrix} \begin{bmatrix} P\_{0}^{x/y} \\ P\_{1}^{x/y} \\ P\_{2}^{x/y} \\ P\_{3}^{x/y} \\ P\_{4}^{x/y} \end{bmatrix}}\_{\mathbf{P}}.\tag{A32}
$$

Therefore, the term *<sup>γ</sup>* **C** (*γ*)<sup>2</sup> *dγ* is computed according to:

$$\begin{split} \int\_{\gamma} \left\| \mathbf{C}'(\gamma) \right\|^2 d\gamma &= \int\_{\gamma} (\mathbf{B}'(\gamma)\mathbf{P})^T (\mathbf{B}'(\gamma)\mathbf{P}) d\gamma \\ &= \int\_0^1 \mathbf{P}^T B'(\gamma)^T B'(\gamma) \mathbf{P} d\gamma \\ &= \mathbf{P}^T M^T \left[ \int\_0^1 \mathbf{T}(\gamma)^T \mathbf{T}(\gamma) d\gamma \right] \text{MP.} \end{split} \tag{A33}$$

Further, note that:

$$\mathbf{T}(\gamma)^T \mathbf{T}(\gamma) = \begin{bmatrix} \gamma^4 & \gamma^3 & \gamma^2 & 0\\ \gamma^3 & \gamma^2 & \gamma & 0\\ \gamma^2 & \gamma & 1 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix} \tag{A34}$$

and, as a consequence:

$$\int\_0^1 \mathbf{T}(\gamma)^T \mathbf{T}(\gamma) d\gamma = Q = \begin{bmatrix} 1/5 & 1/4 & 1/3 & 0 \\ 1/4 & 1/3 & 1/2 & 0 \\ 1/3 & 1/2 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \tag{A35}$$

Hence, for the simplest case of a single B-spline segment, it is known that:

$$\int\_{\gamma} \left\| \mathbf{C}'(\gamma) \right\|^2 d\gamma = \mathbf{P}^T \underbrace{M^T Q M}\_{R\_1} \mathbf{P}.\tag{A36}$$

The easiest way to extend this technique to a B-spline with *n* segments is to consider the modified vector **<sup>T</sup>**(*γ*) = [(*<sup>γ</sup>* <sup>−</sup> *<sup>i</sup>*)2,(*<sup>γ</sup>* <sup>−</sup> *<sup>i</sup>*), 1, 0] *<sup>T</sup>*, where *<sup>i</sup>* <sup>=</sup> *<sup>γ</sup>*, according to the notation introduced in Section 2.3. Then, since (*γ* − *i*) ∈ [0, 1], one can compute individually, for each segment, intermediate matrices *R<sup>i</sup>* <sup>1</sup>, according to (A36). Due to the locality property of B-splines, it is possible to "stack" these intermediate matrices to form the final matrix *R*1. An analogous rationale can be applied to compute *R*2.
