Limit cycle oscillation suppression controller design and stability analysis of the periodically time-varying flapping flight dynamics in hover

The periodically time-varying forces make the equilibrium state of Beihawk, an X-shaped flapping-wing aircraft, to be a periodic limit cycle oscillation. However, traditional controllers based on averaging theory fail to suppress this oscillation and the derived stability result may be inaccurate. In this study, a period-based method is proposed to design the oscillation suppression controller, locate the corresponding cycle and analyze its stability. A periodically time-varying wing-tail interaction model is built and Discrete Fourier Transform is applied to adapt the model for controller design. The harmonics less than quintuple flapping frequency account for more than 96 percent of the total harmonics and are reserved to present a concise model. Based on this model, active disturbance rejection controller (ADRC) is designed and its Extended State Observer can observe the disturbance to suppress the oscillation. Poincaré map is introduced to convert the stability analysis of the cycle to a fixed point. A multiple shooting method is adopted to locate several points on the cycle and the map is obtained by calculating the submaps between the adjacent points with the Floquet theory. The located points are proved to be accurate compared with the numerical solved cycle and the stability analysis result of the cycle is verified by the dynamic evolution. Compared with the State Feedback Controller, the ADRC performs better in suppressing the limit cycle oscillation and eliminating the attitude control error. The oscillation suppression is meaningful in maintaining a stable flight and capturing high quality images.


Abstract
The periodically time-varying forces make the equilibrium state of Beihawk, an X-shaped flapping-wing aircraft, to be a periodic limit cycle oscillation. However, traditional controllers based on averaging theory fail to suppress this oscillation and the derived stability result may be inaccurate. In this study, a period-based method is proposed to design the oscillation suppression controller, locate the corresponding cycle and analyze its stability. A periodically time-varying wing-tail interaction model is built and Discrete Fourier Transform is applied to adapt the model for controller design. The harmonics less than quintuple flapping frequency account for more than 96 percent of the total harmonics and are reserved to present a concise model. Based on this model, active disturbance rejection controller (ADRC) is designed and its Extended State Observer can observe the disturbance to suppress the oscillation. Poincaré map is introduced to convert the stability analysis of the cycle to a fixed point. A multiple shooting method is adopted to locate several points on the cycle and the map is obtained by calculating the submaps between the adjacent points with the Floquet theory. The located points are proved to be accurate compared with the numerical solved cycle and the stability analysis result of the cycle is verified by the dynamic evolution. Compared with the State Feedback Controller, the ADRC performs better in suppressing the limit cycle oscillation and eliminating the attitude control error. The oscillation suppression is meaningful in maintaining a stable flight and capturing high quality images.
Keywords Wing Á Tail interaction Á Periodically time-varying dynamics Á Poincaré map Á Floquet theory Á Active disturbance rejection controller

