A crosswind perturbed PVTOL aircraft robust global control solution for the trajectory tracking problem

Abstract

This study solves the trajectory-tracking control problem for a crosswind-perturbed crewless PVTOL aircraft. The proposed globally feedback-based solution consists of two parts: one makes the aircraft’s horizontal and vertical position accelerations convenient, ensuring no singularities in the closed-loop system. To accomplish the control task, two controllers were designed: a continuous controller to ensure the aircraft tracks the trajectory, and a discontinuous integral-adaptive sliding-mode controller to counteract the destabilizing crosswind effect. The first controller was designed using a simple backstepping approach, with a saturation function replacing linear control to avoid singular points. The second controller was designed using Lyapunov functions and Barbalat’s lemma. The solution’s performance was tested through numerical simulations, yielding convincing results.

Introduction

Even when crewless aircraft have very complex nonlinear dynamics, we need to control them to make them carry out at least some easy maneuvers, for instance, reaching a point and staying stable. However, designing such a controller is a challenging task. Fortunately, a few decades ago, Hauser et al. ¹ introduced the simplified Planar Vertical Take Off and Landing (PVTOL) aircraft to model, as its name implies, the take off and landing of an airplane restricted to move in a vertical plane. This simplified model embodies several characteristics and behaviors found in quadrotors, helicopters, or the Harrier Jump Jet,^2–4 alleviating the complexities found in actual aircraft to some extent. Also, the PVTOL is an underactuated system; it has two controllers and three degrees of freedom. All these features make it an ad-hoc benchmark for testing control techniques.

In general, regulation and trajectory tracking are the main problems that have been addressed by the control community when designing control laws for the PVTOL aircraft, with Hauser et al. ¹ and Martin et al. ⁵ being the pioneering works. Related to the trajectory tracking problem, using the feedback linearization approach, a controller to stabilize a non-minimum phase approximation of the PVTOL system is designed in Hauser et al. ¹ In Martin et al., ⁵ the authors propose a nonlinear control to accomplish the trajectory tracking problem for a version of the PVTOL. To this end, the authors took advantage of the aircraft’s flatness property, ensuring bounded zero dynamics. The latter is thanks to applying the inverse trajectory feedforward and a state tracker. In this line, in Fliess et al., ⁶ the authors combined the PVTOL flatness property and trajectory planning using a Bézier spline to accomplish smooth translational maneuvers. Once again, exploiting the PVTOL flatness property, the authors of Sira-Ramfrez ⁷ designed an iterative off-line algorithm to approximately estimate the flat output reference trajectories in terms of the desired non-minimum phase output trajectories. Recently, in Bonilla et al., ⁸ based on implicit systems techniques, the authors introduced a novel exact linearization procedure applied to control a PVTOL aircraft. The idea behind this technique is to split the nonlinear state representation into a basic rectangular representation, that is, the system is linearized, and an auxiliary nonlinear algebraic equation. This representation can be applied to more sophisticated nonlinear dynamic models. Using an input-output feedback linearization technique, a linear controller is proposed in Bonilla et al., ⁹ provided state availability. Then, a tracking control law was obtained using pole placement and regular perturbation techniques, which was numerically evaluated using a PVTOL model. We mentioned that above, briefly reviewed studies only ensure local stability, and the corresponding control laws are considerably more complex than the ones proposed to solve regulation. A class of sliding-mode-based controllers, known as integral and terminal formulations, for the twin-rotor platform is introduced in Irfan et al. ¹⁰ This approach systematically designs robust control mechanisms to reduce chattering and address the nonlinearities and strong dynamic coupling characteristic of unmanned aerial vehicle systems. In Saleem et al., ¹¹ the authors proposed a novel robust position tracking controller for Vertical Take Off and Landing Unmanned Aerial Vehicles based on a Linear Quadratic Integral (LQI) framework enhanced with a nonlinear rate-varying integral (RVI) compensator. Their approach replaces the conventional integral term in the LQI law with an adaptive mechanism that dynamically scales according to both the tracking-error velocity and the braking acceleration of the reference trajectory. The controller was validated experimentally on a hardware-in-the-loop aero-pendulum testbed equipped with coaxial contra-rotating propellers.

On the regulation side, a backstepping-based approach is developed in Sepulchre et al. ¹² to control the translational dynamics using the cascaded roll subsystem. In Fantoni and Lozano, ¹³ the authors developed an algorithm based on a forward control approach to ensure asymptotic stability using a Lyapunov function. The ideas from this study were used^14–18 to develop solutions for real-world PVTOL applications. For instance, in Zavala-Río et al., ¹⁷ nonlinear combinations of saturation functions were used to limit thrust input to arbitrary saturation limits, thereby developing a global stabilizer controller for a PVTOL aircraft. The Lyapunov theory and saturation functions are used in Sanchez et al. ¹⁹ to propose a global stabilizer control for PVTOL aircraft, which was tested using real-time experiments. A stabilizing control strategy based on Robust Control Lyapunov Functions and Sontag’s universal stabilizing feedback for a crosswind perturbed PVTOL is introduced in Munoz et al. ²⁰ In Castillo et al., ²¹ the authors designed an algorithm that uses the robust control Lyapunov function approach and unknown disturbances rejection for stabilizing a wind-perturbed vertical take-off and landing aircraft. Using two simplified unmanned PVTOL aircraft, the authors of Escobar et al. ²² solved the slung-load transportation problem. To this end, a decentralized passivity-based control scheme was implemented, where the PVTOL’s physical connection to the load served as a means of communication. Based on a feedback linearized PVTOL model that does not cover the entire state space due to a singularity, the authors of Escobar et al. ²³ developed novel controllers that can be adapted to other vehicles. This study uses a feedback linearization-based control that avoids reaching any singularity; the necessary conditions for ensuring asymptotic stability are given. Considering a PVTOL system with a constant force applied to the horizontal axis, a simplified stabilizing control approach was developed in Lozano et al. ²⁴ Under some suitable conditions, this approach ensures exponential stability. In Villaseñor Rios et al., ²⁵ a different extended linear state observer and active disturbance rejection control-based approach is developed to control a nonlinear underactuated PVTOL aircraft system with an inverted pendular load. Using the ideas found in Teel’s paper, ²⁶ the authors of Lozano et al. ²⁷ introduce a modified nested saturation control to stabilize the model of the PVTOL system.

In our humble opinion, the above briefly reviewed PVTOL aircraft control-related studies are the closest to the solution we present in this manuscript; we suggest the interested reader consult the work, ²⁸ where a well-developed and structured evaluation of popular control algorithms for the PVTOL can be found.

