Stabilization of an Unconventional Large Airship When Hovering

Featured Application

Large capacity airships (LCA) began the last straight line before the start of their commercial flights during this decade. Many countries are interested in the tremendous potential of these devices and the growing number of applications that could be assigned to them. These include, for example, the transport of logs from areas that are difficult to access, the transport of wind turbine blades to areas that are often remote and far from infrastructure, and the loading and unloading of container ships on the high seas in regions where there are no adequate structures to accommodate these giants of the seas. Additionally, the list goes on.

Abstract

In this paper, we present the stabilization of an unconventional unmanned airship above a loading and unloading area. The study concerns a quad-rotor flying wing airship. This airship is devoted to freight transport. However, during the loading and unloading phases, the airship is very sensitive to squalls. In this context, we present in this paper the dynamic model of the airship, and we propose a strategy for controlling it under the effects of a gust of wind. A feedforward/feedback control law is proposed to stabilize the airship when hovering. As part of the control allocation, the non-linear equations between the control vectors and the response of the airship actuators are highlighted and solved analytically through energy optimization constraints. A comparison with classical numerical algorithms was performed and demonstrated the power and interest of our analytic algorithm.

1. Introduction

Interest in airships has boomed in this century after a long hibernation due to various factors. This interest has focused on missions that differ from those entrusted to airships in the first half of the 20th century. They relate, in particular, to the use of stratospheric airships [1,2] as surveillance or communication satellites or even the use of large airships for freight transport. Much research has been carried out on the latter subject. Let us quote, for example, the AIRLANDER manufactured by the company Hybrid Air Vehicles, the most successful of these large airships, which made more or less conclusive demonstration flights. The latest version of this airship showed weaknesses during low-speed flight in the approach phase, ultimately leading to a crash in 2016. Other examples of airships from the turn of the century can be seen in [3,4].

The most prominent objectives for these airships are the transport of heavy loads (mobile homes, wind turbine blades, field hospitals, etc.) to places lacking infrastructure, or the unloading of container ships in the absence of a port, and many other applications where other means of air transport would be inefficient or too expensive.

However, in order to achieve this potential, several problems must be resolved, in particular, the stabilization of these machines in the presence of gusts of wind, mainly when hovering during the critical phases of loading and unloading.

Several research laboratories are also seized with this interesting subject. Among the airships developed by laboratories, we note the MC500 (Figure 1) with a capacity of 500 m³ developed by the interuniversity group DIRISOFT in France and which will serve as a support for our developments.

The MC500 is original by the unconventional shape of its hull and the plethora of novelties that characterize it. The advancement of digital aerodynamics, control techniques, and on-board electronics have made it possible to design airships of unconventional shapes, see, for example, [5,6]. The use of an unconventional shape was guided by the search for an optimal way to capture sunlight using photovoltaic films. This requires having large “flat surfaces”, hence the choice of a flying wing.

For the MC500, the choice of the quadcopter configuration with steerable rotors is dictated by a concern for the redundancy of the actuators in order to better control the machine with unusual dimensions and so that they can stabilize quickly and resist a gust of wind detected. It should be emphasized that this tends to complicate the task of developing stabilization control laws that would be compatible with the capacities of the actuators.

Among the new features of the airship, the designers planned to equip it with a LIDAR. This device represents a step forward in carrying out adequate checks to minimize the drift of these bulky objects under the effects of gusts of wind. LIDAR is an optical measurement technique based on the analysis of the properties of a light beam returned to its emitter [7]. LIDAR measures the wind speed upstream of the flying object in order to provide an alert in the event of dangerous turbulence being detected. This technology offers the possibility of using new control concepts for airships, such as feedforward control algorithms.

In this sense, in this article, we present a feedforward control; the objective being to minimize the drift of the airship during these atmospheric turbulences.

Several control approaches have been proposed for these aerial objects, such as linear control [8,9], feedback linearization [9,10], adaptive control [11], LQR based controllers [12], or backstepping control [13]. In reference to [11], two control approaches to analyzing the stability and robustness of the AURORA self-propelled airship are presented. The first is based on linearized models of the airship, and the second is based on nonlinear control techniques such as dynamic inversion and backstepping. The established nonlinear control laws were tested against a simulated wind disturbance. The backstepping approach is extended in [14] to include input saturation and wind effect compensation. A wind estimator was presented. In [15], active disturbance rejection control (ADRC) based on a nonlinear extended state observer (ESO) has been applied for the trajectory tracking of a stratospheric airship in the horizontal plane. The simplified model of the airship in this plane is transformed into a two single input and single output subsystem using the input-output linearization method. The two equivalent linear systems are separately processed in order to use the observer (ESO) in real time to estimate the model parameters and to attenuate the wind force that affects the airship. In the same field of trajectory tracking, Yang [16] chose to use sliding mode control. In his study, he suggested a sliding manifold composed of predefined nonlinear functions.

The previous control algorithms are based on the robustness technique and the observers to compensate for the uncertainties of the model and the effect of an external perturbation. However, the presence of a strong wind causes significant effects on the airship, and stabilization is difficult to obtain when the wind speed exceeds the limit required by the robustness of the control. During a measured disturbance, the use of advance controls allows us to compensate for the effects of an exogenous signal before it influences the system. The minimization of the effects of a perturbation on the output to be controlled is treated as a problem of decoupling between the output of the system and the exogenous input considered.

This problem has been solved using differential geometric analysis tools [17,18]. An alternative formulation based on the relative degree has been introduced in [19,20]. The authors have shown that the coupling problem is solved if the relative degree of the perturbation is greater than or equal to that of the output.

We applied this last methodology to our airship and we demonstrated the power of this technique and its interest in the control and stabilization of these large flying machines.

A study on a robust method is in progress. It is based on the model-free control theory developed by Fliess [21]. This technique is freed from the dynamic model and its uncertainties and seems promising for the stabilization of our airship. The results of this method will be published shortly.

A problem closely related to stabilization but little traced in the literature concerns the allocation of control. In order to prove that the established control can be implemented on the machine, it would be necessary to verify the ability of the actuators to apply them. The airship actuators, represented by the angles of orientation and the forces produced by the rotors, actually have saturation limits. For the MC500 airship, the subject of this study, the actuators are twelve, but the “virtual” controls are six. The connection between the “virtual” controls and the actuators is thus described by a strongly nonlinear rectangular system, which makes solving this complex system a challenge. Usually, the control allocation is processed numerically. In this context, we cite the pseudo-inverse method used to solve the problem of unconstrained control allocation [22], while the linear quadratic optimization method with constraints is used in [23,24,25]. These methods use numerical techniques to find the minimum of a quadratic cost function. However, the question of the applicability of these algorithms in real time in the face of an external disturbance arises. In fact, when a gust of wind occurs, the robustness margin of the virtual controls can be exceeded, so the objective function to be minimized cannot admit a minimum, and there will be a risk of instability of the airship.

This is why, in this study, we propose an original algorithm based on algebraic equations determined analytically according to energy arguments. This new algorithm will allow the airship to be piloted even in bad weather conditions. A comparative study with a numerical algorithm is presented.