Introduction
The superior maneuverability characteristics and high flight efficiency of bio-inspired flapping-wing aircraft have recently become a hot area of research [1,2]. However, the constant dynamic oscillation resulting from their inherent flapping propulsion mechanism makes it difficult to hover steadily like a rotorcraft. And this will be more serious in large aircraft because of the scale effect and thus existing hoverable flapping-wing aircraft are mostly insect-like [3][4][5]. So, this paper tackles the challenging point of the oscillation suppression control of a large flappingwing aircraft, Beihawk weighed 1.2 kg with 1.5 m wingspan [6], to enable a stable hover.
The hover maneuverability of flapping-wing aircraft can be attributed to their unsteady high-lift mechanisms which cannot be explained by traditional steady aerodynamic theories for fixed-wing aircraft due to their low Reynolds number regime and unsteady characteristics [7]. Great achievements have been made in the unsteady aerodynamic study of flapping flight in the past decades [7][8][9][10], among which, Dickinson et al. [7] summarized these mechanisms as delayed stall, rotational lift, added mass and wake capture. They investigated the aerodynamic force generation of a fruit fly by using a dynamically scaled robotic model and various instantaneous force measuring experiments were conducted. The obtained force data was fitted by mathematical empirical formulas to approximate these unsteady mechanisms. Owing to the high precision and intuitive analytical expression for the unsteady mechanisms, it has been widely used in the optimization of wing kinematics [11], dynamic modeling [12] and controller design [13].
The aerodynamic forces in flapping-wing aircraft are periodically time-varying and they are embedded into the dynamic equations of motion to constitute a periodically time-varying model. Such model represents a typical non-autonomous system, whose state equation is a function of time. Therefore, such controller design and stability analysis are challenging tasks [14]. To address this point, the assumption that the flapping frequency is much larger than the body's natural frequency is introduced and averaging theory is thus applied to simplify such a model into a timeaveraged (autonomous) one [15]. Averaging theory bridges the gap of controller design and stability analysis between the flapping-wing and fixed-wing aircraft and Deng et al. conducted a series of flapping flight control studies with this foundation. Linear quadratic Gaussian optimal controller [15], nonlinear controllers including a neural adaptive controller [16], a geometric controller [17] and a reinforcement learning control policy [18] were separately designed.
Basic attitude stability augmentation [16] and trajectory tracking [15,19] were firstly achieved, and then some extreme maneuvers [18] could be performed with the aid of these controllers. Wood et al. [20] adopted this theory as well and adaptive controllers were designed to cope with the uncertainties arisen from manufacturing imperfections in flapping flight with great success. Both Deng's and Wood's aircraft are micro vehicles with high flapping frequency and the body dynamics is mainly affected by the averaged forces. While for large ones, the frequency is strictly limited by the large inertial force of the wing and averaging operation may be invalid. Although great achievements have been made under the averaging theory, there still exist some debates about the valid range of this theory [21]. Taha et al. have made significant progress in proving the limitation of averaging theory and proposing more general theories for flapping wing aircraft [22]. Based on their aerodynamic-dynamic interaction model [23], they proved the need of high order averaging to understand the dynamics of flapping flight [24]. They also revealed a hidden vibrational stabilization mechanism using geometric control theory with the time-varying model [25,26]. Obviously, these interesting findings are beyond the research scope of the averaging theory. Besides, averaging theory converts the inherent dynamic oscillation into a point in flapping flight and controllers designed based on the averaged model behave badly in oscillation suppression.
Even more, the typical time-invariant controller stability analysis methods such as the Lyapunov method and Eigenvalue method based on the averaged model may fail to evaluate the practical periodic stability of flapping flight. This is because the equilibrium state of flapping flight is no longer a point but a periodic orbit. This orbit should not be averaged since only by including this periodic orbit can we explain the large attitude oscillations and the stall phenomena at some extreme instantaneous flapping positions in large flapping-wing aircraft. Many efforts have been done to consider this periodic property in flapping flight as fully as possible. Taylor et al. [27] firstly built a periodically time-varying flapping flight model of locust and redefined the stability issues as the asymptotic limit cycle stability in the phase space. Later researchers enriched this model by including the coupled electromechanical-aeroelastic dynamics [28], multi-body dynamics [29,30] and coupled unsteady aero-flight dynamics [31]. Based on the approximate linear periodically time-varying model of a flappingwing aircraft, Su et al. [32] performed the stability analysis using the transition matrix in Floquet theory. Nevertheless, this method is strictly confined by the linear model assumption. Dietl et al. [33] further improved this approach by introducing a shooting algorithm to search for the trimmed limit cycle of the open-loop dynamics and proved its instability. Yet they did not design controllers to stabilize the aircraft due to the constraints of the model complexity. Banazadeh et al. [34] studied the influence of parameters variation on the stability of the time-periodic flapping wing structure and indicated that control study was the work to be done. In summary, present periodic studies are not overall and systematic. They either concentrated on periodic dynamic modeling, or tried to evaluate the inherent periodic stability of the flapping creatures or the open-loop flapping-wing aircraft. And to date, there is no such research which exploits the related periodic theories to solve the practical control problem in flapping-wing aircraft to obtain a stable limit cycle and minimize its oscillation as far as possible.
The periodic issues are especially evident for the pitching attitude of the target aircraft in hover, so it is rare for such a large flapping-wing aircraft as Beihawk to hover. The main issues can be concluded as two aspects. First, in addition to the periodic thrust force for weight balance in hover, the X-shaped flapping wing also generates some periodic disturbing torque which would exacerbate the attitude oscillation. Second, the flapping wing generates thrust force through accelerating the surrounding air and a strong induced flow is generated. This flow is the source of the inflow over the tail in hover stage and possesses the same periodic time-varying principle as the thrust force, which would reduce the tail control effect in oscillation suppression. The previous work [6] has modeled this flow and the disturbance, and designed simple PD controller to validate this model. The designed controller could maintain the aircraft in hover and the pitching attitude oscillation was not further analyzed. Indeed, these oscillatory behaviors inevitably degrade the quality of images obtained and cause frequent malfunctions in mounted onboard microelectronics that are vulnerable to external vibrations. Therefore, it is necessary to build a both accurate enough and easily analyzed periodic model, design periodic controllers to suppress oscillation based on this model and prove the stability of the controllers from the periodic limit cycle viewpoint.
The objective of this paper is to investigate an easily implemented periodically time-varying model and design an efficient limit cycle suppression controller with the wing-tail interaction model. In particular, the stability of the designed controller is analyzed by a period-based method to consider the limit cycle characteristic in flapping flight. A pitching attitude model of Beihawk including unsteady effects is firstly presented. Specially, the flapping-wing induced flow over the tail is considered to give a wing-tail interaction model. Discrete fourier transform (DFT) is performed on the aerodynamic forces and the significant harmonic contents are incorporated into the dynamic model to present a concise model. The model-free state feedback controller (SFC) and the model-based active disturbance rejection controller (ADRC) are separately designed to suppress the oscillation. The stability of the controller is proved by the stability of the corresponding equilibrium limit cycle. A multiple shooting method based on Levenberg-Marquardt Algorithm (LMA) is utilized to locate the cycle and its stability is described with the characteristic multipliers by introducing the Poincaré map and Floquet theory.
The rest of this paper is organized as follows. Section 2 presents the periodically time-varying wing-tail interaction attitude model of Beihawk. Section 3 describes the DFT of the model and the attitude controller design. Section 4 reports the period-based stability analysis method of the designed controllers. Section 5 presents the limit cycle oscillation suppression effects. Section 6 summarizes the conclusions and proposes some potential applications.