In this work, we propose a globally feedback-based solution for trajectory tracking control of a crewless PVTOL aircraft subject to undesired crosswind effects. To this end, we applied a backstepping version with a saturation function. The idea is to force the horizontal and vertical aircraft dynamics to predefined suitable acceleration references in their respective directions. These accelerations allow the PVTOL to avoid the singular angular position. Afterward, we derive two control laws: a continuous one for trajectory tracking and an integral adaptive sliding mode controller to counteract the unknown bounded crosswind effect. To design the continuous controller, we do not use the feedback linearization method or the differential flatness approach; instead, we propose using an implicit Immersion and Invariance method-based controller. To assess the effectiveness of our control approach, we designed and implemented three numerical simulations, one of which compares against two well-established related control solutions. In every case, we obtained satisfactory results. Even when the regulation problem has been exhaustively studied, the singularities have been neglected, opposite to what we do in this work, in which we give them the corresponding protagonism, because in actual applications, they cannot be ignored.

We organized the rest of our study as follows. We introduced the PVTOL dynamical model, the problem statement, and the control strategy in Section 2. We designed and implemented the numerical simulations needed to assess the effectiveness of our control approach in Section 3. Finally, we presented the concluding remarks in Section 4.

• Notation: In what follows, we use the following notation:

• The symbol $\mathrm{sgn}[x]$ , with $x\in \mathfrak{R}$ , refers to the signum function of a real number and is defined as:

$$\mathrm{sgn}[x]=\{1\mathrm{if}x>0;-1\mathrm{if}x<0;[-1,1]\mathrm{if}x=0\}.$$

• The symbol $\mathrm{sgn}[y]$ , with $y\in {\in \mathfrak{R}}^n$ , refers to the signum function of a vector of real numbers and is defined as:

$$\mathrm{sgn}[y]={\left[\begin{array}{ccc}\multicolumn{1}{c}{\mathrm{sgn}[y_1]}& \dots \hfill & \mathrm{sgn}[y_n]\hfill \end{array}\right]}^T.$$

Dynamical model

The PVTOL aircraft is an underactuated system with two control inputs and three degrees of freedom. This aircraft operates within a two-dimensional plane and is equipped with two independent motors that generate both force and torque on the vehicle (see Figure 1). The dynamics of the PVTOL aircraft are described by the following model^20,29:

$$\begin{array}{c}\multicolumn{1}{c}{m\stackrel{··}{X}=-\sin \theta (f_{1_T}+f_{2_T})+\epsilon \cos \theta (f_{1_T}-f_{2_T})L;}\\ \multicolumn{1}{c}{m\stackrel{··}{Y}=\cos \theta (f_{1_T}+f_{2_T})+\epsilon \sin \theta (f_{1_T}-f_{2_T})L-\mathrm{mg};}\\ \multicolumn{1}{c}{J\stackrel{··}{\theta }=(f_{1_T}-f_{2_T})L;}\end{array}$$

(1)

where $X$ and $Y$ represent the positional coordinates of the aircraft’s center of mass within a fixed inertial reference frame, while $\theta$ denotes the aircraft’s inclination relative to the horizontal axis, $m$ is the total mass of the aircraft, $g$ is the gravitational acceleration, $L$ is the distance between the rotor and the aircraft’s center of mass, and $J$ is the moment of inertia. We denote the forces applied to the aircraft by $f_{1_T}=f_1+f_{w_1}$ and $f_{2_T}=f_2+f_{w_2}$ , where $f_1$ and $f_2$ are the thrust the rotors produce. At the same time, $f_{w_1}$ and $f_{w_2}$ are the crosswind forces (see for instance, Munoz et al. ²⁰ ), we can find a detailed treatment of the crosswind effect). We denote the coupling between the rolling moment and the lateral acceleration by $\epsilon$ .

• Remark 1. From a practical point of view, $\epsilon$ is not very significant because physically, $\epsilon <<1$ . Besides, the parameter $\epsilon$ is difficult to estimate.

Figure 1.

The PVTOL aircraft.

To facilitate the subsequent developments, we normalized the system described in equation (1) by using the following variable transformations:

$$\begin{array}{c}x=\frac{X}{g};y=\frac{Y}{g};\hfill \\ u_{T_1}=\frac{f_1+f_2}{\mathrm{mg}};u_{T_2}=\frac{(f_1-f_2)L}{J};\hfill \\ w_1=\frac{f_{w_1}+f_{w_2}}{\mathrm{mg}};w_2=\frac{(f_{w_1}-f_{w_2})L}{J},\hfill \end{array}$$

and considering $\epsilon =0$ it is easy to show that the model (1) can be written as:

$$\begin{array}{c}\stackrel{··}{x}=-\sin \theta (u_{T_1}+w_1);\hfill \\ \stackrel{··}{y}=\cos \theta (u_{T_1}+w_1)-1;\hfill \\ \stackrel{··}{\theta }=u_{T_2}+w_2.\hfill \end{array}$$

(2)

Now, we introduce the control goal of this study.

• Control Problem: Let us consider the PVTOL model (2) under the crosswind effect and the following assumptions:

• A1: The whole system state is measurable.

• A2: The crosswind forces are bounded as $\left|w_i\right|<\overline{w}$ , for $i=\{1,2\}$ , for some $\overline{w}>0$ .

• Then, the control objective consists of designing a controller such that the tracking errors:

$$\begin{array}{ccc}\multicolumn{1}{c}{e_{T_x}=x-x_r}& \mathrm{and}\hfill & e_{T_y}=y-y_r\hfill \end{array}$$

(3)

• fulfill:

$$\begin{array}{c}\multicolumn{1}{c}{\underset{t\to \mathrm{\infty }}{\lim }\left|e_{T_x}\left(t\right)\right|=0,\underset{t\to \mathrm{\infty }}{\lim }\left|e_{T_y}\left(t\right)\right|=0,}\end{array}$$

while $\theta$ is bounded. The terms $x_r$ and $y_r$ are smooth, bounded desired references with bounded time derivatives up to the fourth order.

Limitations

• The proposed solution requires system position and velocity measurements, as established by assumption A1, and parameters $m$ , $L$ , and $J$ must be known.

• In a context where the system velocity is not available, we can overcome this inconvenience using, for instance, a super-twisting observer like the one proposed in Davila et al. ³⁰ and Cardenas and Aguilar, ³¹ the generalized proportional integral (GPI) observer introduced in Sira-Ramírez, ³² or algebraic observation methods like the ones introduced in Cortés-Romero et al. ³³ and Sira-Ramírez et al. ³⁴ All these alternatives have finite time convergence, and the latter relaxes assumption A1.

We finish this section by expressing the control model (2) in its compact form as:

$$\stackrel{·}{\chi }=F\left(\chi \right)+B\left(\theta \right)(u_T+w),$$

(4)

where $\chi =[q^T,p^T]\in D\subset {\in \mathfrak{R}}^6$ and:

$$\begin{array}{c}\begin{array}{ccc}\multicolumn{1}{c}{q=\left[\begin{array}{c}x\hfill \\ \multicolumn{1}{c}{y}\\ \multicolumn{1}{c}{\theta }\end{array}\right]}& p=\left[\begin{array}{c}\stackrel{·}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{y}}\\ \multicolumn{1}{c}{\stackrel{·}{\theta }}\end{array}\right]\hfill & F\left(\chi \right)=\left[\begin{array}{c}p\hfill \\ \multicolumn{1}{c}{-e_2}\end{array}\right]\hfill \end{array}\hfill \\ \multicolumn{1}{c}{B\left(\theta \right)={\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & -S_{\theta }\hfill & C_{\theta }\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]}^T,}\end{array}$$

(5)

where $e_2=[0,1,0{]}^T$ , $S_{\theta }=\sin \left(\theta \right)$ , and $C_{\theta }=\cos \left(\theta \right)$ . Please observe that $F\left(\chi \right)\in {\in \mathfrak{R}}^{6\times 1}$ , $B\left(\theta \right)\in {\in \mathfrak{R}}^{6\times 2}$ , $u_T$ , $w\in {\in \mathfrak{R}}^{2\times 1}$ and $\mathrm{rank}\left(B\left(\theta \right)\right)=2$ for all $\theta \in \mathfrak{R}$ . To solve the control problem subject to matching perturbations, we split the controller $u_T$ as follows:

$$u_T=u_N+u_I=\left[\begin{array}{c}{u}_{N_1}\hfill \\ \multicolumn{1}{c}{u_{N_2}}\end{array}\right]+\left[\begin{array}{c}{u}_{I_1}\hfill \\ \multicolumn{1}{c}{u_{I_2}}\end{array}\right],$$

(6)

where the nominal controller $u_N$ solves the trajectory tracking problem, and the integral controller $u_I$ prevents the undesirable effects of the endogenous disturbances.

Challenges posed by the PVTOL system

• We must point out that in the case of $\epsilon \ne 0$ and $f_{w_1}=f_{w2}=0$ , we can straightforwardly demonstrate, using certain nonlinear transformations with respect to the inputs $f_{1_T}$ and $f_{2_T}$ , that the internal dynamics of the original system (1) is given by:

$$ϵ\stackrel{··}{\theta }=k\sin \theta$$

with $k>0$ constant. Evidently, the above equation describes the dynamics of an undamped pendulum, for which it is well known that its two equilibrium points are not asymptotically stable. In conclusion, the system is a nonminimum phase system. ³⁵ Consequently, the system may exhibit an undershoot or an initial response opposite to the desired direction. As we mentioned before, the main control goal of our approach is avoiding these singular points.

• If we desire to solve the trajectory tracking problem when $w=0$ and $\epsilon =0$ , based on feedback linearization or flatness-based control approaches, we need to transform the system (2) into the following configuration ²³ :

$$\left[\begin{array}{c}{x}^{\left(4\right)}\hfill \\ \multicolumn{1}{c}{y^{\left(4\right)}}\end{array}\right]=M(u_1,\theta )\left[\begin{array}{c}{\stackrel{··}{u}}_{T_1}\hfill \\ \multicolumn{1}{c}{u_{T_2}}\end{array}\right]+O(\theta ,\stackrel{·}{\theta },u_T,{\stackrel{·}{u}}_T)$$

(7)

where $\det \left(M\right)=u_{T_1}0>$ and $O(*)$ is a suitable nonlinear term. Under this condition, the dynamic control input $u_{T_1}$ is required to maintain a constant sign to allow the implementation of a second-order stabilizing controller $u_{T_1}$ . Besides, it is recommended to select the initial conditions $u_{T_1}\left(0\right)>0$ and ${\stackrel{·}{u}}_{T_1}\left(0\right)\ge 0$ to avoid singularities of the system (7).

itemize

• Comment 1. If the parameter $\epsilon$ is known with high accuracy, we can introduce a convenient change of coordinates as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\overline{x}=x-\epsilon \cos \theta }\\ \multicolumn{1}{c}{\overline{y}=y-\epsilon (\cos \theta -1)}\end{array}$$

which transforms the system (2) in

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{\overline{x}}=-\sin \theta {\overline{u}}_{T_1};}\\ \multicolumn{1}{c}{\stackrel{··}{\overline{y}}=\cos \theta {\overline{u}}_{T_1}-1;}\\ \multicolumn{1}{c}{\stackrel{··}{\theta }=u_{T_2}.}\end{array}$$

where ${\overline{u}}_{T_1}=u_{T_1}+\epsilon {\stackrel{·}{\theta }}^2$ , and $w=0$ . That is, we once again obtain our control model (2). This provides the justification for omitting the parameter $\epsilon$ in our control model.

Control strategy

In the first stage, our control strategy proposes a nominal control law $u_N$ to solve the trajectory tracking problem for the unperturbed system (4); that is, we assume $w=0$ . Afterward, we consider the case where $w\ne 0$ . To this end, we will apply the integral slide mode control approach, found in Castanos and Fridman, ³⁶ to design the controller $u_I$ that aims to counteract the undesirable effect of the bounded and matched perturbation $w$ , which represents the crosswind in a real scenario. ²¹ We end this section by introducing the following helpful definition:

• Definition 1: we use ${\sigma }_m[x]:\mathfrak{R}\to \mathfrak{R}$ to refer to a nonlinear saturation function, having the property that ${\sigma }_m[0]=0$ ; ${\sigma }_m[x]x\ge 0$ ; $|{\sigma }_m[x]|\le m$ , and ${\sigma }_m\prime [x]\ge 0$ for all $x\in \mathfrak{R}$ with $m>0$ .

According to this definition, the following properties are fulfilled:

• ${\mathbf{P}}_1$ ) ${\int }_0^x{s}_m[s]\mathrm{ds}>0$ ; for any $x\ne 0$ ;

• ${\mathbf{P}}_2$ ) $|s_m[x+y]-s_m[y]|\le 2\mathrm{m}\mathrm{ax}\{1,m\}\left|y\right|$ ;

• ${\mathbf{P}}_3$ ) $(s_m[x]-s_m[x+y])y\le 0$ , for any $x$ , $y\in \mathfrak{R}$ .

Properties ${\mathbf{P}}_1$ and ${\mathbf{P}}_2$ can be found in Teel, ³⁷ while property ${\mathbf{P}}_3$ is a direct consequence of a simple computation. ³⁸ For instance, we can use ${\sigma }_a[x]=a\arctan \left(x\right)$ , where $a>0$ can be any constant. In the following developments, we use the saturation function.

Nominal control design

Following the ideas found in Fantoni and Lozano ¹³ and Cardenas and Aguilar, ³¹ we solve the trajectory tracking problem by using a simple backstepping approach. However, we use saturation functions instead of simple linear control to avoid the singular points $\theta =\pm \pi /2$ .

First of all, we must note that, from relation (2), we can obtain the following expression:

$$\tan \theta =-\frac{\stackrel{··}{x}}{\stackrel{··}{y}+1},$$

(8)

and

