Adaptive tracking control of flapping wing micro-air vehicles with averaging theory

: An input constrained adaptive tracking controller is designed for flapping micro aerial vehicles, wherein the moving averaging filter is adopted to estimate the averaged states of the system. Specifically, in the outer loop controller, an observer is constructed to estimate the disturbances within the system. Moreover, the constrained thrust is designed to keep the frequency in a proper region so as to meet the requirement of average estimation. Then, a tracking differentiator is used to provide trackable trajectories for the inner loop. Subsequently, a new quaternion-based hybrid attitude tracking controller is designed which successfully deals with high-frequency noises and avoids possible chattering. As supported by mathematical analysis, the proposed control strategy guarantees the uniform ultimate boundedness of the closed-loop system, and it keeps the control torques within the permitted range to meet the application requirement. At last, numerical simulations are carried out to support the validity of the proposed controller, whose results are satisfactory even when the thrust and torques are saturated.


Introduction
In the size of insects and birds, flapping flight possesses amount of advantages [1,2]. Based on this observation, many flapping wing micro vehicles have been invented till now, including Mentor flapping-wing micro air vehicle (FMAV) [3], Delfly [ 4], Nano Hummingbird [5], Harvard Robotic Fly [6] and so on. To make a flapping wing micro vehicle fly with much flexibility, it is of key importance to design and implement a suitable control for the system. In recent works, many control algorithms, such as model-free control [7], adaptive control [8,9], geometric control [10], and some other control methods [11], have been successfully introduced for the control of different flapping wing vehicles.
In the past decade, the problem of stabilising the vehicle at the desired attitudes has long been a popular research topic for the control of flapping wing vehicles. Lin et al. use a modified P-control to stabilise the altitude of their Golden Snitch [12]. Banazadeh et al. design an adaptive sliding mode technique for the attitude and position control of a rigid body, insect-like flapping wing model in the presence of uncertainties [13]. An adaptive tracking controller is proposed in [14] for the Harvard fly. A common feature of these methods is that much effort has been focused on the attitude control of the studied FMAV. In fact, once the fast changing inner loop control (attitude control) is determined, the outer loop control (position control) can be designed according to practical applications with much flexibility. For the position control of an FMAV system [13,14], averaging method proves to be a mainstream way to construct practical controllers [11,15], which approximates the orbit of the oscillating system by the solution of an average system, resulting in more straightforward analysis for the system's stability [16]. Averaging theory seeks a way to approximate the non-autonomous vector field with an autonomous one. That is, with the help of averaging, one may use a non-linear time invariant system to describe the intractable non-linear time varying system [16]. The theory generally applies to numerous classes of time-dependent vector fields, instead of only strictly periodically varying vector fields [17]. More details about averaging theory could be found in [17,18]. However, to utilise averaging method, the frequency of the wing flap needs to be high enough and the flaps' effects need to be relatively stable. Another issue with averaging method is that the averaged vector field cannot be computed explicitly for an FMAV system in practical applications [16,19]. In this paper, we seek no particular autonomous approximation. Oppositely, using a moving average filter, the original time-dependent vector field of FMAV is regarded as a non-autonomous vector field.
A main difficulty for flapping wing tracking is the unavoidable vibration, which is apparently reflected in the sense that both the forces and torques are oscillating. Also, practical flapping wing systems are intrinsically subjected to various input constraints. In fact, to facilitate the application of the averaging method, it is even necessary to set a lower bound on the flapping frequency. Therefore, how to propose a suitable controller to satisfy those constraints is another challenging issue. Considered the aforementioned problems and inspired by the work reported in [20,21], this paper proposes an adaptive tracking control method for an FMAV system, which introduces averaged states as feedback signals and applies some robust auxiliary terms to handle the drifts or disturbances within the system. Meanwhile, considering that the measurement is always oscillating, the controller utilises the estimation generated from a reliable average estimator to construct the controller. To avoid singularity, some specific saturation function has been introduced and the quaternion-based description [22,23] is adopted to construct attitude filters. The performance of the proposed adaptive tracking controller is verified by both rigorous analysis and simulation results.
The contribution of this paper mainly consists of the following three parts: (i) The proposed adaptive position tracking controller guarantees satisfactory control performance even in the presence of various disturbances and input constraints; (ii) a new quaternion-based attitude tracking controller is designed, which is free of singularity and successfully keeps the tracking error within a given bound; (iii) a systematic procedure is proposed, which designs control laws for the outer loop and the inner one separately, and thus brings much flexibility and increases the robustness of the overall system.
The main body of this paper falls into the following four sections: in Section 2, we briefly analyse the dynamics of the flapping wing vehicle and provide the dynamic model of the system. Next, the control laws are designed in Section 3 in details. Stability analysis for the closed-loop system is performed in Section 4. Finally, some simulation results, with corresponding analysis, are presented in Section 5.