Attitude dynamic model of Beihawk
The wing aerodynamics, tail aerodynamics and body dynamics form the model of Beihawk. Three significant high-lift mechanisms which are delayed stall, rotational lift and added mass are considered in the wing model. The rotation of the wing manipulates not only the magnitude of the aerodynamic forces but also the orientation [35]. Taha et al. [26] even pointed a new pitch damping effect due to the rotational contributions. However, the averaged rotational lift may equal to zero in some cases [36], which also indicates the drawback of the averaging theory and the necessity of a periodically time-varying model. The wing generates thrust force through accelerating the momentum of the surrounding air and thus momentum theory [37] is adopted to solve the induced flow velocity based on the thrust force. This flow is the main incoming flow over the tail in hover and a periodic tail model including the wing induced flow is thus established. The wing model and tail model are incorporated into the dynamic model to give the whole periodically time-varying model.

Wing aerodynamics including unsteady effects
The blade element theory with the aerodynamic coefficients identified from experiment is adopted to approximate the unsteady effects. The flapping angle / and twisting angle b are given by In a quasi-steady model, the force equations derived from thin aerofoils translating with constant velocity and angle of attack also apply to time-varying flapping wings [12]. Therefore, the translational aerodynamic forces exerted on per blade element aerofoil are given by where dF tr;N and dF tr;T separately represent the normal force and tangential force of the wing. Dickinson et al. [38] solved the force coefficients due to delayed stall as follows, Fung et al. [39] derived the theoretical force exerted on per blade element aerofoil due to rotational lift as where C rot ¼ 2pð0:75 Àx 0 Þ is the rotational coefficient,x 0 is the dimensionless distance of the longitudinal rotation axis from the leading edge, x is the instantaneous angular velocity of the wing, approximating to _ aðtÞ. It is a pure pressure force and therefore acts perpendicularly to the wing profile [12,35].
In [40], an additional normal force contributed by the added mass effect is given by where _ v 2 is the normal acceleration at the mid chord due to the wing motion.
The total normal force and tangential force exerted on per blade element aerofoil are thus given by, dF N ðt; rÞ ¼ dF tr;N ðt; rÞ þ dF rot;N ðt; rÞ þ dF add;N ðt; rÞ; ð7Þ These forces are rotated into body axis by The up-wing and down-wing are different due to the design of the dihedral angle / 0 (see Fig. 1). Therefore, the flapping and pitching motions of the upwing and down-wing can be listed as follows: The whole wings' instantaneous normal and thrust forces are obtained by integration as follows, The detailed introduction to the layout and description of the parameters of Beihawk can refer to the previous work [6].

Tail aerodynamics considering wing-tail interaction
Momentum theory [37], modeling the circular swept area of the flapping wing as an actuator disk that accelerates the surrounding air, is used to calculate the flapping wing induced flow velocity over the tail. The momentum change of the induced flow generates the thrust force and its velocity at the wing root is given by where T(t) is the thrust force and S d is the stroke plane area. This flow changes along the body and the velocity at x (see Fig. 1 for axis definition) position is given by [41] where R is the diameter of the actuator disk, i.e., b/2. This formula originates from the calculation of the propeller-induced velocity at a distance from the propeller disk [41] and has been proved to be efficient in the calculation of the distribution of the flapping wing induced flow [42]. Beihawk's tail is installed approximately 0.5 m from the wing root, and the induced flow velocity at the tail, i.e., V iT , can then be obtained.
Bernoulli's equation is used to calculate the aerodynamic force acting on the tail as follows, where d is the tail pitch angle. With thin flat-plateairfoil assumption, the lift-curve slope of a finite span three dimensional wing is given by [43] C a l ¼ where AR is the aspect ratio, and s is 0.04 for the designed wing [44].

Periodically time-varying model of Beihawk
The pitch torque generated by the wing can be calculated based on the center of pressure (CP) of the wing, A W , Considering the fact that the flapping amplitude is only 15 deg, the position variation of CP is not significant. Therefore, the CP position of each flapping wing is assumed to be located at the average flap plane of each flapping wing, and at 1/4c position along the chord direction. The overall CP of a pair of flapping wings is then decided at the middle of the two CP point of each flapping wing.
The tail torque is calculated based on the CP of the tail, i.e., A T , where x 0 T is the distance between A T and G along the x axis.
The longitudinal damping of Beihawk is much smaller than the lateral and thus hovering stability is mainly affected by the pitching attitude dynamics which is given by The controlled variable is tail pitch angle and Eq. (18) can be further summarized as follows according to Eq. (14) _ hðtÞ ¼ qðtÞ It is obvious that Eq. (19) is periodically timevarying under the periodically forced terms NðtÞ and TðtÞ. The tail torque is influenced by the thrust force generated by the wings, and thus wing-tail interaction is fully considered. Beihawk relies on the tail to suppress the attitude oscillation by counteracting the disturbing torque generated by flapping wings. However, the tail torque possesses the same periodic timevariation principle as the thrust force. This variation characteristic would reduce the tail control effect if we still adopt traditional controllers which are designed based on constant tail control torque.