$$u_{N_1}=-\stackrel{··}{x}\sin \theta +(\stackrel{··}{y}+1)\theta .$$

(9)

The idea behind our control proposal consists of designing the nominal control $u_{N_1}$ , such that the accelerations in the directions $x$ and $y$ have the following structure:

$$\begin{array}{c}{r}_x(x,t)=-{\sigma }_a[k_p(x-x_r)+k_d(\stackrel{·}{x}-{\stackrel{·}{x}}_r)]+{\stackrel{··}{x}}_r,\hfill \\ \multicolumn{1}{c}{r_y(x,t)=-{\sigma }_b[k_p(x-x_r)+k_d(\stackrel{·}{x}-{\stackrel{·}{x}}_r)]+{\stackrel{··}{y}}_r,}\end{array}$$

(10)

where the positive constants $\{k_p,k_d,a,b\}$ are appropriately selected, such that:

$$b+{\stackrel{··}{y}}_r\left(t\right)<1$$

(11)

holds. In other words, we desire that accelerations in $x$ and $y$ have the following structure:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{x}=-{\sigma }_a[k_p(x-x_r)+k_d(\stackrel{·}{x}-{\stackrel{·}{x}}_r)]+{\stackrel{··}{x}}_r,,}\\ \multicolumn{1}{c}{\stackrel{··}{y}=-{\sigma }_b[k_p(x-x_r)+k_d(\stackrel{·}{x}-{\stackrel{·}{x}}_r)]+{\stackrel{··}{y}}_r.}\end{array}$$

(12)

If the expressions above hold, then the tracking errors (3) also fulfill:

$$\begin{array}{c}{\stackrel{··}{e}}_{T_x}=-{\sigma }_a[k_p{e}_{T_x}+k_d{\stackrel{·}{e}}_{T_x}],\hfill \\ \multicolumn{1}{c}{{\stackrel{··}{e}}_{T_y}=-{\sigma }_b[k_p{e}_{T_y}+k_d{\stackrel{·}{e}}_{T_y}].}\end{array}$$

(13)

We underscore that the equation (13) ensures globally asymptotic and locally exponential convergence to zero (see Supplemental Appendix A), provided that conditions of relations (11) and (12) are satisfied. The relations (8) and (9) suggest us to select $u_{N_1}$ and $\tan \theta$ , such that the following holds:

$$\tan \theta =-\frac{r_x(x,\stackrel{·}{x},t)}{r_y(y,\stackrel{·}{y},t)+1},$$

(14)

and:

$$u_{N_1}=-r_x(x,\stackrel{·}{x},t)\sin \theta +(r_y(y,\stackrel{·}{y},t)+1)\theta ,$$

(15)

where, according to relation (11), $r_y(y,\stackrel{·}{y},t)+1\ne 0$ for all $t$ , helping us to avoid the singular points $\theta =\pm \pi /2$ .

To control the angle $\theta$ and ensuring that condition (14) holds, we design the input $u_{N_2}$ , such that the following error:

$$\xi =\theta +\mathrm{arc}\mathrm{t}\mathrm{an}\left(\frac{r_x(x,\stackrel{·}{x},t)}{r_y(y,\stackrel{·}{y},t)+1}\right),$$

(16)

goes to zero. If $\xi =0$ , then the restriction (14) automatically holds. To this end, we propose the following auxiliary variable:

$$r=\frac{r_x(x,\stackrel{·}{x},t)}{r_y(y,\stackrel{·}{y},t)+1}.$$

(17)

Obtaining the second time derivative of (16) and using the third equation of (2), we have that:

$$\stackrel{··}{\xi }=u_{N_2}-\frac{2{\mathrm{rr}}^{\prime }}{{(1+r^2)}^2}+\frac{r^{″}}{1+r^2}.$$

(18)

Therefore, to ensure that $\xi$ exponentially converge to zero, we propose the controller $u_{N_2}$ as:

$$u_{N_2}=-k_0\xi -k_1\stackrel{·}{\xi }+\frac{2{\mathrm{rr}}^{\prime }}{{(1+r^2)}^2}+\frac{r^{″}}{1+r^2},$$

(19)

where $k_0$ and $k_1$ are strictly positive constants and ¹ (To simplify, we stand for $r_x=r_x(x,\stackrel{·}{x},t)$  and $r_y=r_y(y,\stackrel{·}{y},t)$ .):

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\xi }=\stackrel{·}{\theta }+\frac{r^{\prime }}{1+r^2};}\\ \multicolumn{1}{c}{r^{\prime }=-\frac{r_x{r}_y^{\prime }}{{(1+r_y)}^2}+\frac{r_x^{\prime }}{1+r_y};}\\ \multicolumn{1}{c}{r^{″}=-\frac{2r_x^{\prime }{r}_y^{\prime }}{{(1+r_y)}^2}+\frac{2r_x{\left(r_y^{\prime }\right)}^2}{{(1+r_y)}^3}+\frac{r_x^{″}}{1+r_y}-\frac{r_y^{″}}{{(1+r_y)}^2}.}\end{array}$$

The algebraic expression of the set of variables:

$$\begin{array}{cccccc}\multicolumn{1}{c}{\{r_x^{\prime },}& r_x^{″},\hfill & r_y^{\prime },\hfill & r_y^{″},\hfill & r\prime ,\hfill & r″\}.\hfill \end{array}$$

(20)

can be found in Supplemental Appendix B.

We must note that our control law (15)–(19) avoids singular points. It suggests that our control law ensures global stability in the trajectory tracking for any initial condition, as long as condition (11) holds. Additionally, implementing the controller (19) and the algebraic variables in (20) can be a tedious and error-prone task. To alleviate this issue, we can use the following tracking differentiator^39–41:

$$\begin{array}{c}{\stackrel{·}{\vartheta }}_1={\vartheta }_2\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\vartheta }}_2={\vartheta }_3}\\ \multicolumn{1}{c}{{\stackrel{·}{\vartheta }}_3=-R^3[w_n^3({\vartheta }_1-z)+3\frac{w_n^2{\vartheta }_2}{R}+3\frac{w_n{\vartheta }_3}{R^2}].}\end{array}$$

(21)

For $R_n\to \mathrm{\infty }$ and $w_n>1$ , the following holds as long as $\left|z^k\right|<M$ :

$$\begin{array}{ccc}\multicolumn{1}{c}{{\vartheta }_1=z}& {\vartheta }_2=\stackrel{·}{z}\hfill & {\vartheta }_3=\stackrel{··}{z}.\hfill \end{array}$$

(22)

In other words, we do not need to compute the algebraic expression (20). Particularly, we used the above-introduced tracking derivator (21) to estimate the signals in (20), with the filter values $R=75$ and $w_n=1$ . To counteract the peaking phenomena produced by the system above, we used the clutch function proposed in equation (17) of Luviano-Juárez et al., ⁴² where we instantiate its parameters as ${ϵ}_i=0.5$  s and $q=4$ .