System dynamics analysis
It is worth stressing that flapping wings cause non-negligible transient effects on the vehicles' dynamics. For most practical designs, the vehicle's subcycle movement is neither controllable nor observable, and only macro parameters, for example, flapping frequency, can be adjusted relatively slowly. Parametrisation method is usually used to model the system, which conveniently sets up the relationship between the adjustable parameters and the forces/torques actuated on the vehicle. Later on, these parameters can be then regarded as control inputs to change the state of FMAVs. For example, Deng et al. use this method to parameterise the wing trajectory and then propose a flight controller for biomimetic robotic insects in [24].
Inspired by [24], also referring to [25], the dynamic equation for the studied FMAV system is given as In (1), X [ R 13 refers to the system states, including the position of the vehicle's mass centre p [ R 3 in the inertia coordinate, the speed of the vehicle's mass centreṗ [ R 3 in the inertia coordinate, the unit quaternion q [ Q representing the attitude of the vehicle in the inertia coordinate and the angular rate B v [ R 3 in the body coordinate. The signal g represents gravitational acceleration, e 3 = 001 T , m and J stand for the mass and the inertia matrix of the vehicle. I B R [ SO(3) denotes the rotation from the vehicle body coordinate to the inertia coordinate, which can be calculated according to q. The definition of quaternion multiplication operator ⊗ can be found in [23]. Steady aerodynamic damping force F ad (·) [ R 13 R 3 and toque t ad (·) [ R 13 R 3 depend on the states, which are parameterised as linear functions. Note that R 3 determine the fast oscillating behaviours of the system, where 1 denotes a small positive constant, and F th [ R 3 , t r [ R 3 represent the normalised thrust input and torque input, respectively. In addition, F p and t p are zero mean periodic functions if F th , t r and X are stationary.
Assumption 1: Bounded actuation torque t r is supposed to be uniformly omnidirectionally adjustable around the origin, provided that F th remains higher than a certain magnitude. This assumption is actually reasonable since it only implies that the vehicle's attitude is locally controllable.  [18] presents that the drift between the averaged system and the original one can be offset by a time-varying function of the order O d 1 1 () / 1T () as 1 0. Although it is not practical to determine the exact form of the offset, there do exist bounded choices of the offset which keep the performance of the two systems sufficiently close. In this sense, the offset can be regarded as a part of the disturbances, which will be then addressed in the subsequently designed adaptive tracking controller.

Remark 2: (Control inputs):
Since the relationships between the parameters and the average effect of input thrust are assumed known, F th and t r are directly taken as the inputs of the system, instead of the flapping frequencies and the rudders' steering angles.

Remark 3: (aerodynamic forces torques):
The aerodynamic forces and torques are divided into three pairs: (i) the steady aerodynamic damping F ad and t ad , (ii) the fast oscillating terms F p and t p , (iii) the steady input thrust F th and torque t r . The steady aerodynamics and control inputs present the same form in the following averaged systems, while the averaged effects of F p and t p correspond to the drifts taken as components of disturbances.

Control strategy development
For the studied FMAV system, the structure of the closed-loop system is demonstrated in Fig. 1. As can be seen, the states are filtered by a moving average filter, which estimates the average of the oscillating states. The filtered states are then utilised by both the inner loop controller and the outer loop one to calculate proper input for the FMAV system, wherein the outer loop controller directly determines the thrust input, as well as the desired orientations q d by integrating the yaw signal. Subsequently, the quaternion filter between the two controllers translates the desired orientations q d to a smooth trackable attitude trajectory q f . Thereafter, the constrained torque is produced by the inner loop controller, which enables the vehicle to track the generated trajectory. Along these lines, the control inputs F th , t x , t y and t z vary slowly since the average estimated states are employed as feedbacks, even though the flapping wing vehicle actually vibrates with high frequency. In general, by virtue of the average estimation of the states and the averaging method applied to the system, the complicated FMAV tracking task can be split into the following two sub-tasks: (i) determining the proper thrust direction; (ii) efficiently adjusting the controller when the thrust is relatively far from the desired one, with the input constraints and system disturbances successfully addressed, both of which can be fulfilled by the subsequently designed control strategy.

Outer loop control design
After the first-order averaging and rewriting x 1 as p, and x 2 asṗ, the state space representation of the outer loop dynamics becomeṡ (2) where the diagonal velocity matrix is defined as is the uncertain parameters vector of the translational normalised aerodynamic damping.
To keep the control torque effective and to satisfy the prerequisites of averaging method, it is necessary to set the upper and lower limits for the normalised thrust input u t [ R 3 as u t min ≤ u t ≤ u t max , where u t min and u t max are both positive constants. Since the original oscillating states are replaced by the relatively steady averaged states, the drift of the averaged system or the estimation error will occur unavoidably. To deal with that, we use resultant force disturbance d T [ R 3 to describe the deviation between slow varying part of the oscillating real system and the averaged system, which consists of three components, namely the averaged system drift, the averaging estimation error and the unknown disturbances (see also Remarks 1 and 3 for more details).

Assumption 2:
In order to make the desired trajectory trackable, we assume that the trajectory, as well as its first and second derivatives, is continuous, bounded and known. Moreover, the desired yaw angle and its first and second derivatives are also assumed to be continuous, bounded and known. To facilitate subsequent description, we define the desired trajectory as x 1d , and its derivative as x 2d , and denote the tracking errors as Assumption 3: As indicated by practical systems, the uncertain parameters vector u a is bounded, in the sense that u aj ≤ u aj ≤ū aj , for 0 , u aj ,ū aj with j = x, y, z.
Assumption 4: Without loss of much generality, the resultant force disturbance d T is assumed continuous and derivable.ḋ T is assumed within a relatively small bound in the sense thaṫ After some deep analysis for the system dynamics and referring to [21,26], it is concluded that endowing the saturated controller with the robustness has the equivalent importance as using the dynamics model information and instant feedback does, thus, the control input is designed in the following manner: with positive definite, diagonal matrices K z 1 [ R 3×3 and K z 2 [ R 3×3 as the control gains, where h [ R 3 and z [ R 3 are both auxiliary terms for input constraint and system robustness, respectively. In (3), to ensure that the controller is always within the permitted range, the following three-dimensional saturation function Sat T (·) [ R 3 R 3 is constructed: where a T [ R 3 , u ts , u t min and u t max are positive constants.
The update law of parameter estimationû a is designed aṡ where G d T [ R 3×3 and G u a [ R 3×3 are positive definite, diagonal gain matrices, and the projection function, introduced to keep the estimates within the predefined range, is explicitly formulated as Instead of directly tracking unobservable d T , it is convenient to track the modified observable variant d temp with compensation in other terms, since the deviation between d temp and d T is relative to the parameter estimation error. The specific observer for the resultant force disturbance d T is designed aṡ where and 1 z is a positive definite, diagonal matrix, and J d T takes the form: Note that the last term in (7), related with J d T , is introduced to cancel the intractable absolute term in the Lyapunov candidate derivatives, which will be shown in the subsequent stability analysis part. The dynamic equation of the robust auxiliary term z is defined aṡ where the gain matrix K z [ R 3×3 is positive definite and diagonal. The update law for the auxiliary term h is designed to mitigate the influences caused by past time force constraints: The desired orientation q d is computed according to the thrust orientation I t = u t / u t and the desired yaw angle c d , with the specific method formulated as where Remark 4: (Robustness enhancement): The term z detailed in (10) is seemingly redundant, yet it filters the higher frequency components of the observerd T , which is mostly aroused by the termẋ 2 in d temp .
Due to the existence of term 1 −1 z K z 2 (K z 1 z 1 + K z 2 z 2 ), z plays a similar role as the integration component in traditional PID control.

Inner loop control design
For inner loop control, the following quaternion-based hybrid filter is used to provide smooth attitude trajectory with the illumination of Mayhew's work [23], which relaxes the requirements on the desired trajectory stated in Assumption 1: where m qd1 m qd2 m qd3 m qd4 denoting a matrix function of unit quaternion: where p * stands for the conjugation of p. In (15), k qd is a positive constant and K v fd refers to a positive definite, diagonal matrix, s h (·) is a hysteretic function, which will be detailed below. Based on the attitude dynamics (1) and Assumption 1, after using the first-order average and linearising the aerodynamic damping, the dynamics of B v is calculated as where the torque control gain g T · () is a known monotonously increasing function, and u ta [ R 3 and u t× [ R 3 are uncertain parameters of the normalised damping and the normalised reduced inertia terms, and d t [ R 3 is the total normalised torque disturbance. The inertia matrix is assumed to be diagonal, due to the reason that most FMAVs are designed symmetrical. Moreover, u t [ R 3 stands for the normalised torque input.

Remark 5: (Torque control gain):
In (18), the term g T u t is assumed to be known, while the damping and inertia terms are uncertain. For a practical system, this assumption is reasonable since the torque control gain g T u t can be obtained for a fixed platform, yet the damping and the inertia terms cannot be determined conveniently as they are also influenced by various aerodynamic factors.
Assumption 5: For the studied system, the signals u ta , u t× and d t are assumed bounded as 0 , u taj , u taj ≤ū taj , where j = x, y, z.
With the motivation to obtain smooth desired angular rates, and to avoid unexpected effects of hybrid system jump, the compensation term q j is introduced with the inspiration of command filter [27], which is generated via the following manner: with the positive constant k q j and positive diagonal matrices as the compensator gains. The term d m0 refers to a small positive constant. Additionally, q js and q jv are the scalar and vectorial parts of quaternion q j . It is important to avoid chattering occurring around the point q T f q j ⊗ q = 0, because the real system is oscillating and prone to fall into the chattering range. Inspired by the hysteretic fashion presented in [23], the virtual angular speed a v is constructed as in which m qf 1 m qf 2 m qf 3 m qf 4 T W q * ⊗ q * j ⊗ q f , and m qf W m qf 2 m qf 3 m qf 4 T , and k q , as a positive constant, denotes the virtual angular velocity control gain. The term s h a ()is a one-dimensional hysteretic function as shown in Fig. 2, whose upper half and lower half are s h1 a ()and s h2 a () , respectively. The complete form of the function is shown as Item c i , i = 1, 2 is the function-choosing flag, introduced to identify which part of the hysteretic function the system is currently located on.c is the flag c in last computing cycle. The positive constant d a determines the width of the flat area, which is clearly shown in Fig. 2.
As mentioned previously, we use a second-order filter to shield the negative effects aroused by the hybrid system jumps: Similar as (3), the normalised torque is used to provide the control input with the specific composition formulated as u t = Sat t u to , where the positive diagonal matrix K 4 [ R 3×3 denotes the control gain. Note that ds h is the derivative of the function s h1 with respect to a, if the flow is on the first continuous part. After the jump, the term ds h becomes the derivative of s h2 with respect to a. Thus, there is no infinitive point for the function ds h , even though the function s presents two jumps. For practical reasons, the normalised torque constraint is considered, and the specific saturation function Sat t (·) [ R 3 R 3 takes the following form: The update laws of the uncertain parameters are designed aṡ where b ta , b t× , b dt , G ta , G t× and G dt , all belonging to R 3×3 , are positive diagonal gain matrices, and the projection function follows the previous definition in (6) with bounds shown in Assumption 5. Note that every entry of b dt should be larger than corresponding entry of G −1 dt , which enhances the convergence. The observer is not used for inner loop control because the direct angular acceleration is unobservable.
The term h t [ R 3 in (30) is used to alleviate the influences caused by past time constraint of normalised torque. Similar to (11), the update law of h t is constructed aṡ where s ht is a small positive constant, k h t refers to the positive gain and Du t = Sat t u t0 − u t0 .

Stability analysis
This section presents the stability analysis for the proposed control strategy.
Theorem 1: For the flapping wing vehicle governed by (1), the control law given in (3), (4) and (30), together with the desired attitude (12), ensures that the tracking errors of the position and yaw angle are uniformly ultimately bounded.
Proof: For the sake of clarity, we will present the proof in two parts, namely for the inner loop controller and the outer loop one, respectively. Analysis for the inner loop: To start with, as mentioned in (15), the filter translates the desired orientation to a continuous attitude trajectory as q f (t). Essentially, it is a generator for the desired trajectory. Therefore, it may be more comprehensible to treat the filter as a tracking differentiator of the signal q d , rather than a part of the inner loop controller. The convergence of the filter can be verified in a similar method validating the convergence of (26).
In accordance with the virtual angular speed definition (26), choose the Lyapunov candidate for the first step as From (26) and (36), we know that V q is positive when {} , and V q reaches its maximum right before its jump at q T f q j ⊗ q = d a or . d m0 , otherwise, the actual attitude is in the neighbourhood of the desired attitude.
For the sake of brevity, we only present the analysis for the situations when the system locates on s h1 as an illustrative example. Bearing the definition of the hysterical function (27) and (28) in mind, the stability analysis for the system with the proposed virtual control input can be distributed into three parts.
whose derivative iṡ (38) whose derivative can be calculated aṡ where the fact q T f q j ⊗ q [ −1, 1 [] has been utilised. Choosing its derivative can be obtained aṡ In subsequent analysis, we will consider the jump effect on the system's convergency. Due to the symmetry of the hysteretic function, two jumps generate the same change of V q , which is DV q =−(3/2)d 2 a . Therefore, V q is decreasing even when the system goes through the discontinuous jumps. Owing to the filter, the jump effect does not appear explicitly in the following parts of the system. Consequently, there is no need to apply special analysis for the corresponding Lyapunov functions. Moreover, the analysis is similar with the former analysis when the system is on s h2 , hence, it is omitted here. According to (40) and (42), it can be shown that V q reaches the region V q , 1 − d 2 a in finite time t , (d 2 a /2k q ) + (d 2 a /C q2 ) when the initial condition satisfies q T f q j ⊗ q 0 [ −d a , d a , and then, as indicated by (38), V q converges to the region where m qf 2 ≤ d m0 , exponentially fast, with adjustable bound d m0 defined in (25). The bound d a can be selected to meet different robustness requirements, as long as it is smaller than 2 √ /2.
Using the energy function as the Lyapunov candidate for the second part, we have The parameters estimation errors are defined asũ ta =û ta − u ta , u t× =û t× − u t× andd t =d t − d t , respectively. With the proposed update laws (32)-(35) and control laws (30) and (31), the derivative of (43) transforms intȯ The second term on the left side of (44) is the un-compensated part of last step. Based on previous analysis, we can conclude that g T u t /s 2 ht and g T u t are all positive, h t keeps small in the region h t ≤ s ht , only when the angular velocity error v d − v remains relatively small, since its derivative is bounded. In addition, v is bounded due to the damping term in (18), the bounded torque input and the bounded torque disturbance. Relatively small angular velocity error indicates that the virtual angular speed a v can be obtained with tolerable small errors. Then, we treat the inner loop system as the following hybrid system to obtain thorough comprehension: With the definition of the hybrid system and the choice of the closed sets, the system satisfies assumptions (A0)-(A3) in [28] with inspection. The conditions (VC) and (VD) in [28] are also satisfied accordingly, which ensure the existence of the hybrid system's solution. In this circumstance, q f is bounded and smooth enough to maintain the viability of the flow set. After checking the derivatives of the Lyapunov candidates, we can conclude that the largest invariant set is included in where d m0L and s htL are positive constants satisfying d m0L . d m0 and s htL . s ht . Therefore, according to Corollary 7.7 in [28], the system asymptotically converges to the largest invariant set in W, implying that the attitude tracking error, including the yaw angle error, is uniformly ultimately bounded. In addition, the chattering problem analysis is similar with the one presented in [23], from which we know that the hysteretic function successfully inhibits the chattering at the unstable equilibrium point of q T f q j ⊗ q = 0. Therefore, as long as the filter and the inner loop controller perform fast enough, and the initial orientation error is bounded, the orientation error will be uniformly bounded in any arbitrarily small bounds. This fact states that the virtual control forces are obtained instantly, and the unachieved part can be regarded as a part of force resultant disturbances. Although with different manners, this kind of cooperation problem between inner loop controller and outer loop one is also considered in [29,30].
Analysis for the outer loop: After rigorous inner loop controller stability analysis, it is straightforward to analyse the stability of the outer loop controller.
Based on the outer loop controller design and the system characteristics, we choose the Lyapunov candidate as The coefficients 1 z , K z 1 , K z 2 , G ua and G d T are all positive definite, diagonal matrices. The variabled T is the disturbance observer, and the uncertain parameter estimation error is defined as u a =û a − u a . Differentiating both sides of (48), the derivative of the Lyapunov candidate (48) is obtained: where According to (5) and the inequality the term − (ũ T a b uaũa )/2 + (u T a b ua u a )/2 can be deduced after substitutingũ T a G −1 uau a in (49). When decomposing d T −d T into d T − d temp and d temp −d T and conventionally magnifying the right side of (49), the intractable term d tempj −d Tj G −1 dTjḋTj appears, where j = x, y, z. Recalling (9), we introduce the following inequation to substitute this term with a manageable one. For ∀1 . 0, and ∀g [ R where k p is a constant satisfying k p = e −(k p +1) , whose approximation is chosen as 0.2758, which also meets the inequation (52). Therefore, there exists the following inequality: which facilitates the aforementioned substitution and the subsequent cancelation.
According to Chapter 4 in [31], regarding ge 3 + u t + d T as the input of the system (2), then due to the existence of the damping term, the system B vu t , d T is input-to-state stable, implying that B v is uniformly ultimately bounded with bounded u t and d T . Therefore, based on Assumptions 2 and 4 and the above analysis, there exist positive constants C ua , C k1 and C d T satisfying For clarity, we perform the substitutions and transformations of the derivative (49) in two steps. Firstly, based on the proposed update laws (5)-(10) and the control law (3), also using (53) and (54), we havė The second step considers the update law of h. Using the update law given in (11), the derivative becomeṡ On the other hand, we select the positive constant k Vo such that According to the previous analysis, we obtain the following inequation:V which states that where V o0 refers to the initial condition, and C o /(2k Vo ) indicates the final convergence value. Therefore, lim t 1 V o t ()≤(C o /(2k Vo )).

Simulation results and analysis
Simulations are carried out in order to verify the validity of the proposed controller. According to the FMAV in [3], we set the system parameters as J = diag 0.03, 0.03, 0.02 () kg m 2 , m = 0.5 kg, u a = 110 .4 T s −1 and u ta = 110 .5 T s −1 .
The main controller gains are chosen as The way of tuning controller gains follows a basic guideline that the quaternion filter and the inner loop controller performs much faster than the outer loop one. Yet, the specific choice for the control gains depends on the performance of the control. In order to prevent dramatic change of control inputs, excessively small s h , s ht , d a , d m0 are not recommended. However, smaller bounds usually generate smaller tracking errors. As a result, this is a tradeoff between robustness and accuracy.
Case 1 (Tracking sinusoidal trajectory): Here we select the original desired trajectory as x d = 1 sin (t), y d = 1 sin (t), z d = 2.5 sin (t), c d = p. A second-order filter is used to smooth the desired trajectory. In order to show the input constraint phenomenon, their feasible regions are set within a narrow scope as u t [ 8, 12 () , t x [ −10, 10 () , t y [ −10, 10 () , t z [ −15, 15 () . In order to validate the robustness of the controller, we apply zero mean random noises to the acceleration layer of rotation and translation, whose variances are 0.02 and 0.2, respectively. The obtained simulation results are shown in Figs. 3 and 4. From Figs. 3 and 4, it can be seen that the thrust force reaches its upper bound repeatedly, and the X, Y position errors increase at the same time. However, the system's tracking errors remain bounded, as indicated by the simulation results. The system meets no singularity, even though the yaw angle reaches 180 degrees.
Case 2 (Tracking linear trajectory): Desired trajectories of ramp functions are tested in this case. The controller parameters remain the same as those in Case 1. We choose the original desired trajectory as x d = 3t, y d = 3t, z d = 3t, c d = (p/2)t. Based on the simulation results shown in Figs. 5 and 6, we conclude that the vehicle is capable to track the ramp function with the proposed controller. It can be also seen from the simulation results that saturation occurs for the thrust force, yet its influences are tolerable due to the specific design of the control law. Moreover, although the noise has more significant influences on the torques than on the thrust force, the torques remain bounded and the tracking error of the yaw angle keeps relatively small.

Case 3 (Robustness test):
To examine the capability of handling large angular error, Case 2 is re-implemented with a new initial condition that the vehicle is upside down. After analysing the results shown in Fig. 7, it can be drawn that the angular controller remains effective even when the angular error is relatively large. The system gradually converges to the similar trajectory as Case 2.
It can be extracted from the simulation results that the fast oscillating instant thrust force is successfully approximated by a changing sinusoidal function with the form of F ins = k F v 2 m sin u p + k F v 2 m ,u p = v m , whose average value is modulated as F th . Using quaternion-based hybrid attitude tracking