Attitude controller design
The complex aerodynamic formulas in flapping flight make it impossible to give a concise aerodynamic model which can be easily embedded into the dynamic model to constitute an appropriate model for controller design. Therefore, present researchers either directly control the aerodynamic forces without considering the aerodynamic model, or simplify the complex model by employing direct averaging theory. Apparently, both assumptions neglect the most important time-varying features of unsteady aerodynamics.
In this section, we adapt the model built in Sect. 2 to controller design by taking discrete fourier transform (DFT) to the instantaneous aerodynamic forces. The transformed force is composed of infinite discrete fourier series: the fundamental harmonic frequency is the same as the flapping frequency, the direct component corresponds to the averaged value, and the harmonics far above the body's natural frequency will not cause significant changes on the body dynamics and thus can be neglected. The above operation makes the original model convenient for controller design and stability analysis, and reserves the periodically time-varying features as well.
Based on the transformed model, two controllers are designed. The model-free state feedback controller (SFC) is designed to primarily achieve the stability augmentation of the hover maneuver. Advanced model-based controller, active disturbance rejection controller (ADRC) is further designed to suppress the oscillation. The extended state observer (ESO) in ADRC can make full use of the built periodic wingtail interaction model to observe not only the model state but also the disturbance as the extended state. The ESO has been proved to be efficient in disturbance observation and the observed disturbance can then be counteracted in the control signal design.

Main fourier series representation of the model
The normal force and thrust force generated by the flapping motion are functions of the periodic wing kinematics and can thus be presented by the fourier terms, i.e., ½a n cosð2pfntÞþb n sinð2pfntÞ where n is an integer and f is the flapping frequency. Any non-harmonic spectral content can be treated as noise. Set the hover flapping frequency as 8 Hz and we can obtain the corresponding aerodynamic forces over one whole flapping period. One thousand uniformly distributed data points in a period are selected for DFT and the results are presented below, The amplitude-frequency characteristic diagram (Fig. 2a) indicates that the harmonics below 4f (include) are the main harmonic content of the thrust force. Besides, the harmonic above 4f is much higher than the body's natural frequency and the body dynamics is mainly influenced by the averaged value. The contribution of these harmonics can be overlooked, since the averaged value equals zero. Therefore, the main harmonics are chosen to approximate the original thrust force and the result is shown in Fig. 2b. It is nearly the same as the original force.
The same operation is carried out for the normal force and the result is shown in Fig. 3. We also select harmonics below 4f (include) as the main harmonic content of the normal force. The transformed thrust force and normal force are represented as ½a Tn cosð2pfntÞþb Tn sinð2pfntÞ; ð21Þ The percentage of each harmonic content in the thrust force and normal force is calculated separately and the results are listed in Table 1. It can be easily seen that the harmonics below 4f (include) account for more than 98% of the whole harmonics in the thrust force and 96% in the normal force, which also proves the rationality of our approximation. This approximation method has also been validated under other frequencies and here only the results of the hovering frequency are presented.
The main fourier series representation of Eq. (19) is written as Suppose Equation (24) is also a non-autonomous system, or more specifically a periodically time-varying system and the built wing-tail interaction model is reflected in bðtÞ, which is a time-varying term along with the periodic flapping motion.