Active compensation of matching perturbation

Inspired by Castanos and Fridman ³⁶ and Rubagotti et al., ⁴³ we address the crosswind disturbance rejection assuming $\overline{w}$ is unknown. To this end, we design this control goal by introducing the following convenient sliding manifold

$${\mathrm{\Delta }}_{\chi }(t,t_0)={\chi }_{\mathrm{\Delta }}-{\int }_{t_0}^t(F\left(\chi \left(\tau \right)\right)+B\left(\theta \left(\tau \right)\right)u_N\left(\tau \right))d\tau ,$$

(24)

where ${\chi }_{\mathrm{\Delta }}=\chi \left(t\right)-\chi \left(t_0\right)$ and (24) is the difference between the actual and the unperturbed ideal solutions. Notice that $s(\chi ,t)\in {\in \mathfrak{R}}^2$ and $\mathrm{rankB}\left(\theta \right)=2$ . Then, computing the time derivative of (23), we obtain:

$$\stackrel{·}{s}(\chi ,t)=\stackrel{·}{\theta }\frac{{\mathrm{dB}}^T\left(\theta \right)}{d\theta }{\mathrm{\Delta }}_{\chi }(t,t_0)+B^T(\stackrel{·}{\chi }-F\left(\chi \right)-B\left(\theta \right)u_N).$$

Substituting the value of:

$$\stackrel{·}{\chi }=F\left(\chi \right)+B\left(\theta \right)(u_T+w),$$

given in (4), considering that $B^T\left(\theta \right)B\left(\theta \right)=I_2$ , and using the definition in (6), we obtain after some algebraic manipulation:

$$\stackrel{·}{s}(\chi ,t)=u_I+w+\stackrel{·}{\theta }\frac{{\mathrm{dB}}^T\left(\theta \right)}{d\theta }{\mathrm{\Delta }}_{\chi }(t,t_0).$$

(25)

Following the ideas introduced in Roy et al., ⁴⁴ we propose the adaptive controller:

$$u_I=-(\hat{\lambda }+\alpha )\frac{s}{\parallel s\parallel }-\beta s-\stackrel{·}{\theta }\frac{{\mathrm{dB}}^T\left(\theta \right)}{d\theta }{\mathrm{\Delta }}_{\chi },$$

(26)

where $s$ is computed according to (23), the constants $\alpha >0$ and $\beta >0$ , and the estimated $\hat{\lambda }$ evolves as:

$$\stackrel{·}{\hat{\lambda }}=\parallel s\parallel ,$$

(27)

with $\hat{\lambda }\left(0\right)>0$ . First of all, because we must note that $\hat{\lambda }\left(0\right)>0$ , then $\hat{\lambda }\left(t\right)>0$ , for all $t\ge 0$ . Substituting (26) in (25), we obtain:

$$\stackrel{·}{s}=-(\hat{\lambda }+\alpha )\frac{s}{\parallel s\parallel }-\beta s+w$$

(28)

Now, to perform the convergence analysis, we use the following Lyapunov candidate function:

$$V=\frac{1}{2}{s}^{T}s+\frac{1}{2}{(\hat{\lambda }-\overline{w})}^2.$$

(29)

Recall that $\parallel w\left(t\right)\parallel \le \overline{w}$ . Now, calculating the time derivative of (29) around the trajectories of (28) and (27), we obtain the following:

$$\stackrel{.}{\stackrel{·}{V}=s^T(-(\hat{\lambda }+\alpha )\frac{s}{\parallel s\parallel }-\beta s+w)+(\hat{\lambda }-\overline{w})\parallel s\parallel }.$$

(30)

Using $\parallel w\left(t\right)\parallel \le \overline{w}$ , is easy to see that (30) can be upper bounded by:

$$\stackrel{·}{V}\le -\beta \parallel s{\parallel }^2-\alpha \parallel s\parallel .$$

(31)

Please note that $\stackrel{·}{V}\le 0$ . However, $V$ is positive definite, bounded from below, and a non-increasing function relying on the arguments $s$ and $\hat{\lambda }-\overline{w}$ , implying that $\hat{\lambda },$ and $s$ belong to $L_{\infty }$ . To ensure the convergence of $s$ , we apply the Barbalat lemma. Then, integrating the equation (31) we can obtain the following inequality:

$$\beta {\int }_0^t\parallel s\left(\xi \right){\parallel }^{2}d\xi \le V\left(0\right)<\infty .$$

(32)

From the above, we conclude that $s\in L_2\cap L_{\infty }$ . Additionally, from the equation (28), we can conclude that $\stackrel{·}{s}$ is bounded. According to the lemma of Barbalat, ⁴⁵ we can ensure that $s\to 0$ as long as $t\to \mathrm{\infty }$ .

It is worth mentioning that if $s\to 0$ , then the proposed adaptive controller (26) and (27) effectively counteracts the effect of $w$ because it is inside the span of $B\left(\theta \right)$ (see Castanos and Fridman ³⁶ ). In practice, the evolution of $\hat{\lambda }$ , given in (27), can grow a lot because $s\to 0$ when the time $t$ is very large. We can avoid this undesirable phenomenon by proposing the $\lambda$ estimation evolution as:

$$\stackrel{·}{\hat{\lambda }}=\parallel s\parallel -\gamma \hat{\lambda },$$

with $0<\gamma <<1$ and $\hat{\lambda }\left(0\right)>0$ . Under these circumstances, we can ensure that $s$ can be small enough according to the size of $\gamma$ . In fact, we can achieve a stability analysis following the ideas found in Roy et al., ⁴⁴ ensuring that $s$ can be bounded by:

$$\parallel s\parallel \le \sqrt{\frac{\gamma {\overline{w}}^2}{2\beta }}\forall t\ge T_1,$$

where $T_1$ is the time when signal $s$ enters the confinement region; if we select the parameter control $\gamma$ sufficiently small, then the confinement ratio becomes very small, implying that the perturbation $w$ can be almost eliminated.

We summarize the above discussion in the following proposition, which is the main contribution of this study.

• Proposition 1. Let us consider the PVTOL model (4), under assumptions A1 and A2 and the restriction (11), in closed-loop with the total controller (6), composed by the nominal feedback control law $u_N$ defined in (15) and (19), where the set of constants $\{k_p,k_d,k_0,k_1\}$ are selected such that their corresponding characteristic polynomials:

$$p(s)=s^2+k_{d}s+k_p\mathrm{and}p(s)=s^2+k_{1}s+k_0,$$

are Hurwitz, together with the adaptive integral controller $u_I$ defined in (26) and (27). Then, the tracking errors globally asymptotically converge to zero.

In general, the control constants are selected such that their corresponding characteristic polynomials have the following structure:

$$p\left(s\right)=s^2+2w_e\varphi s+w_e^2,$$

where $w_e\ge 1$ and $\varphi =\{1,1.2\}$ .