We organized this paper as follows: in Section 2 a dynamic model of the airship MC500 is proposed, in Section 3 we precisely present the stabilization strategy developed and a control algorithm based on a feedforward/feedback technique. In Section 4 we present the analytical algorithm developed for the control allocation that we compare with other numerical techniques. Finally, in Section 5, numerical simulations are presented.

2. Modeling

2.1. Kinematic Description

For the kinematic description we used a Galilean reference frame fixed to the ground $R_0=(O,X_{0,}{Y}_0,Z_0)$ and a mobile reference frame $R_m=(O,X_{m,}{Y}_m,Z_m)$ linked to the airship at its center of mass G. As we can see in Figure 2, the airship MC 500 is a quadrotor. Each rotor can swivel around two axes: β_i around the Y_m axis (Figure 2), and γ_i around an axis Z_Ri, orthogonal to Y_m and which coincides in the beginning with the Z_m axis. The value of the thrust of the rotor i is denoted by $\Vert {\mathrm{F}}_{\mathrm{i}}\Vert$, and its position in R_m is P_i, such as: P₁ = (a, b₁, c)^T; P₂ = (a, −b₁, c)^T; P₃ = (−a, b₃, c)^T; P₄ = (−a, −b₃, c)^T.

The orientation of the airship is described through three rotations defined by the classical Euler angles (yaw $\psi$, pitch $\theta$, and roll $\varphi$). The rotation matrix between the mobile frame R_m and the fixed reference frame R₀ is then given by:

$$J_1^T=\left(\begin{array}{ccc}c\theta .c\psi & c\theta .s\psi & -s\theta \\ s\varphi .s\theta .c\psi -c\varphi .s\psi & s\varphi .s\theta .s\psi +c\varphi .c\psi & \mathrm{s}\varphi .c\theta \\ c\varphi .s\theta .c\psi +s\varphi .s\psi & c\varphi .s\theta .s\psi -s\varphi .c\psi & \mathrm{c}\varphi .c\theta \end{array}\right)$$

where: $c\hspace{0.17em}\theta =cos\theta \hspace{0.33em}\hspace{1em}and\hspace{0.33em}\hspace{0.33em}s\hspace{0.17em}\varphi =sin\varphi \hspace{0.33em}$

Using the rotation matrix $J_1$, the expression of the linear speed in the fixed reference frame ${\dot{\eta }}_{\hspace{0.17em}1}={[{\dot{x}}_{\hspace{0.17em}0},\hspace{0.33em}{\dot{y}}_0,\hspace{0.33em}{\dot{z}}_0]}^{\hspace{0.17em}T}$ as a function of the speed expressed in the mobile frame ${\nu }_{\hspace{0.17em}1}={[\hspace{0.33em}u,v,w\hspace{0.33em}]}^{\hspace{0.17em}T}$ is given by:

$${\dot{\eta }}_1=J_1.{\nu }_1$$

(1)

where ${\eta }_1={[x_0, y_0, z_0]}^T$ is the position vector of the origin of the mobile frame R_m expressed in the fixed frame $R_{\hspace{0.17em}0}$.

On the other hand, the angular speed of the airship ${\nu }_2^T=\left[p, q, r\right]$ expressed in the mobile reference frame is the combination of the angular speeds around the three axes of yaw, pitch, and roll. It can be written related to ${\dot{\eta }}_2^T=\left[\dot{\varphi }, \dot{\theta }, \dot{\psi }\right]$ as:

$${\nu }_2=\left(\begin{array}{ccc}1& 0& -\mathrm{s}\theta \\ 0& \mathrm{c}\varphi & \mathrm{s}\varphi \hspace{0.17em}\hspace{0.17em}\mathrm{c}\theta \\ 0& -\mathrm{s}\varphi & \mathrm{c}\varphi \hspace{0.33em}\mathrm{c}\theta \end{array}\right).\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)$$

(2)

Conversely:

$${\dot{\eta }}_2=J_2.{\nu }_2$$

(3)

where the transformation matrix J₂ is represented by: $J_2=\left(\begin{array}{ccc}1& \mathrm{s}\varphi \tan \theta & \mathrm{c}\varphi \tan \theta \\ 0& \mathrm{c}\varphi & -\mathrm{s}\varphi \\ 0& \frac{\mathrm{s}\varphi }{\mathrm{c}\theta }& \frac{\mathrm{c}\varphi }{\mathrm{c}\theta }\end{array}\right).$

It is important to mention that the parametrization by the angles of Euler presents a singularity for $\theta =\frac{\pi }{2}+k\pi$; however, this configuration is not reachable for the airship.

The kinematics of the airship can then be written as:

$$\left(\begin{array}{c}{\dot{\eta }}_1\\ {\dot{\eta }}_2\end{array}\right)=\left(\begin{array}{cc}{J}_1& 0\\ 0& J_2\end{array}\right)\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\end{array}\right)$$

(4)

2.2. Dynamics

In a compact form, the whole dynamic system of the airship can be expressed as follows [26]:

$$M.\dot{\nu }=\tau +Q_G$$

(5)

With $\nu =\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\end{array}\right)$ as the velocity vector and $M=\left(\begin{array}{cc}{M}_{TT}& 0\\ 0& M_{RR}\end{array}\right)$ as the mass matrix. M_TT and M_RR are, respectively, the translation and the rotation part of the mass matrix of the airship. We denote by M_ij the components of the 6 × 6 matrix M.

Let us note that this mass matrix comprises the terms of inertia of a solid body, but also the terms of added masses. The acceleration of a moving body in a fluid creates a resistance force exerted by the latter, proportional to the acceleration of the solid. This force can therefore be likened to an increase in the mass of the solid, hence the notion of added mass. The values of these added masses are important and not negligible in the case of submarines or airships [27]. The computation of these added masses can be done thanks to experiments, or analytically as was done, for example, for the MC500 [28].

In the right side, $\tau =\left(\begin{array}{c}{\tau }_1\\ {\tau }_2\end{array}\right)$, τ₁, and τ₂ are, respectively, the external forces and torques, including the rotors effects, the weight (m.g), the buoyancy B, the wind Force F_v, and the aerodynamic lift (F_L) and drag (F_D). While $Q_G=\left(\begin{array}{c}-{\nu }_2\wedge \left(M_{TT}{\nu }_1\right)\\ -{\nu }_2\wedge \left(M_{RR}{\nu }_2\right)-{\nu }_1\wedge \left(M_{TT}{\nu }_1\right)\end{array}\right)$ are the gyroscopic forces and torques, we will denote by Q_Gi its components, $\wedge$ being the vector product.

When combining with the kinematic relations (4), the global dynamic model in its developed form becomes:

$$\left\{\begin{array}{c}{\dot{x}}_0=c\psi .c\theta .u+(-s\psi .c\varphi +c\psi s\varphi s\theta )v+(s\psi s\varphi +c\psi c\varphi s\theta )w\hfill \\ {\dot{y}}_0=s\psi c\theta .u+(c\psi c\varphi +s\psi s\varphi s\theta )v+(-c\psi s\varphi +s\psi c\varphi s\theta )w\hfill \\ {\dot{z}}_0=-s\theta .u+s\varphi .c\theta .v+c\varphi c\theta w\hfill \\ \dot{\varphi }=p+s\varphi .tan\theta .q+c\varphi .tan\theta .r\hfill \\ \dot{\theta }=c\varphi .q-s\varphi .r\hfill \\ \dot{\psi }=\frac{s\varphi }{c\theta }q+\frac{c\varphi }{c\theta }r\hfill \\ M_{11}\dot{u}={\alpha }_1+Q_{G1}+F_{v_1}\hfill \\ M_{22}\dot{v}={\alpha }_2+Q_{G2}+F_{v_2}\hfill \\ M_{33}\dot{w}={\alpha }_3+Q_{G3}+F_{v_3}\hfill \\ (M_{46}^2-M_{66}{M}_{44})\dot{p}=M_{46}{\alpha }_6+M_{46}{Q}_{G6}-M_{66}{\alpha }_4-M_{66}{Q}_{G4}\hfill \\ M_{55}\dot{q}={\alpha }_5+Q_{G5}\hfill \\ (M_{46}^2-M_{66}{M}_{44})\dot{r}=M_{46}{\alpha }_4-M_{44}{\alpha }_6+M_{46}{Q}_{G4}-M_{44}{Q}_{G6}\hfill \end{array}\right.$$

(6)

where F_vi are the components of the wind applied forces and α_i the components of the control vector α. The different characteristics of the actuators will be developed in detail in Section 4.

3. Stabilization Strategy

3.1. Mathematical Tools

In some applications, exogenous perturbation is known in advance. This is often the case for the airship. The MC500 is equipped with LIDAR sensors that can measure the force of a wind gust. From the control point of view, the purpose of using this information is to reduce or minimize the influence of the disturbance that affects the airship. To deal with the problem of minimizing the impact of a wind gust, we propose the use of feedforward controls. This control vector is built using techniques of differential geometry. The diagram of this method is described in Figure 3.

The feedforward/feedback method controls a set of nonlinear systems which are as follows:

$$\left\{\begin{array}{l}\dot{x}=f(x)+g(x)u+{\displaystyle \sum _{i=1}^m}{p}_i{d}_i\\ y=h(x)\end{array}\right.$$

(7)

Here, $x\in R^n$ represents the state system, y is the output system, $u\in R^m$ the input system, $g(x)={\left[g_1(x),\cdots ,g_m(x)\right]}^T$, $g_i(i=1,\cdots ,m)$ is an m-dimensional vector field, $f(x)\in R^n$ is a vector field, and $h(x)={\left[h_1(x),\cdots ,h_m(x)\right]}^T$, $h_i(i=1,\cdots ,m)$ is a sufficiently smooth scalar function, d_i the disturbance inputs, and p_i are vectors detailed in Appendix A.

Definition 1.

If we assume that the vector fields h and f are sufficiently smooth, the Lie derivative of h with respect to f can be defined as:

$$L_{f}h={\displaystyle \sum _{i=1}^n}\frac{\partial h_i}{\partial x_i}{f}_i(x)$$

(8)

Definition 2.

In the system of Equations (7), if the two ensuing conditions are satisfied for all$x\in R^n$in the vicinity of an equilibrium point x_e [17]:

$$L_{g_j} L_f^k h_i\left(x\right) = 0$$

(9)

The (m × m) matrix of decoupling D is not singular at x = x_e, with:

$$D(x)=\left(\begin{array}{ccc}{L}_{g_1}{L}_f^{r_1-1}{h}_1(x)& \dots & L_{g_m}{L}_f^{r_1-1}{h}_1(x)\\ L_{g_1}{L}_f^{r_2-1}{h}_2(x)& \dots & L_{g_m}{L}_f^{r_2-1}{h}_2(x)\\ \dots & \dots & \dots \\ L_{g_1}{L}_f^{r_m-1}{h}_m(x)& \dots & L_{g_m}{L}_f^{r_m-1}{h}_m(x)\end{array}\right)$$

(10)

Then, the system (7) is said to have a relative degree:

$$r={\displaystyle \sum _{i=1}^m\hspace{1em}{r}_i}$$

(11)

The indices I, j, and k verify: $1\le \hspace{0.33em}j,\hspace{1em}i\le m,\hspace{1em}1\le k\le r_i-1$. r_i are the relative degree of the outputs. Examples of the derivatives $L_{gj} L_f^k h_i\left(x\right)$ are presented later in Section 3.2.1.

Definition 3.

The relative order ρ_ji of the output y_i of Y with respect to the disturbance input d_i is defined [19,20] as the smallest integer for which: $L_{p_i} L_f^{{\rho }_j^i-1} h_j\hspace{0.33em}\ne \hspace{0.33em}0$

Proposal 1.

The system of Equation (7) is feedback linearizable if it exists as a function h ∈ ℝ^m sufficiently smooth such that the system has the relative degree r = n = dim(x), x is the state.

We will apply this result to demonstrate that the model of the airship MC500 subjected to a wind gust is linearizable. We can thus apply the input-output linearization method to define a control vector that anticipates the effect of a wind gust pre-detected by LIDAR.

3.2. Scheme of the Control Vector

In order to stabilize the airship when hovering in the vicinity of a loading and unloading point around a desired state $\mathsf{\Psi }={\left[x_d, y_d, z_d, {\varphi }_d, {\theta }_d, {\psi }_d\right]}^T$, and to anticipate the effect of a gust of wind, we will consider as a new output the error ${\left({\xi }_i\right)}_{1\le i\prec 6}$ defined by: ${\xi }_i=y_i-{\mathsf{\Psi }}_i$ with: $Y={\left[x_0, y_0, z_0, {\varphi }_{}, {\theta }_{}, {\psi }_{}\right]}^T$. The objective is to make the output converge towards the desired state, in other words, to make the error ${\xi }_i$ converge towards zero.

Without losing generality, we assume that $F_{{\nu }_2}=F_{{\nu }_3}=0$ (the gust of wind comes only from the X axis without creating any torque).

The system of errors associated with the system in Equation (6) could be written in the same form given by Equation (7) as:

$$\left\{\begin{array}{l}\dot{X}=f(X)+g(X)\alpha +p_1{F}_{v_1}\\ Y=h=[x_0,y_0,z_0,\phi ,\theta ,\psi {]}^t\end{array}\right.$$

(12)

With: X the system of error defined as: $X=\left[{\xi }_1, {\xi }_2, {\xi }_3, {\xi }_4, {\xi }_5, {\xi }_6, u, v, w, p, q, r\right]$, $\alpha ={\left[{\alpha }_1, {\alpha }_2, {\alpha }_3, {\alpha }_4, {\alpha }_5, {\alpha }_6\right]}^T$ the control vector, and $p_1$, f and g are vectors of ℝ¹² (see Appendix A for more details).

The input-output linearization method is summarized in two steps (see Figure 4). The first is to transform the system into a decoupled linear system using the techniques of differential geometry. The second is to construct a control vector by using the theory of linear control to ensure the stability of the new system.

3.2.1. Relative Degree of Outputs Associated with the Control Vector

Let r be the relative degree of the system (12).

To determine the relative degree ${\left(r_j\right)}_{1\le j\le 6}$ corresponding to the output y_j, this is derived last until at least one input appears.

In their work, Daoutidis [19,20] proposes the following expression for ${\xi }_j^{\left(r_j\right)}$:

$${\xi }_j^{(r_j)}=L_f^{r_j}{h}_j+{\displaystyle \sum _{i=1}^6}{L}_{g_i}{L}_f^{r_j-1}{h}_j{\alpha }_i+{\displaystyle \sum _{i=1}^6}{L}_{p_i}{L}_f^{r_{j1}-1}{h}_j{d}_1$$

(13)

Here, we have: $L_{g_i} h_1=0$

Then:

$${\dot{\xi }}_1=L_f{h}_1=f_1$$

(14)

The first derivative of ${\xi }_{\hspace{0.17em}1}$ does not include any of the controls. Therefore, another derivation of ${\xi }_{\hspace{0.17em}1}$ is necessary:

$${\ddot{\xi }}_1=L_f^2{h}_1+L_{g_1}{L}_f{h}_1{\alpha }_1+L_{g_2}{L}_f{h}_1{\alpha }_2+L_{g_3}{L}_f{h}_1{\alpha }_3+L_{p_1}{L}_f{h}_1{F}_{v_1}$$

(15)

With:

$$L_f^2{h}_1=\frac{\partial f_1}{\partial \varphi }{f}_4+\frac{\partial f_1}{\partial \theta }{f}_5+\frac{\partial f_1}{\partial \psi }{f}_6+\frac{\partial f_1}{\partial u}{f}_7+\frac{\partial f_1}{\partial v}{f}_8+\frac{\partial f_1}{\partial w}{f}_9$$

$$L_{g_1}{L}_f{h}_1=\nabla f_1.g_1=\frac{c\psi c\theta }{M_{11}}; L_{g_2}{L}_f{h}_1=\nabla f_1.g_2=\frac{-s\varphi +c\varphi s\theta }{M_{22}}$$

$$L_{g_3}{L}_f{h}_1=\nabla f_1.g_3=\frac{s\varphi +c\varphi s\theta }{M_{33}}; L_{g_j}{L}_f{h}_1=0, j=4,5,6$$

$\nabla$ is the gradient. We notice that the second derivative of ${\xi }_1$ is written as a function of the control. Hence, the degree $r_1$ associated with the output is equal to 2.

In the same way, the other relative degrees are computed, and we thus obtain:

$$r_1=r_2=\hspace{0.33em}\cdots =r_6=2$$

The total relative degree of the system is then r = 12, we conclude that the nonlinear model of the airship in Equation (11) satisfies the conditions in proposal 1. We can then deduce that the system in Equation (11) is feedback linearizable. Therefore, a state of feedback control and a diffeomorphism that transforms the nonlinear system into an equivalent linear system exist.

3.2.2. Relative Degree Associated with Disturbances

The relative degree of the output ${\left(r_j\right)}_{1\le j\le 6}$ and that of the disturbance ρ_j1 play a fundamental role in characterizing the influence of the perturbation on the outputs. Two cases may arise [19]:

• ${\rho }_{j_1}\succ r_j$ the disturbance does not directly affect the output, so the stabilization of the airship is ensured by a feedback control law.
• ${\rho }_{j_1}\le r_j$ the disturbance affects the output. We use a feedforward/feedback control to stabilize the airship and anticipate the effect of this exogenous signal.

The relative degree ρ_j1 of the disturbance associated with the output y_i is then determined using definition 3:

• for the output, y₁ = x₀, y₂ = y₀, y₃ = z₀, we have:

  $$\left\{\begin{array}{c}{L}_{p_1}{h}_1=L_{p_1}{h}_2=L_{p_1}{h}_3=0\hfill \\ L_{p_1}{L}_f{h}_1=\frac{c\psi c\theta }{M_{11}}\hfill \\ L_{p_1}{L}_f{h}_2=\frac{-s\theta }{M_{11}}\hfill \\ L_{p_1}{L}_f{h}_2=\frac{s\psi \hspace{0.17em}c\theta }{M_{11}}\hfill \end{array}\right.$$

  (16)

  then:

  $${\mathsf{\rho }}_{11}={\mathsf{\rho }}_{21}={\mathsf{\rho }}_{31}={\mathsf{\rho }}_2=0 {\mathrm{r}}_2=2$$

  (17)

For the outputs y₄ = ϕ, y₅ = θ and y₆ = ψ, we have:

$$L_{p_1}{L}_f^{-1+{\rho }_{j_1}} h_i=0\hspace{1em}\hspace{1em}\forall i=4,5,6$$

(18)

We then deduce that:

$${\rho }_{j_1}\succ r_j\hspace{1em}\hspace{1em}\forall j=4, 5, 6$$

(19)

On the basis of the two previous findings, the outputs x, y, and z are affected both by the disturbance and by the input in the same way. Therefore, to anticipate the effect of the exogenous signal, the feedforward control is necessary. We will combine it with the feedback control to realize the control schema.

However, for the other degrees of freedom (y, z, ϕ, θ, and ψ ), we noticed in Equation (19) that the effect of the inputs is preponderant compared to the effect of the disturbance of the outputs. This indicates that we do not need a feedforward control to minimize the disturbance effect.

3.2.3. Coordinate Change and New Linear System

Following the recommendations of Isidori [17], we can define a diffeomorphism Φ. This diffeomorphism allows us to transform the nonlinear system into another nonlinear system, as we can see below, in order to determine the control vector which linearizes the system. This diffeomorphism is given by:

$$\left(\begin{array}{c}{x}_0\\ c\psi c\theta .u+\left(-s\psi c\varphi +c\varphi s\theta \right)v+\left(s\psi s\varphi +c\varphi s\theta \right)w\\ y_0\\ s\psi c\theta .u+\left(c\varphi +s\varphi s\theta \right)v+\left(-c\varphi +s\varphi s\theta \right)w\\ z_0\\ -s\theta .u+s\varphi .c+c\varphi .c\\ \varphi \\ p+s\varphi .\tan \theta .q+c\varphi .\tan \theta .r\\ \theta \\ c\varphi .q-s\varphi .r\\ \psi \\ \frac{s\varphi }{c\theta }+\frac{c\varphi }{c\theta }r\end{array}\right)$$

(20)

Using the new coordinates, by applying the diffeomorphism, the system in Equation (12) can be written as follows:

$$\left(\begin{array}{c}{\dot{\xi }}_1\\ {\ddot{\xi }}_1\\ {\dot{\xi }}_2\\ {\ddot{\xi }}_2\\ {\dot{\xi }}_3\\ {\ddot{\xi }}_3\\ {\dot{\xi }}_4\\ {\ddot{\xi }}_4\\ {\dot{\xi }}_5\\ {\ddot{\xi }}_5\\ {\dot{\xi }}_6\\ {\ddot{\xi }}_6\\ \end{array}\right)=\left(\begin{array}{c}{L}_f{h}_1\\ L_f^2{h}_1+{\displaystyle \sum _{i=1}^3}{L}_{g_i}{L}_f^2{h}_1{\alpha }_i+L_{p_1}{L}_f{h}_1{F}_{v_1}\\ L_f{h}_2\\ L_f^2{h}_2+{\displaystyle \sum _{i=1}^3}{L}_{g_i}{L}_f^2{h}_2{\alpha }_i\\ L_f{h}_3\\ L_f^2{h}_3+{\displaystyle \sum _{i=1}^3}{L}_{g_i}{L}_f^2{h}_3{\alpha }_i\\ L_f{h}_4\\ L_f^2{h}_4+{\displaystyle \sum _{i=4}^6}{L}_{g_i}{L}_f^2{h}_4{\alpha }_i\\ L_f{h}_5\\ L_f^2{h}_5+{\displaystyle \sum _{i=4}^6}{L}_{g_i}{L}_f^2{h}_5{\alpha }_i\\ L_f{h}_6\\ L_f^2{h}_6+{\displaystyle \sum _{i=4}^6}{L}_{g_i}{L}_f^2{h}_6{\alpha }_i\end{array}\right)$$

(21)

By combining the equations containing the control vector α_i and Equation (16), it becomes:

$$\left(\begin{array}{c}{\ddot{\xi }}_1\\ {\ddot{\xi }}_2\\ {\ddot{\xi }}_3\\ {\ddot{\xi }}_4\\ {\ddot{\xi }}_5\\ {\ddot{\xi }}_6\end{array}\right)=\underset{\mathsf{ƛ}}{\underbrace{\left(\begin{array}{c}{L}_f^2{h}_1\\ L_f^2{h}_2\\ L_f^2{h}_3\\ L_f^2{h}_4\\ L_f^2{h}_5\\ L_f^2{h}_6\end{array}\right)}}+D\underset{\alpha }{\underbrace{\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\\ {\alpha }_4\\ {\alpha }_5\\ {\alpha }_6\end{array}\right)}}+F_{v_1}.\underset{\mathsf{\Omega }}{\underbrace{\left(\begin{array}{c}\frac{c\psi c\theta }{M_{11}}\\ \frac{s\psi c\theta }{M_{11}}\\ \frac{-s\theta }{M_{11}}\\ 0\\ 0\\ 0\end{array}\right)}}=\hspace{0.33em}\hspace{0.33em}\mathsf{ƛ}+D\alpha +F_{v_1}.\mathsf{\Omega }\hspace{0.33em}\hspace{0.33em}=\aleph$$

(22)

3.2.4. The Linearizing Control Law

According to Equation (22), we can deduce the linearizing control α. This control vector will ensure the decoupling and can be written as:

$$\alpha =D^{-1}\left(\aleph -\mathsf{ƛ}-F_{v_1}.\mathsf{\Omega }\right)$$

(23)

Note that the linearization would only be possible if the decoupling matrix D is invertible. However, the decoupling matrix D given by:

$$D=\left(\begin{array}{cc}\mathsf{\Pi }& 0\\ 0& \mathsf{\Gamma }\end{array}\right)$$

(24)

is invertible, and we have:

$$D^{-1}=\left(\begin{array}{cc}{\mathsf{\Pi }}^{-1}& 0\\ 0& {\mathsf{\Gamma }}^{-1}\end{array}\right)$$

(25)

(For the expressions of Π, Γ, ${\mathsf{\Pi }}^{-1}$, and ${\mathsf{\Gamma }}^{-1}$ see Appendix B).

Replacing Equation (23) in Equation (22), the equivalent system becomes linear and completely decoupled. It can be written as follows:

$$\left(\begin{array}{c}{\ddot{\xi }}_1\\ {\ddot{\xi }}_2\\ {\ddot{\xi }}_3\\ {\ddot{\xi }}_4\\ {\ddot{\xi }}_5\\ {\ddot{\xi }}_6\end{array}\right)=\left(\begin{array}{c}{\dot{z}}_2^1\\ {\dot{z}}_2^2\\ {\dot{z}}_2^3\\ {\dot{z}}_2^4\\ {\dot{z}}_2^5\\ {\dot{z}}_2^6\end{array}\right)=\left(\begin{array}{c}{\aleph }_1\\ {\aleph }_2\\ {\aleph }_3\\ {\aleph }_4\\ {\aleph }_5\\ {\aleph }_6\end{array}\right)$$

(26)

Finally, we obtain the following six decoupled linear sub-systems with their outputs:

$$\left\{\begin{array}{c}{\dot{z}}_1^i=z_2^i\\ {\dot{z}}_2^i={\aleph }_i\\ {\xi }_i=z_1^i\end{array}\right.$$

(27)

The vector $\aleph$ is designed according to the control objectives. In general, for a tracking trajectory problem, the expressions of ${\aleph }_i$ are:

$${\aleph }_i={\omega }_i^{\left(r_i\right)}+{\sigma }_{r_i-1}\left(y_i^{\left(r_i-1\right)}-{\omega }_i^{\left(r_i-1\right)}\right)+\dots +{\sigma }_1\left(y_i-{\omega }_i^{}\right)$$

(28)

${\omega }_{\hspace{0.17em}\hspace{0.17em}i}^{\hspace{0.17em}\hspace{0.17em}\left(\hspace{0.17em}{r}_{\hspace{0.17em}i}\right)}$ are reference trajectories. In our case these are functions in steps $\left({\dot{\omega }}_i^{\left(r_i\right)}={\ddot{\omega }}_i^{\left(r_i\right)}=0\right)$.

${\sigma }_{r_i-1}$ are the gains. If the ${\sigma }_i$ are chosen so that the polynomial $s^{r_i}+{\sigma }_{r_i-1}.s^{r_i-1}+\cdots {\sigma }_2.s+{\sigma }_1=0$, the polynomial is said to be a Hurwitz polynomial (possesses roots with negative real parts), then it can be shown that the error ξ_i = y_i − Ψ_i tends towards zero.

We have chosen the following poles P = {−10, −4} for the system in Equation (28) to be stable, then the gains of the control ${\aleph }_i$, which are solutions of the equation $s^2+{\sigma }_2.s+{\sigma }_1=0$, are given by: σ₁ = 40 and σ₂ = 14.

This values are needed in the numerical simulations to assess the performance of the proposed controller.

Previous developments show that it is possible to anticipate the effect of a gust of wind on the airship. The feedforward control that we have just established is a key part of this objective. Indeed, this control ensures the accelerated attenuation of the effect of the wind disturbance when this is known in advance, and it gives the control vector a certain robustness with respect to any modeling errors.

We have also demonstrated the following result:

Proposal. Consider the Multiple-Input Multiple-Output (MIMO) system of the airship MC500 described by Equation (12).

We define ${\left(r_j\right)}_{1\le j\le 6}$ as the relative degree of the outputs and ρ₁₁ as the relative degree of the disturbance F_ν1 associated with the outputs y₁ = x₀.

It is assumed that the decoupling matrix D is invertible, then the control vector U given in Equation (23):

• anticipates the effect of a gust of wind.
• ensures the stability of the airship in the vicinity of a target state.

Note:

• We have restricted ourselves to a single excitation force along the X axis to limit the size of the equations. However, to take into account the effect of other force and torque components, the same steps described above would be followed.
• The airship is operated by electric motors having sufficient degrees of freedom. As mentioned previously, the airship is over-actuated. It will thus be possible to anticipate the effect of all types of wind gusts and thus prevent against the main destabilizing element of the airship (the wind).

It should be remembered that the latter is the main factor in delaying the expansion of the use of the airships as flying cranes or as means of transport of heavy loads.

4. Control Allocation

We will designate the command established previously by virtual command because it acts macroscopically on the overall motion of the airship. However, in order to make these global setpoints achievable, it is necessary to couple them to the real actuators and to check that the latter are not saturated. This part is often obscured in the literature; however, we cannot ignore it in this study for the case of an airship having a very high inertia. We mentioned in Section 2 that the actuators of the airship MC500 are the thrust forces and the swivel angles of the propellers. As can be seen in Figure 5, to close the control loop, we need to express the non-linear control α_i in terms of the thrust values $‖F_{\hspace{0.17em}i}‖$ and inclination angle values ${\beta }_{\hspace{0.17em}i}$ and ${\gamma }_{\hspace{0.17em}i}$. It is therefore a question of control allocation.

As an input, we have a control vector with six components corresponding to the six degrees of freedom of the airship, and as an output we have the twelve degrees of freedom of the actuators. It would therefore be a question of solving a rectangular system which is done by means of numerical inversions, this can be penalizing for a real time control. We present, in what follows, an original analytical technique, based on energy concepts, which allows us to solve the previous problem quickly and precisely. Finally, we will compare the results obtained with those obtained from other known techniques.

We would like to emphasize that, in this study, we are dealing with the case of the airship when hovering. We will therefore neglect the aerodynamic forces of lift and drag.

However, the buoyancy B_u, which is an ascending vertical force opposing the weight (m.g) and which is one of the main characteristics of airships, is well taken into account. This force is applied on the center of volume of the airship, and has the coordinates (0, 0, z_B).

The relation between the actuators and the control α_i is described by the following rectangular nonlinear system:

$$\left\{\begin{array}{c}{\alpha }_1={\displaystyle \sum _{k=1}^4}∥F_i∥c{\gamma }_{{}_i}c{\beta }_{{}_i}-(mg-B_u).s\theta \hfill \\ {\alpha }_2={\displaystyle \sum _{k=1}^4}∥F_i∥s{\gamma }_{{}_i}+(mg-B_u).s\varphi .c\theta \hfill \\ {\alpha }_3={\displaystyle \sum _{k=1}^4}∥F_i∥c{\gamma }_{{}_i}s{\beta }_{{}_i}+(mg-B_u).c\varphi c\theta \hfill \\ {\alpha }_4=c{\displaystyle \sum _{k=1}^4}∥F_i∥s{\gamma }_{{}_i}+b_1(∥F_1∥c{\gamma }_{{}_1}s{\beta }_{{}_1}-∥F_2∥c{\gamma }_{{}_2}s{\beta }_{{}_2})+b_3(∥F_3∥c{\gamma }_{{}_3}s{\beta }_{{}_3}-∥F_4∥c{\gamma }_{{}_4}s{\beta }_{{}_4})+B_u{z}_{B}s\varphi .c\theta \hfill \\ {\alpha }_5=-c{\displaystyle \sum _{k=1}^4}∥F_i∥c{\gamma }_{{}_i}c{\beta }_{{}_i}+a(∥F_4∥c{\gamma }_{{}_4}s{\beta }_{{}_4}+∥F_4∥c{\gamma }_{{}_3}s{\beta }_{{}_3}-∥F_1∥c{\gamma }_{{}_1}s{\beta }_{{}_1}-∥F_2∥c{\gamma }_{{}_2}s{\beta }_{{}_2})+B_u{z}_{B}s\theta \hfill \\ {\alpha }_6=b_1(∥F_1∥c{\gamma }_{{}_1}c{\beta }_{{}_1}-∥F_2∥c{\gamma }_{{}_2}c{\beta }_{{}_2})+b_3(∥F_3∥c{\gamma }_{{}_3}c{\beta }_{{}_3}-∥F_4∥c{\gamma }_{{}_4}c{\beta }_{{}_4})+a(∥F_4∥s{\gamma }_{{}_4}+∥F_3∥s{\gamma }_{{}_3}-∥F_1∥s{\gamma }_{{}_1}-∥F_2∥s{\gamma }_{{}_2})\hfill \end{array}\right.$$

(29)

Remember that we have twelve actuators for the airship, while the control vector has six dimensions. The resolution of this nonlinear system admits an infinity of solutions.

To solve this connection problem, we introduce the following intermediate parameters for each rotor i:

$$l_i=‖F_i‖c{\gamma }_{i}c{\beta }_i; m_i=‖F_i‖s{\gamma }_i ; n_i=‖F_i‖c{\gamma }_{i}s{\beta }_i$$

(30)

The parameters of the actuators are deduced from these intermediate variables by the following relationships:

$$‖F_i‖=\sqrt{l_i^2+m_i^2+n_i^2}; {\beta }_i=\arctan \left(\frac{n_i}{l_i}\right) \mathrm{and} {\gamma }_i=\arctan \left(\frac{m_i}{c{\beta }_i.l_i}\right)$$

(31)

It would therefore suffice to express these intermediate variables as a function of the control vector in order to be able to deduce therefrom a relation between the actuators and the control.

By introducing the variables $l_i$, $m_i,$ and $n_i$ into the system (29), we obtain the following matrix system:

$$C.U_1={\alpha }^{\prime }$$

(32)

With:

$$C=\left(\begin{array}{c}1 1 1 1 0 0 0 0 0 0 0 0\\ 0 0 0 0 1 1 1 1 0 0 0 0\\ 0 0 0 0 0 0 0 0 1 1 1 1\\ 0 0 0 0 c c c c b_1 b_1-b_3-b_3\\ -c-c-c-c 0 0 0 0-a-a a a\\ b_1-b_1 b_3-b_3-a-a a a 0 0 0 0\end{array}\right)$$

(33)

$U_1={\left(l_1\dots l_4, m_1\dots m_4, n_1\dots n_4\right)}^T$

$${\alpha }^{\prime }=\left(\begin{array}{c}{\alpha }_1+\left(mg-B_u\right)s\theta \\ {\alpha }_2-\left(mg-B_u\right)s\varnothing c\theta \\ {\alpha }_3-\left(mg-B_u\right)c\varnothing c\theta \\ {\alpha }_4-B_u{z}_{B}s\varnothing c\theta \\ {\alpha }_5-B_u{z}_{B}s\theta \\ {\alpha }_6\end{array}\right)$$

(34)

We can use the pseudo-inverse method to find the values of, $l_i$, $m_i$, and $n_i$. However, this technique sometimes gives physically unacceptable values incompatible with the capacities of the actuators so that the stabilization of the airship is compromised. We have therefore decided to analytically solve the linear system in Equation (29).

4.1. Analytical Approach

The system in Equation (29), according to the new variables $l_i$, $m_i,$ and $n_i,$ can be divided into two subsystems that will be separately investigated:

$$\left\{\begin{array}{ccc}{\displaystyle \sum _{k=1}^4}{n}_i={{\alpha }^{\prime }}_3\hfill & \hfill & \hfill \\ b_1(n_1-n_2)+b_3(n_3-n_4)={{\alpha }^{\prime }}_4-c{{\alpha }^{\prime }}_2\hfill & \hfill & \hfill \\ a(n_3+n_4-n_1-n_2)=c{{\alpha }^{\prime }}_1+{{\alpha }^{\prime }}_5\hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill \end{array}\right.$$

(35)

$$\left\{\begin{array}{ccc}{\displaystyle \sum _{k=1}^4}{l}_i={{\alpha }^{\prime }}_1\hfill & \hfill & \hfill \\ {\displaystyle \sum _{k=1}^4}{m}_i={{\alpha }^{\prime }}_2\hfill & \hfill & \hfill \\ b_1(l_1-l_2)+b_3(l_3-l_4)+a(m_4+m_3-m_1-m_2)={{\alpha }^{\prime }}_6\hfill & \hfill & \hfill \end{array}\right.$$

(36)

The first and third equation of the system in Equation (35) give:

$$n_3+n_4=\frac{{{\alpha }^{\prime }}_3}{2}+\frac{c{{\alpha }^{\prime }}_1+{{\alpha }^{\prime }}_5}{2a}$$

(37)

We initially imposed this choice:

$$n_3=n_4=\frac{{{\alpha }^{\prime }}_3}{4}+\frac{c{{\alpha }^{\prime }}_1+{{\alpha }^{\prime }}_5}{4a}$$

(38)

Substituting h₃ and h₄ by their values in the system in Equation (36), one can obtain:

$$n_1=\frac{{{\alpha }^{\prime }}_3}{2}+\frac{{{\alpha }^{\prime }}_4-c{{\alpha }^{\prime }}_2}{2b_1}-2n_3$$

(39)

Additionally, then:

$$n_2={{\alpha }^{\prime }}_4-n_1-n_3$$

(40)

By combining the second and the third equation of the system in Equation (36), this gives:

$$b_1(l_1-l_2)+b_3(l_3-l_4)+2a(m_4+m_3)={{\alpha }^{\prime }}_6+a{{\alpha }^{\prime }}_2$$

(41)

We decided to add these two conditions according to the logical operation of the actuators:

$$\mathrm{and}: m_3=m_4=\frac{{{\alpha }^{\prime }}_2}{4}$$

(42)

$$b_1(l_1-l_2)=b_3(l_3-l_4)=\frac{{{\alpha }^{\prime }}_6}{2}$$

(43)

These assumptions allowed us to retrieve the expressions of the missing variables through the following relations:

$$\left\{\begin{array}{l}{m}_1=m_2=\frac{{{\alpha }^{\prime }}_2}{4}\hspace{1em};\hspace{1em}{l}_2=l_4=\frac{{{\alpha }^{\prime }}_1}{4}-(\frac{1}{8b_1}+\frac{1}{8b_3}){\alpha }^{\prime }\\ \\ l_1=\frac{{{\alpha }^{\prime }}_6}{2b_1}+l_2\hspace{1em};\hspace{1em}{l}_3=\frac{{{\alpha }^{\prime }}_6}{2b_3}+l_4\end{array}\right.$$

(44)

The block diagram shown in Figure 6 describes the original algorithm established to define the relationship between the actuators of the MC500 and the controllers. The algorithm is based on algebraic relations which will be solved analytically.

The proposed approach is based on energy principles. Indeed, and in the case of the hovering flight of a propeller-driven quadcopter, there are several solutions. However, the solution that induces the minimum energy consumption is the one that balances the thrust of the four rotors. To this goal, we have imposed some constraints to be as close as possible to the equilibrium configuration of the different thrusts whenever the operating conditions allow it.

As a comparative analysis, we present in the next paragraph a numerical approach establishing the control allocation.

The comparison of the results of these two methods will be performed in the chapter devoted to numerical simulations.

4.2. Numerical Approach

To check this approximate approach, we compare our results here with a numerical method. We used the fixed-step gradient method to compute the solution numerically. This algorithm is iterative based on the minimization of a cost function $J\left(U_1\right)$.

The principle of the algorithm is to generate a vector series $U_{1,k}$ from an arbitrary point $U_{1,0}$ such that the cost value of the function $J\left(U_1\right)$ decreases with each iteration, i.e.,:

$$J\left(U_{1,k+1}\right)<J\left(U_{1,k}\right)$$

For the fixed-step gradient method, the vector $U_{1,k}$ is updated in this way:

$$U_{1,k+1}=U_{1,k}-\frac{1}{2}\mu \frac{\partial J\left(U_{1,k}\right)}{\partial U_1}$$

(45)

The index k denotes the iteration and $\mu$ is a positive constant parameter.

To solve the linear system in Equation (32) we propose the use of the cost function J defined by:

$$J\left(\tau \right)={\displaystyle \sum }_{i=1}^6{\left({{\alpha }^{\prime }}_i-{\tau }_i\right)}^2$$

(46)

The terms of the quadratic function penalize the error $e_i={{\alpha }^{\prime }}_i-{\tau }_i$ between the control vector, the vector forces ${\tau }_1,$ and the torques ${\tau }_2$(see Equation (29)).

We can see that when the quadratic function is minimal, the terms $e_i$ tend toward zero, then ${\tau }_i$ tend to ${{\alpha }^{\prime }}_i$.

Hence, the intermediate parameters obtained, l_i, m_i, and n_i, are the components of the vector that minimize the cost function.

By using the matrix relations in Equation (32)–(34), the cost function J can be written as:

$$J\left(U_1\right)={\displaystyle \sum }_{i=1}^6{\left({{\alpha }^{\prime }}_i-C_i{U}_1\right)}^2$$

(47)

With $C_i$ is the ith line of the matrix C.

The gradient of J is:

$$\nabla J\left(U_1\right)=-2{\displaystyle \sum }_{i=1}^6{C}_i^t\left({{\alpha }^{\prime }}_i-C_i{U}_1\right)$$

(48)

Then, Equation (45) becomes:

$$U_{1,k+1}=U_{1,k}+2\mu {\displaystyle \sum }_{i=1}^6{C}_i^t\left({{\alpha }^{\prime }}_i-C_i{U}_{1,k}\right)$$

(49)

This method has the advantage of being easy to implement. Unfortunately, the convergence of this method is generally slow. The steps for obtaining the numerical solution are detailed in Algorithm 1:

Algorithm 1 fixed-step gradient

Step 1: Initialization k = 0: Choice of $l_{i,0},m_{i,0},n_{i,0}$ and μ, ε > 0 Step 2: Iteration k $\left(i=1,\cdots ,4\right)$ ${\left[l_{i,k+1},m_{i,k+1},n_{i,k+1}\right]}^t={\left[l_{i,k},m_{i,k},n_{i,k}\right]}^t+ 2\mu {\displaystyle \sum _{i=1}^6}{C}_i^t\left({{\alpha }^{\prime }}_i-C_i{\left[l_{i,k},m_{i,k},n_{i,k}\right]}^t\right)$ Step 3: Stop criterion: $‖{\left[l_{i,k+1},m_{i,k+1},n_{i,k+1}\right]}^t,-,{\left[l_{i,k},m_{i,k},n_{i,k}\right]}^t‖<\epsilon$

Note:

• In order to implement the proposed algorithm, the convergence of it is an important factor to take into account. The gradient algorithm that minimizes the quadratic function in Equation (46) converges if: $0\prec \mu \prec \frac{2M_2}{M_1^2}$ with $M_1=\max \left({\lambda }_i\right); M_2=\min \left({\lambda }_i\right)$, Where ${\left({\lambda }_i\right)}_{1\le i\le 12}$ are the eigenvalues of the symmetric matrix $C_1=2{\displaystyle \sum _{i=1}^6{C}_i^T{C}_i}$.

For both the analytical and the numerical methods, we have succeeded in stabilizing the airship. However, during a short transient phase, the actuators’ capacities were exceeded. In this first approach, we have introduced a saturation block assuming that the propellers are in saturation abutment. Further studies are underway to solve this problem using constrained optimization techniques. This will be dealt with in future work.

5. Numerical Results

We present here some examples demonstrating the utility of the proposed formulation in different cases.

For this, we used data from the MC500 airship, such as:

The volume V = 500 m³; z_G = 0.5 m; a = 2.5 m; b₁ = 5.4 m; b₃ = 6.5 m; c = 2 m;

The added mass matrix and the inertia are given by:

M₁₁ = 583 kg; M₄₆ = 160 kg.m²; µ = 0.01; ε = 10⁻³; M₂₂ = 620 kg; M₃₃ = 687 kg; M₄₄ = 9413 kg.m²; M₅₅ = 10,456 kg.m²; M₆₆ = 18,700 kg.m²;

5.1. Ideal Case: No Wind Disturbance

In the first step we check the efficiency of the control vector and the performances of the proposed algorithms to solve the connection problem envisaged under ideal conditions, i.e., the wind speed is zero.

As a reference example, we have chosen to deal with the scenario of the inclination of the airship in roll and yaw and its displacement with respect to an original position. Here, we do not consider the effect of the wind. Only the feedback control comes into action.

We can see in Figure 7 and Figure 8 that the control manages to stabilize the positions x₀, y₀, z₀ and the angles ϕ, θ, ψ of the initial state X_e = [10; 8; 5; 0.3; 0.2; 0.1] to the desired state X_d = [0; 0; 0; 0; 0; 0].

We can also see that the stabilization was carried out in a quick manner, less than 2 s for the analytical method, and between 4 and 10 s for the gradient method. This confirms the power of the analytic algorithm developed, not to mention the time saving CPU compared to that of the numerical method.

The behaviors of the actuators in response to this command are shown in Figure 9, Figure 10 and Figure 11. This concerns the thrust forces F_i, the tilt angles β_i, and the swing angles γ_i. It will be noted that, after the stabilization of the airship, the thrust is equitably distributed over the four rotors, the swing angles γ_i become equal to zero, and the tilt angles β_i are equal to 90°. This effectively characterizes the optimum stationary stability position. This therefore highlights the interest of the analytical methodology followed.

5.2. System Analysis in the Presence of Wind

In this section, we consider two scenarios of the airship expecting acceleration due to a gust of wind pre-detected by a LIDAR. In one case, the wind is regular, and in the other, the wind is turbulent (sinusoidal). We will apply the combined feedforward/feedback control developed in order to minimize the drift of the airship under the effect of this gust of wind and thus be able to save a lot of energy necessary for the displacement of the airship due to the great inertia of the latter. This also helps to minimize the possible sway of a load suspended by the airship.

The wind is supposed to arrive along the X axis and produce an average force of 10⁴ N (see Figure 12).

It can be seen in Figure 13 that the airship is subject to a large drift with respect to its equilibrium position if it is subjected solely to feedback control (red dash). However, with the application of the feedforward/feedback control previously developed, the reduction of this drift is well demonstrated, especially in the case of the use of the analytical algorithm of transformation of the controls to the actuators (blue dash).

To make the experiment more realistic and to demonstrate the robustness of our algorithm, we have developed our control on the basis of a constant wind step with an intensity of 10⁴ N, whereas, for the force applied to the airship, we allowed a sinusoidal variation of the wind intensity (see Figure 12). In spite of the relatively large difference between this applied force and the force on which the feedforward/feedback control was based, we note (in Figure 14) that our feedback control is efficient.

The response of the actuators to the disturbance due to the wind is just as effective as in the previous scenario when the feedforward control is applied, in particular, when the analytical algorithm is used. The airship stabilizes quickly after having undergone low amplitude oscillations. This drift is even less when using the analytical control allocation algorithm.

We must emphasize that the established control is not considered to be robust, and can suffer from modeling errors. We tested its robustness through the modification of the terms of the mass matrix which were increased by 50%. The results are presented in Figure 15. It is noted that the control no longer manages to stabilize the flying machine with such deviations in the modeling. Concretely, this control remains valid as long as the errors in the mass matrix do not exceed 30%.

6. Discussion

The different scenarios presented here have proven the effectiveness of the feedforward control to support feedback control. The latter, although being effective for stabilizing the airship in the case of small disturbances (scenario 1), shows its limits during strong gusts of wind. It was shown in scenario 2 (constant strong wind) or scenario 3 (strong wind with sinusoidal amplitude) that the help of the feedforward control developed makes it possible to greatly attenuate the airship’s drift if the disturbance is precisely predicted (scenario 2) or hovers around a predicted value (scenario 3).

The analytical control allocation algorithm developed also contributed to the improvement of the stabilization of the machine compared to purely numerical algorithms of reference.

7. Conclusions

As part of the development of large airships, we have contributed, through this paper, to the modeling and control of an unconventionally shaped airship. These airships are designed to load and unload freight at altitude. Our objective focused on the stabilization of the flying machine in this critical phase of hovering under the effect of an exogenous disturbance, and especially under the effect of a gust of wind. We then established a control algorithm based on a feedforward/feedback technique. The use of the feedforward control was made possible by the use of LIDAR, which allow for anticipating the amplitude and the speed of the wind gusts. The combined feedforward and feedback linearizing control developed is based on the differential geometric analysis and the relative degree.

The great sensitivity of these flying machines to the wind and their great inertia has pushed the designers of the airship studied during this work to multiply the actuators and their degrees of freedom. This over-actuation has complicated the control allocation problem. For this, we have established a fully analytical algorithm with some assumptions for solving this problem. Our algorithm has proved its efficiency in comparison with other classical numerical methods, and especially its superiority in the face of these same algorithms in terms of minimizing the computation time.

In future works, the feedforward control vector algorithm will be applied to the trajectory tracking.