Limit cycle oscillation suppression controller design
In this section, two controllers are separately designed to suppress the attitude oscillation of Beihawk, the basic SFC and the advanced ADRC. The SFC is designed as follows, where x 1d and x 2d are separately the expected value of x 1 and x 2 , k p and k d are the corresponding feedback gains. SFC is simply designed without including model information. It can only achieve stability augmentation and cannot well reject the disturbing torque generated by the wings. The tail control effect is noticeably reduced under this traditional SFC because of the time-varying induced flow. Therefore, it is urgent to design model-based controllers to take full advantage of the built model. ADRC was firstly proposed by Han [45] and has been widely used in the industrial control [46] and quadrotor control [47] due to its excellent disturbance rejection capacity. It is employed here with the following superiorities. First, it can take full advantage of the wing-tail interaction model built in Sect. 2.2 to achieve more accurate disturbance observation and better disturbance rejection. Second, it can observe the pitching attitude velocity and thus the attitude velocity sensors can be replaced, which would reduce the weight of the aircraft and prolong the endurance. The ADRC is implemented as follows.
In fact, in the context of feedback control, dðtÞ in Eq. (24) is something to be overcome by the control signal. Treating dðtÞ as an additional state variable, x 3 ¼ dðtÞ, and let _ dðtÞ ¼ GðtÞ with GðtÞ unknown, the original plant Eq. (24) is now described as which has been proved to be always observable by its proposer [45]. Then, the periodic ESO can be constructed as follows, where the system output x 1 and the control signal uðtÞ are the inputs to ESO, and the output of the ESO,x 1 ,x 2 andx 3 are respectively the observation of x 1 , x 2 and dðtÞ. l 1 , l 2 and l 3 are the gains of the observer. There is no doubt that better observation results can be obtained with the help of the built tail model bðtÞ.
Note that the important information dðtÞ, provided by the ESO, allows the control law uðtÞ ¼ ðu 0 À dðtÞÞ=bðtÞ to reduce the plant in Eq. (24) to a cascade integral form of The sufficient condition for this reduction is bðtÞ 6 ¼ 0. Therefore, the ADRC can utilize the built wing-tail interaction model to compensate the plant and the compensated model can be easily controlled by making u 0 a function of the tracking error and its derivative, i.e., a SFC. Based on the above description, the periodic control signal of the ADRC in our paper is designed as 4 Stability analysis of the designed controller The designed controllers in Sect. 3 are also periodically time-varying since it is based on the periodically time-varying model. Considering this periodically time-varying feature, there is no equilibrium state where all states remain constant as the equilibrium point in fixed-wing aircraft. Therefore, general stability analysis methods for fixed-wing aircraft cannot be directly utilized to evaluate the stability of flappingwing aircraft. The equilibrium state of the flapping wing aircraft is a periodic orbit and is defined as a limit cycle. The stability of the system can be proved by proving the stability of the limit cycle. The defined limit cycle is stable if any trajectory starting near the limit cycle converges on it. If any trajectory starting near the limit cycle diverges from it, the limit cycle is unstable.
Our method converts the stability of a limit cycle in the flapping wing aircraft to the stability of a fixed point in a specific hyperplane by introducing the Poincaré map. The stability of the fixed point is determined by the characteristic multipliers of the Poincaré map's Jacobian matrix. To solve the characteristic multipliers of the original system, a variation equation of the cycle is introduced and the system is transformed into a linear periodically time-varying system. Floquet theory applied for linear periodically time-varying system is utilized to solve the monodromy matrix and this matrix can be regarded as a Poincaré map. Floquet multipliers of this matrix are thus proved to be coincident with the multipliers of the aforementioned map. Finally, the key theorem on the stability of the periodically time-varying flapping flight model is given.
The exact location of the limit cycle is prerequisite for the monodromy matrix. A multiple shooting method is introduced to locate the limit cycle by shooting several points on it. Levenberg-Marquardt Algorithm (LMA) is adopted to iterate these solutions until these points converge on the circle. Two adjacent points can be related by the transition matrix and all these transition matrixes are combined to form the monodromy matrix. The LMA can not only reduce the reliance on the initial value but also obtain the corresponding monodromy matrix simultaneously in the last iteration. The Floquet multipliers can thus be solved and used to analyze the stability.

Location and stability analysis of the limit cycle
The model of Beihawk, Eq. (24), can be written as a general periodically time-varying model, where f ðt þ T; xÞ ¼ f ðt; xÞ and T is the flapping period, x represents the states of Beihawk. The equilibrium state of Beihawk is a n dimensional limit cycle in the phase plane. Assume the limit cycle is L 0 and a point x 0 2 L 0 . For convenience, we transform the model of Beihawk as follows, with coordinates ða; xÞ 2 S 1 Â R n (a ¼ tðmodÞT). After the transformation, Eq. (31) now has n ? 1 dimensional periodic solutions ða; xÞ.
Definition 1 For Eq. (31), we define the following n dimensional hyperplane along the flapping cycle a, x 2 R n is the coordinate on R. xða; x 0 Þ is a solution to Eq. (32) with initial value x 0 on the interval a 2 ½0; T. Based on the defined hyperplane for the flapping wing model, the corresponding Poincaré map is introduced (shown in ?tic=?>Fig. 4), It is obvious that this map describes the states evolution of the aircraft after a flapping period. If the point x starting close enough to x 0 can finally converge to x 0 after several Poincaré maps, the fixed point x 0 is stable and the limit cycle L 0 is stable as well. After the introduction of the Poincaré map, the stability of the limit cycle L 0 is now rephrased as the stability of the fixed point x 0 .

Lemma 1 ([48])
The Jacobian matrix of the Poin-care´map P taken at the point x 0 is defined as, Its eigenvalues are called characteristic multipliers. If the characteristic multipliers have complex modulus less than one, the fixed point is stable.