It is worth mentioning that, maybe, the fundamental integral sliding mode control advantages lie in its ability to achieve almost-complete disturbance rejection from the beginning, in contrast to conventional sliding mode schemes where robustness is only ensured after the reaching phase is achieved (the interested reader can find a detailed treatment of this topic in Javaid et al. ⁴⁶ ).

Numerical simulations

In this section, we design three numerical experiments: two are control tasks consisting of tracking an ellipse-shaped trajectory, and the third is a point-to-point translation tracking a Bezier polynomial. This last experiment allows comparing the proposed control approach with two well-established solutions found in Castillo et al. ²¹ and Cardenas and Aguilar, ³¹ which are based on slide mode and Sontag’s formula, respectively. To accomplish these experiments, we used the following setup. The control gains are set as:

$k_p$	$k_d$	$k_0$	$k_1$	$a$	$b$	$\alpha$	$\beta$
1.44	1	1.44	1	1	0.7	0.5	1

For simplicity and without loss of generality, we conducted all experiments using normalized coordinates. Then, we chose the following normalized functions to simulate the crosswind:

$$\begin{array}{c}{w}_1=\frac{f_{w_1}+f_{w_2}}{\mathrm{mg}}=0.3+0.3{\sin }^2\left(\frac{t}{3}\right);\hfill \\ \multicolumn{1}{c}{w_2=\frac{L(f_{w_1}-f_{w_2})}{\mathrm{Jg}}=0.25+0.4\cos \left(\frac{t}{3}\right)\sin \left(\frac{t}{3}\right).}\end{array}$$

and we used as the saturation function $\sigma [x]=\tanh \left(x\right)$ , and we estimated the corresponding time derivatives of signals in (20) using the filter (21). The numerical simulations were carried out using MATLAB R2014a on a 64-bit personal computer running Windows 8, with an Intel(R) Core i5-4440 CPU @ 3.1 GHz and 8.00 GB of RAM. Finally, we will introduce the performance indices related to the tracking error:

$$\begin{array}{c}\multicolumn{1}{c}{\begin{array}{cc}\mathrm{IAE}={\int }_0^t\left|e\right|d\tau \hfill & \mathrm{ITAE}={\int }_0^t\tau \left|e\right|d\tau \hfill \\ \multicolumn{1}{c}{\mathrm{MSE}=\frac{1}{t}{\int }_0^t\parallel e\parallel d\tau }& \mathrm{ISE}={\int }_0^t\parallel e{\parallel }^{2}d\tau \hfill \end{array}}\\ \multicolumn{1}{c}{\left|e\right|=\left|e_{T_x}\right|+\left|e_{T_y}\right|+\left|{\stackrel{·}{e}}_{T_x}\right|+\left|{\stackrel{·}{e}}_{T_y}\right|+\left|\theta \right|+\left|\stackrel{·}{\theta }\right|}\\ \multicolumn{1}{c}{\parallel e\parallel =\sqrt{e_{T_x}^2+e_{T_y}^2+{\stackrel{·}{e}}_{T_x}^2+{\stackrel{·}{e}}_{T_y}^2+{\theta }^2+{\stackrel{·}{\theta }}^2}}\end{array}$$

• Experiment 1: The reference trajectory that the PVTOL must track is defined as:

$$\begin{array}{cc}\multicolumn{1}{c}{x_r=\sin (t/4);}& y_r=1+0.5\cos (t/4).\hfill \end{array}$$

We set the PVTOL initial conditions as:

$$\begin{array}{ccc}{x}_0=0.2;\hfill & {\theta }_0=1.6\mathrm{rad};\hfill & y_0=0.2;\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_0=0;}& {\stackrel{·}{\theta }}_0=0;\hfill & {\stackrel{·}{y}}_0=0,\hfill \end{array}$$

with $\widehat{{\lambda }_0}=1$ . We show the outcome of this experiment in Figures 2 and 3. Figure 2 shows the ellipse phase space and the positions of variables $x$ , $y$ , and $\theta$ . There, we can see that, after 15 s, the controller can nullify the crosswind effect, and the horizontal and vertical positions effectively track the reference trajectory. Please note that we observe overshoots during the first 15 s because we set the initial PVTOL angular position far from the origin. Also, in the third graph of this figure, we see that the angular position exhibits an oscillatory behavior very close to the origin.

Figure 2.

Proposed control strategy performance when tracking an ellipse-shaped trajectory, with the initial condition set at $x\left(0\right)=0.2$ , $y\left(0\right)=0.2$ , and $\theta \left(0\right)=2$  rads.

Figure 3.

Controllers $u_{T_1}$ and $u_{T_2}$ behavior in the presence of crosswind when tracking an ellipse-shaped trajectory.

In other words, our control approach ensures the whole-state stability, and the crosswind effect on the output variables $x$ and $y$ is almost imperceptible. In Figure 3, we show the controllers $u_{T_1}$ and $u_{T_2}$ , which correspond to the exerted effort needed to accomplish the task. As with the position behavior, we can also observe overshoots in both controllers, which vanished after 10 s. Finally, in Figure 4, we show the Integral of Absolute Error (IAE), Integral Time Absolute Error (ITAE), Mean Squared Error (MSE), and Integral Squared Error (ISE) indices. In this figure, we can see that MSE converges to $0.39$ after $90$ seconds, and ISE converges to $10.79$ after $12$ seconds. The IAE and ITAE indices exhibit increasing behavior due to their accumulative nature. Additionally, we observe that the IAE and ITAE indices increase due to the ellipse-tracking maneuver. We can explain it because the variables $\theta$ and $\stackrel{·}{\theta }$ exhibit slight oscillatory behavior. In other words, variables $r_x$ and $r_x$ are different from zero. On the other hand, the first and second time derivatives of $r$ produce a slight numerical error when computing the controller $u_{N_2}$ .

• Experiment 2: Using the same setup as in the previous experiment, except for the angular position, which we set to ${\theta }_0=2$ rad, this experiment had the PVTOL perform point-to-point translation tracking a Bezier’s polynomial. In Figure 5, we present the experiment obtained results related to the control task, where we can see that the PVTOL effectively accomplishes the point-to-point translation from the initial position $q_i=(x_i=0,y_i=0)$ to the final $q_f=(x_f=1,y_f=2)$ during the time interval $[t_i=0,t_f=20]$ . We specified the reference trajectories $x_r$ and $y_r$ using the following well-known Bezier’s polynomial ⁴⁷ :

$$\begin{array}{c}\delta (t_i,t_i,t_f)=0;\hfill \\ \multicolumn{1}{c}{\delta (t_f,t_i,t_f)=1,}\end{array}$$

where a finite time of time derivatives of this interpolating polynomial was set to be equal to zero at the initial and final maneuvers times, $t_i$ and $t_f$ . We use the following functions to describe the reference trajectory:

$\begin{array}{c}{x}_r=x_i+[x_f-x_i]\delta (t,t_i,t_f);\hfill \\ \multicolumn{1}{c}{y_r=y_i+[y_f-y_i]\delta (t,t_i,t_f),}\end{array}$

that suitably interpolate between the initial and final values of displacement. In the mentioned figures, we can see that the aircraft effectively translates from the position $q_i$ to position $q_f$ , even when $q_i$ is far away from the initial position, particularly ${\theta }_0$ . In Figure 6, we show the behaviors of control actions of $u_{T_1}$ and $u_{T_2}$ . There, we can see that the controller $u_{T_2}$ exhibits overshoots due to the angular position initial distance. This behavior vanished after 12 s, a moment when both controllers were almost at steady state. Complementarily, both controllers counteract the destabilizing effects of the crosswind-induced matching perturbation.

Figure 4.

The corresponding IAE, ITEA, MSE, and ISE performance indexes for the ellipse-shaped trajectory.

Figure 5.

Proposed control strategy performance when tracking a Bézier’s spline-shaped trajectory, with the initial condition set at $x\left(0\right)=0.2$ , $y\left(0\right)=0.2$ , and $\theta \left(0\right)=2$  rads.

Figure 6.

Controllers $u_{T_1}$ and $u_{T_2}$ behavior in the presence of crosswind when tracking a Bézier’s spline-shaped trajectory.

Finally, in Figure 7, once again we show the IAE, ITAE, MSE, and ISE indices for the point-to-point translation task. This time, we can see that the four indices converge, mainly because the angular position and velocity converge to zero. Besides, the variable $r$ and its corresponding time derivative converge to zero for this particular trajectory.

• Experiment 3. Once again, using the same setup as in the previous numerical experiments, we ran a comparison between our control proposal and the two well-established robust control approaches, one based on Robust Control Lyapunov Functions (RCLF) and Sontag’s universal stabilizing feedback, and the other based on output feedback sliding mode, respectively introduced in Munoz et al. ²⁰ and Cardenas and Aguilar ³¹ The control task consists of tracking the following piece-wise trajectory:

$$\begin{array}{cc}\multicolumn{1}{c}{x_r=\{\begin{array}{c}1,t\le 25s,\hfill \\ \multicolumn{1}{c}{-1,t>25s.}\end{array};}& y_r=\{\begin{array}{c}1.5,t\le 25s,\hfill \\ \multicolumn{1}{c}{2.5,t>25s.}\end{array},\hfill \end{array}$$

and we set the origin as the initial condition, but we set the angular position to $\theta \left(0\right)=0.5\mathrm{rad}$ . From now on, we use the letters O, M, and C to, respectively, refer to our controller (O), the controller proposed by Munoz et al. ²⁰ (M), and the controller introduced by Cardenas and Aguilar ³¹ (C). Finally, we chose the same control parameters for O as in the previous experiments, while M and C used the same parameters as in their corresponding publications.

Figure 7.

The corresponding IAE, ITEA, MSE, and ISE performance indexes for PVTOL point-to-point translation tracking a Bezier’s polynomial task.

In Figure 8, we show the simulation outcomes, where we can see O, M, and C converge to coordinate $x$ at 17, 20, and 4 s, respectively; at 20 s, the trajectory reference changes from 1 to −1 m, and, once again, O, M, and C converge to coordinate $x$ at 40, 30, and 40, respectively. On the other hand, M never tracks the reference trajectory coordinate $y$ because needs more time, C tracks it close, and O converges at 15 s. We observe a similar behavior after the trajectory change from 0 to 2.3 m. Finally, the three control strategies converge to coordinate $y$ , where C converges in less time. However, overshoots of a large magnitude can be observed. In the same figure, we see that the velocity at coordinate $x$ is reached by C after 3 s, and O and M do it at 16 and 18 s, respectively; the three strategies reached the velocity at coordinate $y$ close to 15 s, being C the one who exhibits the largest overshoots. Finally, the angular velocity is reached by O and M almost immediately, while C needed close to 3 s, exhibiting, once again, overshoots at the beginning and when the reference changes abruptly. We omitted the controllers’ plots because they exhibit chattering, making it difficult to interpret them and obtain conclusions.

Figure 8.

Comparison of the closed-loop response of O, M, and C when performing a piece-wise reference trajectory.

We note that this comparison provides only a rough idea of how effective our control solution is compared to the other two approaches. Due to the nature of our proposal, the comparison is unfair because we included an integral action that completely counteracts the external perturbation from the beginning of the experiment. On the other hand, the previously introduced solutions in the literature neglect the crosswind effect when the aircraft’s position is located at the singular position $\pi /2$ or when $u_{T_1}=0$ .

Finally, in Table 1, we present the corresponding statistical performance indices to quantify how our proposal compares with the other two methods; we observe that ours performs better.

• Remark 2. Given the theoretical nature of our control approach, we note that it is a common practice to design and propose control approaches under idealized mathematical models and assumptions; these approaches can then be numerically validated for performance, as we did in this study. In fact, we considered unmodeled matching perturbations. In real-world applications, implementers must counteract the undesirable effects caused by position sensor noise. It can be accomplished using a low-pass filter, allowing the measurement of unavailable velocities, which can then be evaluated using, for instance, an exact Levant differentiator, or a suitable version of the Kalman filter.^48,49Also, limitations due to the actuators’ saturation have to be addressed. A straightforward solution to this issue can be convenient trajectory planning. On the other hand, the unknown and unmodeled non-matching dynamics can negatively affect the actuators’ performance. To overcome this inconvenience, we can introduce an $H_{\infty }$ technique. A detailed analysis of these scenarios is beyond the scope of our solution, as already mentioned, due to its theoretical nature. Although our approach may ensure nominal stability and asymptotic convergence of the tracking error to zero, in the presence of unmodeled matching perturbations, its performance is considerably reduced, leading to increased tracking error, oscillations around the reference trajectory, or, in the worst case, loss of stability. For instance, if the integral control action fails, the performance degrades; as long as the external perturbation is sufficiently small, the system tends, as mentioned, to oscillate around the planned trajectory.

• Remark 3. The overshoots and oscillatory transients observed in the control inputs are mainly due to two factors: large initial tracking errors and the use of the tracking differentiator (20) to compute the first and second time derivatives of the function $r$ (see 16). As shown in the simulations, these behaviors are transient and vanish as the system approaches the vicinity of zero. In practice, this can degrade the control strategy’s performance due to the inherent saturation of the PVTOL motors. Concerning chattering, although the integral sliding-mode controller contains a discontinuous term, we can mitigate it in practice by using smooth saturation functions. Furthermore, in our implementations, chattering is reduced by replacing the discontinuous sign function with $\mathrm{sgn}[x]=x/(\left|x\right|+{10}^{-3})$ , a well-established technique in sliding-mode control. This claim is consistent with the simulation results, which show no high-frequency oscillations in the control inputs. In general, using integral slide mode control produces a lower frequency compared with the first-order slide mode control method.

Table 1.

Performance index for $t=50$  s.

Method	IAE	ITAE	MSE	ISE
O	40.55	630.90	0.41	28.55
M	44.01	682.72	0.57	100.12
C	60.84	891.53	0.73	64.99

Conclusions

We designed a global, feedback-based control solution for trajectory tracking of a crosswind-perturbed crewless PVTOL aircraft. This solution applies a backstepping version in conjunction with a saturation function. The controller works as follows: to prevent the PVTOL from reaching a singular angular position, the accelerations in the horizontal and vertical aircraft positions are set to convenient values. The latter enables us to obtain a continuous controller for trajectory tracking and an integral adaptive sliding-mode controller for canceling the unknown bounded crosswind effect. We designed a continuous controller using an implicit version of the Immersion and Invariance method, without using feedback linearization or differential flatness. We numerically tested the performance of our control solution across three simulations, comparing it with two other well-established related control solutions. The numerical simulation outcomes allowed us to claim that, in general, our proposal performance can equate with those of previously introduced approaches, even when, in our case, the time derivatives of signals (20) were computed using a tracking differentiator to avoid solving a large number of algebraic expressions.

Stability proof

Let us write the equation of (13) in space form as

$$\begin{array}{c}{\stackrel{·}{e}}_1=e_2\hfill \\ {\stackrel{·}{e}}_2=-{\sigma }_m[k_p{e}_1+e_2]\hfill \end{array}$$

(33)

where for simplicity, we fix $k_d=1$ , and $m>0$ . Let us use the following Lyapunov function

$$V={\int }_0^{e_1}{\sigma }_m\left[k_{p}s\right]\mathrm{ds}+\frac{1}{2}{e}_2^2,$$

whose time derivative around the trajectories of (33) leads to:

$$\stackrel{·}{V}=-e_2({\sigma }_m[k_p{e}_1+e_2]-{\sigma }_m\left[k_p{e}_1\right])$$

(34)

Evidently, using ${\mathbf{P}}_3$ , we have that $\stackrel{·}{V}(e_1,e_2)\le 0$ . As $V$ is a strictly positive and non-increasing function, we can conclude that $e_1$ and $e_2$ are bounded. Now, from LaSalle’s invariance principle, ⁴⁵ all solutions of the system (33) converge to the largest invariant set in:

$$\mathrm{\Omega }=\{(e_1,e_2)\in {\in \mathfrak{R}}^2:\stackrel{·}{V}=0\}.$$

(35)

According to (34), we have that $e_2({\sigma }_m[k_p{e}_1+e_2]-{\sigma }_m\left[k_p{e}_1\right])=0$ if and only if $e_2=0$ , meaning that inside the invariant set (35) holds that ${\stackrel{·}{e}}_2=0$ and $e_1=\overline{e}$ , with $\overline{e}$ constant. From the second equation of (33), we have that $0=-{\sigma }_m\left[k_p\overline{e}\right]$ . That is, $\overline{e}=0$ . Hence, according to the invariant theorem of LaSalle, we conclude that $e_1\to 0$ and $e_2\to 0$ as long as $t\to \mathrm{\infty }$ . In other words, the system is globally asymptotically stable. Using simple linearization of the system (33), we can also conclude that the system is locally exponentially stable.

Expression of the variables, $r^{\prime }$ , $r_x^{\prime }$ , $r_x^{″}$ , $r_y^{\prime }$ and $r_y^{″}$

For the following developments, we use the following notation: according to the definition (11), we introduce the following variables

$$\begin{array}{cc}\multicolumn{1}{c}{r_x=-{\sigma }_a[X]+{\stackrel{··}{x}}_r}& r_y=-{\sigma }_b[Y]+{\stackrel{··}{y}}_r\hfill \end{array}$$

where $X=k_p(x-x_r)+k_d(\stackrel{·}{x}-{\stackrel{·}{x}}_r)$ and $Y=k_p(y-y_r)+k_d(\stackrel{·}{y}-{\stackrel{·}{y}}_r)$ . From this definition, it is easy to see that:

$$\begin{array}{c}\multicolumn{1}{c}{r_x^{\prime }=-{\sigma }_a^{\prime }[X](k_p(\stackrel{·}{x}-{\stackrel{·}{x}}_r)+k_d(\stackrel{··}{x}-{\stackrel{··}{x}}_r))+x_r^{\left(3\right)}}\\ \multicolumn{1}{c}{r_x^{″}=-{\sigma }_a^{″}[X]{(k_p(\stackrel{·}{x}-{\stackrel{·}{x}}_r)+k_d(\stackrel{··}{x}-{\stackrel{··}{x}}_r))}^2-}\\ \multicolumn{1}{c}{{\sigma }_a^{\prime }[X](k_p(\stackrel{··}{x}-{\stackrel{··}{x}}_r)+k_d(x^{\left(3\right)}-x_r^{\left(3\right)}))+x_r^{\left(4\right)}}\end{array},$$

and:

$$\begin{array}{c}\multicolumn{1}{c}{r_y^{\prime }=-{\sigma }_b^{\prime }[Y](k_p(\stackrel{·}{y}-{\stackrel{·}{y}}_r)+k_d(\stackrel{··}{y}-{\stackrel{··}{y}}_r))+y_r^{\left(3\right)}}\\ \multicolumn{1}{c}{r_y^{″}=-{\sigma }_b^{″}[X]{(k_p(\stackrel{·}{y}-{\stackrel{·}{y}}_r)+k_d(\stackrel{··}{y}-{\stackrel{··}{y}}_r))}^2-}\\ \multicolumn{1}{c}{{\sigma }_b^{\prime }[Y](k_p(\stackrel{··}{y}-{\stackrel{··}{y}}_r)+k_d(y^{\left(3\right)}-y_r^{\left(3\right)}))+y_r^{\left(4\right)}}\end{array},$$

where:

$$\begin{array}{cc}\multicolumn{1}{c}{x^{\left(2\right)}=-S_{\theta }{u}_{N_1};}& x^{\left(3\right)}=-\stackrel{·}{\theta }{C}_{\theta }{u}_{N_1}-S_{\theta }{\stackrel{·}{u}}_{N_1;}\hfill \\ \multicolumn{1}{c}{y^{\left(2\right)}=-1+C_{\theta }{u}_{N_1};}& y^{\left(3\right)}=-\stackrel{·}{\theta }{S}_{\theta }{u}_{N_1}+C_{\theta }{\stackrel{·}{u}}_{N_1;}\hfill \end{array}$$

where ${\stackrel{·}{u}}_{N_1}$ is computed as:

$${\stackrel{·}{u}}_{N_1}=-r_x^{\prime }\sin \theta -r_x\stackrel{·}{\theta }\theta -(r_y+1)\stackrel{·}{\theta }\theta +r_y^{\prime }\theta .$$