Remark 1 ([48])
The characteristic multipliers are independent of the hyperplane choice for the Poincaré map and the fixed point, because they are a property of the limit cycle. The next problem to be addressed is how to solve the characteristic multipliers of the flapping wing model.
A linear periodically time-varying system with period T 0 is introduced The eigenvalues k 1 ; k 2 ; :::; k n of UðT 0 ; 0Þ are called Floquet multipliers.
Base on the above, the key theorem on determining the stability of the periodically time-varying flapping wing aircraft model is presented.
Theorem 1 For the periodically time-varying flapping wing model of Beihawk, Eq. (30), with a limit cycle L 0 , if all Floquet multipliers corresponding to the cycle have complex modulus less than one, the limit cycle is stable.
Proof Let x ¼ uðtÞ denotes the periodic solution to the system Eq. (30). Represent a solution to Eq. (30) in the form where n is a deviation from the periodic solution. Then, the system can be linearized around the cycle, where GðtÞ ¼ of ox x¼uðtÞ . The original periodically time-varying system Eq. (30) has been reduced to a linear periodically time-varying system Eq. (39) locally around the cycle.
From the definition of the monodromy matrix in Definition 2, we can find that it describes the map of a flapping period, and thus is an option of the Poincaré map specified in Definition 1 for the flapping wing aircraft. Therefore, the Floquet multipliers of UðT; 0Þ coincide with the multipliers of the Poincaré map P. According to Lemma. 1, if the Floquet multipliers have complex modulus less than one, the fixed point is stable and the corresponding limit cycle is stable. Hence, Theorem 1 holds.h The monodromy matrix of the cycle L 0 can be obtained only after exactly locating the limit cycle. The location of the limit cycle is apparently more complicated than that of the equilibrium point, which also hinders the relevant research.
In this paper, a multiple shooting method which was developed in [50] is introduced to locate the cycle. After obtaining the exact points position, the state transition matrix between the adjacent points is calculated and all these transition matrices can be combined to form the corresponding monodromy matrix. It is evident that this division can significantly Fig. 4 The defined Poincaré map associated with the cycle for Beihawk improve the accuracy of the solution by depicting a more detailed cycle.
The multiple shooting method is implemented as follows. The circle is divided into m non-empty subintervals 0 ¼ s 0 \s 1 \:::\s m ¼ T and the corresponding coordinate position of the points on the circle are respectively x 0 ; x 1 ; :::; x mÀ1 ; x m ; x j 2 R n (shown in Fig. 5). Note that the point x j is composed of several states ½x j 1 ; x j 2 ; . . .; x j n T . The time intervals are defined as D j ¼ s jþ1 À s j ; j ¼ 0; 1; :::; m À 1 and fðx j Þ s jþ1 s j denotes the integral results from s j to s jþ1 with the initial point value x j .
The Poincaré map is also divided into several submaps correspondingly based on its semigroup property [51] P ¼ P mÀ1 Á ::: where P j : R j ! R jþ1 ; j ¼ 0; 1; 2; :::; m À 1, and R j is the hyperplane on the point x j . The fundamental state transition matrix also shares the semigroup property, UðT; 0Þ ¼ UðT; s mÀ1 ÞUðs mÀ1 ; s mÀ2 Þ:::Uðs 1 ; 0Þ: The multiple shooting problem in this study can be formulated as finding the pointsx ¼ ½x 0 ;x 1 ; Á Á Á ;x mÀ1 ;x m ;x 2 R ðmþ1ÞÁn that can better fit the points x ¼ ½x 0 ; x 1 ; Á Á Á ; x mÀ1 ; x m on the cycle. The concrete computation process is implemented as follows.
The error vector is firstly defined as ::::::: For convenience, Eq. (42) is rearranged by defining the left side as the target error function g ¼ ½ g 0 g 1 Á Á Á g m T ; g 2 R ðmþ1ÞÁn and the right side as the error vector e ¼ ½ e 0 e 1 Á Á Á e m T ; e 2 R ðmþ1ÞÁn , The target is now transformed into minimizing the sum of the squares of the errors.
The LMA solves the optimization problem by iteration: initiated at the starting pointx 0 , the method produces a series of solutionsx 1 ;x 2 ; . . . until the error condition is satisfied. The required Jacobian matrix in the LMA is defined as where ofðx j Þ s jþ1 s j ox j 2 R nÂn ; j ¼ 0; 1; . . .; m À 1 and is computed by a numerical method where e is the differentiation step andê is the unit vector.
The advantage of the LMA is its state updating policy and the state augmentation dx is calculated as follows, where k is a non-negative and adaptive damping parameter. Detailed calculation process is presented in the Appendix. When the trajectory eventually converges on the limit cycle, the norm of the error vector is approximately zero. Thus, The monodromy matrix of the whole cycle can be obtained based on Eq. (41), ox mÀ2 Á ::: All of the elements in Eq. (48) can be found in Eq. (44), which presents that the monodromy matrix does not need to be solved separately but is obtained naturally after the location of the limit cycle. The eigenvalues of the monodromy matrix can then be solved and applied to the stability analysis of the system.

Stability analysis results of the designed controllers
The stability of the designed controller in Sect. 3.2 is analyzed by the limit cycle location and stability analysis method proposed in Sect. 4.1.
The SFC is firstly analyzed and the located limit cycle is shown in Fig. 6. The triangles in red are points located by the proposed method and the blue curve describes the practical limit cycle by the numerical integration of the practical system dynamics. It is evident that all points are on the cycle and thus the method can give a precise location result. There is no doubt that an accurate limit cycle location result contributes to the solution of the monodromy matrix and the derived multipliers.
The corresponding calculated Floquet multipliers (circles in red) are plotted in Fig. 7. These multipliers are inside the unit cycle, and thus the designed SFC is stable. A grid of constant damping factors from zero to one in steps of 0.2 and natural frequencies from zero to p in step of p/5 is also plotted over the current axis. It can be seen that the damping factor of the system is significantly improved compared with the open loop and the aircraft can maintain a more stable hover flight. Here we mainly focus on whether the controller is stable or not. Meanwhile, the dynamic characteristics can be further improved by adjusting related parameters of the aircraft and the controller with the grid graph of the multipliers in Fig. 7.
To further verify the stability analysis result by the proposed theorems, we separately choose a random point inside the limit cycle and a random point outside the limit cycle as the initial states and analyze their practical dynamic evolution. Their convergent trajectories are shown in Fig. 8 and the located limit cycle in Fig. 6 is replotted as well. Obviously both the inside one and the outside one converge to the limit cycle. Therefore, the located limit cycle is stable and the proposed method is feasible.
In the stability analysis of the ADRC, the stability issues of the ESO should also be included. The ADRC including the ESO can be presented as follows, Fig. 6 The located points and the practical limit cycle under the SFC (Triangles in red are the located points and blue curve describes the practical limit cycle) Fig. 7 The position of the Floquet multipliers under the SFC in z-plane This system has five states including the practical states and the observed states. In this method, the observer and the controller can be easily listed together and the stability of the whole system can be obtained. This also presents the advantage of our method: the stability issues of any controllers which can be expressed as a periodically time-varying differential system like Eq. (30) can be analyzed via this method.
These overall five states in Eq. (49) cannot be plotted in one coordinate system. Thus, the practical states and the observed states are separately plotted in two coordinate systems shown in Figs. 9 and 10. The results also indicate that this method can give an accurate location of the limit cycle.
The corresponding Floquet multipliers (circles in red) are also plotted in Fig. 11. The five multipliers inside the unit cycle also prove the stability of the ADRC.
The stability analysis result of the ADRC is also tested and the corresponding dynamic evolutions are shown in Figs. 12 and 13. The trajectories on both sides of the limit cycle converge on the limit cycle, which also confirms the stability analysis method proposed in this paper.
In this section, the proposed method is utilized to analyze the stability of the periodically time-varying flapping-wing aircraft under the SFC and the ADRC without the prior assumption of the averaging theory. Fig. 8 The convergent trajectories of the point inside the limit cycle (in blue) and outside the limit cycle (in red), and the located limit cycle (dash line in black) under the SFC Fig. 9 The located points and the limit cycle of the practical states in ADRC (Triangles in red are the located points and blue curve describes the practical limit cycle) Fig. 10 The located points and the limit cycle of the observed states in ADRC (Triangles in red are the located points and blue curve describes the practical limit cycle) The analysis results are validated by simulation and proved to be accurate. Thus, the controller design and stability analysis in flapping-wing aircraft is no longer limited by the range of flapping frequency. This method can also give a more detailed description of the time-varying features in flapping-wing aircraft, like the interesting limit cycle oscillation, which cannot be observed in the averaged system.

Flight control results
After ensuring the stability of the designed controllers, we separately use them to control Beihawk in hover.
For testing, the following standard 4-state longitudinal flight dynamics are considered, where u is the velocity along the body, w is the velocity normal to the body, T is the thrust force along the body, N is the normal force normal to the body, M is the pitching moment including the pitching moment generated by the wing, M W , and the pitching moment generated by the tail, M T . The calculation formulas of the above forces and moments have been introduced in Sect. 2. Besides, the contribution of the body motion to the aerodynamic loads which is proposed by Taha et al. [23] is also added to our aerodynamic model to capture the essential dynamical features of flapping flight.
The designed attitude controller is tested on this model. The velocity controllers adopt the traditional PID controller and the control parameters are tuned to obtain a satisfactory result. The final control results are shown as follows.

State feedback controller
The parameters of the SFC are set to be k p ¼25 and k d ¼ 10 to obtain a satisfactory control effect. The desired velocities, pitching attitude and attitude velocity are set to be zero to stay hover status and the control effects are shown in Figs. 14, 15. Fig. 11 The position of the Floquet multipliers under the ADRC in z-plane Fig. 12 The convergent trajectories of the point inside the limit cycle (in blue) and outside the limit cycle (in red), and the located limit cycle (dash line in black) of the practical states in ADRC Fig. 13 The convergent trajectories of the point inside the limit cycle (in blue) and outside the limit cycle (in red), and the located limit cycle (in black) of the observed states in ADRC It can be seen that Beihawk experiences large velocity oscillations along the body direction in hover, which is due to the significant periodic change of the thrust force. While the oscillation of the normal velocity, w, is small. This can be attributed to the design of the X-shaped flapping wings and the normal forces of the symmetrical flapping wings almost counteract. The results are consistent with the thrust force and normal force calculation results shown in Figs. 2 and 3. With the SFC of the tail, the pitching attitude of the aircraft can achieve the hover stability augmentation (shown in Fig. 15). This indicates that the tail control mechanism is important for a large flapping wing aircraft to achieve stable hover. However, large attitude oscillations still exist in both the pitching attitude and pitching attitude velocity. One reason for the oscillation is that the tail control effect is significantly reduced under the periodically timevarying induced flow over the tail if we still adopt traditional SFC. The other reason is due to the disturbing torque generated by the flapping wings. The aerodynamic-dynamic interaction effect included in the testing model would also influence the tail control effect and the disturbing torque. Therefore, it is necessary for us to design advanced attitude controller by including this periodic induced flow over the tail and compensating the disturbing torque.

Active disturbance rejection controller
The parameters k 1 and k 2 of the ADRC are designed the same as the SFC to highlight the superiority of the ESO in the ADRC. The parameters of the ESO are designed according to the bandwidth parameterization results in [52]: where x 0 is the observer bandwidth. Here x 0 is designed to be 1000 considering the maximum 10 Hz flapping frequency. The control effects and the observed results are shown in Figs. 16, 17, 18. From Fig. 16, we can observe that the velocity oscillation along the body axis under the ADRC is comparable to that of the SFC. Since this oscillation is mainly due to the time-varying feature of the thrust force and the velocity controller still adopts the traditional PID controller. While the normal velocity oscillation under the ADRC is smaller than that of the SFC and this can be attributed to the better pitching attitude oscillation suppression effect of the ADRC which can be seen in Fig. 17. This is no doubt that better attitude oscillation suppression effect can contribute to the velocity oscillation suppression and the attitude stabilization is of great significance to the  Fig. 17, we can clearly observe that the attitude oscillation is significantly suppressed and the attitude of Beihawk tends to be more stable. Both Figs. 17 and 18 indicate that the ESO in ADRC can not only accurately observe the system state but also give a satisfactory observation of  the disturbance although the varying frequency of the generated disturbance is considerably high, the same as the flapping frequency. The high accuracy of observation can be owing to the built accurate wingtail interaction model. By compensating the observed disturbance, the ADRC can effectively suppress the attitude oscillation. The compensation of the disturbance also relies on the built wing-tail interaction model. That is to say, an accurate tail model is necessary for the implementation of ADRC and better oscillation suppression effect can be obtained with a more accurate tail model. Finally, as we all know, the limit cycle oscillation is the inherent feature of flapping-wing aircraft. It cannot be entirely counteracted and can only be suppressed. The ADRC can better suppress the oscillation amplitude than the SFC and give a satisfying attitude control result. The Mean Value (MV) and Root Mean Square (RMS) of the error under both controllers are separately calculated as listed in Table 2 and the results quantitively confirm the above conclusions.
Both the SFC and the ADRC can achieve the stability augmentation of the flapping-wing aircraft, Beihawk, in hover. This furtherly confirms the validity of the proposed periodically time-varying stability analysis method. The ADRC performs better in attitude limit cycle oscillation suppression than the SFC. Because it can make full use of the built wingtail interaction model to observe the disturbance with the aid of the ESO and compensate the observed disturbance in the control signal design process. This also emphasizes the importance of an accurate tail model in the attitude controller design.

Conclusion
The widely used averaging theory fails to reveal the period limit cycle properties in flapping flight and may lead to erroneous results on the stability analysis. In this study, the critical periodically time-varying feature of the flapping-wing aircraft Beihawk is considered in dynamic modeling, controller design and stability analysis. This consideration facilitates the unique limit cycle oscillation study in flapping flight. To present an easily implemented periodically timevarying model, DFT is applied to the original model and the transformed model is found to be of 96 percent fidelity by reserving the harmonics less than quintuple flapping frequency. Based on this model, two controllers are separately designed to suppress the inherent oscillation. The stability of the controller is analyzed in terms of the unique equilibrium limit cycle which is apparently more complex than the equilibrium point. Poincaré map is introduced to simplify this cycle to an equivalent fixed point. A multiple shooting method is employed to locate several points on the cycle and calculate the submaps between the adjacent points. All these submaps are combined to form the Poincaré map and this division markedly improves the accuracy of the solved whole map. The characteristic multipliers associated with the Poincaré map can give an estimation of the fixed point's stability. The located points are discovered to be exactly distributed on the practical limit cycle solved by the numerical integration of the dynamics. The stability analysis results characterized by the multipliers coincide with the practical dynamic evolution represented by the convergent trajectories inside and outside the cycle. The proposed method can accurately describe the periodic feature and the stability of flapping flight and relevant conclusions are no longer constrained by the range of flapping frequency.
The designed controllers are used to control Beihawk after ensuring their stability. Both controllers can achieve stability augmentation while the ADRC performs better in suppressing the oscillation. This superiority of the ADRC can be attributed to the built wing-tail interaction model which is fully utilized to aid disturbance observation in the ESO and disturbance compensation in the control signal design. This also indicates the importance of tail control mechanism in the attitude stabilization and the necessity of an accurate tail model. The limit cycle oscillation suppression presents powerful research potential after realizing the stability augmentation. It enables Beihawk to achieve a more stable flight and to capture higher resolution images while carrying a camera. The proposed period-based stability analysis method is of high flexibility and it can be easily extended to analyze the stability of other controllers in flapping flight. Moreover, it has potential applications in analyzing the influence of parameters variation on the damping and frequency features of flight dynamics which is described by the position of the characteristic multipliers in the grid graph.