## Abstract

With the goal of operating a biologically inspired robot autonomously outside of laboratory conditions, in this paper, we simulated wind disturbances in a laboratory setting and investigated the effects of gusts on the flight dynamics of a millimetre-scale flapping-wing robot. Simplified models describing the disturbance effects on the robot's dynamics are proposed, together with two disturbance rejection schemes capable of estimating and compensating for the disturbances. The proposed methods are experimentally verified. The results show that these strategies reduced the root-mean-square position errors by more than 50% when the robot was subject to 80 cm s^{−1} horizontal wind. The analysis of flight data suggests that modulation of wing kinematics to stabilize the flight in the presence of wind gusts may indirectly contribute an additional stabilizing effect, reducing the time-averaged aerodynamic drag experienced by the robot. A benchtop experiment was performed to provide further support for this observed phenomenon.

## 1. Introduction

Foreseeable applications of small flying robots in our everyday lives, ranging from city courier services to automated aerial construction [1], have driven research advances in the area of *micro aerial vehicles* (MAVs) as seen in numerous examples [2,3]. Among these, multi-rotor vehicles have gained considerable popularity, thanks to simple mechanical designs and well-understood dynamic properties.

Biologically inspired MAVs are another platform of active research [4–8]. Flapping-wing devices are of particular interest owing to their manoeuvrability, as exemplified by their natural counterparts. Over the past decades, there have been several examples of flight-worthy flapping-wing robots, ranging from metre scale [4] to centimetre [5,6] and millimetre scale [9,10]. Inspired by flying insects, in [10,11], researchers have developed and demonstrated stable unconstrained flight of a fly-sized robot (shown in figure 1). Those demonstrations reflect the culmination of research in mesoscale actuation and manufacturing technology [11,12], as well as in understanding the control, stability and the aerodynamics of flapping-wing flight [8,13–15].

Despite having achieved stable flight, the flapping-wing robot in figure 1 is still tethered for power, sensing and control. Limited by the payload capacity (the robot can carry no more than 30–40 mg in addition to its current weight of 80 mg), sensing and control computation are executed offboard.

Thus far, we have demonstrated hovering flight [10], basic flight manoeuvres [15], perching [16] and realization of an acrobatic pre-planned flight trajectory such as landing on a vertical surface [17]. In parallel, research activities are being carried out to enable this vehicle to operate autonomously outside of laboratory conditions. Multiple sensors with suitable power and mass have been explored, ranging from a microelectromechanical inertial measurement unit [18] to a small footprint, onboard vision sensor inspired by insect ocelli [19]. Lightweight and efficient power electronics have also been designed [20]. In preparation for outdoor operations, another important requirement is the ability to overcome the effects of external disturbances such as wind gusts. An ability to robustly fly and navigate in the presence of wind disturbance would potentially allow the flapping-wing robot to autonomously operate in diverse natural environments, which necessitates better understanding of complex interactions between the robot locomotion and the real world. To facilitate the process, controlled experiments in the laboratory are essential as seen in various examples involving the study of robot locomotion in complex environments [21,22]. Simplified principles can then be extracted from emerging physical phenomena, allowing researchers to understand and develop robots with enhanced performance across diverse environments.

In this work, we seek to stabilize and control the position of an insect-scale flapping-wing robot in the presence of wind gusts, using only position and orientation feedback provided by an external motion capture system. To realize control in the presence of gusts, we first seek to acquire a better understanding of the impact of gusts on flapping-wing flight. To compensate, the robot is required to estimate the force and torque effects caused by the wind disturbances and apply corrective commands to stabilize its attitude and maintain the desired position without the *a priori* information on magnitude and direction of the wind.

### 1.1. Previous work on disturbance rejection in micro aerial vehicles and flapping flight

To date, there has been a number of efforts made to control aerial vehicles in windy environments. Examples include attempts that involve aircraft harvesting energy from atmospheric phenomena such as thermals and updrafts, and efforts to stabilize and control different types of aerial vehicles in wind gusts [23–26]. In this context, the airspeed of fixed-wing vehicles is typically much greater than the wind speed, and the stability and control of the vehicles are not compromised. In contrast, smaller gliders with relatively low nominal airspeed are likely to suffer from atmospheric turbulence. This requires the flight controller to take into account the effects of wind and compensate to ensure stability. For rotary-wing vehicles, early research on this topic was based on the rudimentary understanding of the effects of airspeed on thrust and drag generated by the propellers. More recently, attempts have been made to estimate and model how the direction and magnitude of wind affects thrust generation and vehicle drag [27]. On small rotary platforms, observer-based methods have been employed to control and actively reject disturbances [28,29].

The effects of wind disturbances on flapping flight are different from fixed-wing and rotary-wing flight in nature owing to unsteady aerodynamics. Literature related to wind gusts and flapping flight could be categorized into either the study of flapping-wing aerodynamics in the presence of gusts or turbulence [30–34] or the study of stability in forward and lateral flight [35–40]. Research focusing on aerodynamics typically employs computational fluid dynamics (CFD) or flow visualization to quantify instantaneous (as opposed to stroke averaged) flow around the wings and deduce corresponding aerodynamic properties such as lift and drag coefficients. For instance, Jones & Yamaleev [32] inspected turbulent flow from unsteady compressible Reynolds-averaged Navier–Stokes equations on a centimetre-scale wing and revealed that frontal gust has a very strong effect on the instantaneous thrust force. In [33], Fisher experimented on a flapping robotic wing (at *Re* = 43 000) at the geometric angle of attack in the range of 0–12° using a pressure measurement system and found that the average pressure distributions were significantly affected by turbulence. Based on the results from a revolving and translating dynamically scaled wing, Dickson & Dickinson [30] proposed a modified quasi-steady model for accurately predicting instantaneous force produced during forward flight. While these studies provide valuable insights into the aerodynamics of flapping flight, they do not directly attempt to quantitatively model the effects of airspeed on average lift and drag as a function of incoming flow direction and magnitude, which is essential for flight stability and control purposes.

On the other hand, in the context of flight stability, researchers have attempted to quantify changes in lift and drag on flapping wings in forward [35–38,40] and lateral flight [35,39]. In two studies [35,36], researchers have proposed that the wings act as a source of drag that is linear in body velocity or air speed owing to the averaged drag on two half-strokes. The linear frontal drag assumption was employed in the experimental study on *Drosophila* flight in [37] and showed a reasonable agreement with the data. For flapping wings experiencing lateral wind disturbances, limited experimental evidence on lateral drag was given in [35] from a wind tunnel test, suggesting it could be linear in air speed similar to the case of frontal drag. In multiple examples [38–40], additional lift and drag from forward or lateral motion were calculated from the Navier–Stokes equation at a specific airspeed and used as a baseline for stability analysis without providing an explicit relationship to the airspeed as provided by [35,36].

### 1.2. Challenges, contribution and outline

To demonstrate stable flight for the insect-scale flapping-wing robot in figure 1 under the influence of wind gusts, we require an understanding of the force and torque contribution from wind gusts on the dynamics of the robot and design or modify the flight controller to compensate for disturbances and retain flight stability.

As outlined above, existing literature describing the aerodynamics of wind disturbances for flapping-wing flight predominantly emphasizes instantaneous forces at a few specific operating conditions obtained from experiments [30,33,40,41] or fundamentally derived using CFD methods [32,38,39]. While these methods could provide a comprehensive understanding of the interaction between flapping wings and moving fluids, the lack of a tractable analytical model and the high computational cost of CFD render them unsuitable for flight control purposes. In addition, the assumption of linear frontal drag in [35,36] has not been extensively verified in practice. Very limited information is available on the effects on lift generation of the robot in figure 1 or for the scenario where flapping wings encounter wind gusts from the lateral direction.

Based on qualitative observation from preliminary flight experiments and previous work, in §2, we propose simple dynamic models that capture the essence of the effects from wind disturbances on flight dynamics of the robot. It is important to note that the models introduced herein are for flight control purposes, that is, we do not anticipate that they would produce results as accurate as those obtained from CFDs. However, we benefit from the simplicity.

Even with the proposed models, the robot still lacks the knowledge of the model coefficients, i.e. the current wind speed and direction. In §3, we discuss two possible control strategies, initially presented in [42], for the robot to estimate and compensate for wind effects in flight [17]. The proposed models and control methods are verified in a series of flight experiments in §4, where the robot is subject to wind gusts with constant, sinusoidal and naturalistic profiles with an instantaneous maximum speed of approximately 80 cm s^{−1}. The results presented here are an expansion of our preliminary research from [42]. In §5, we employed an identification technique to analyse the collected experimental data. The results indicate an unexpected influence of wind gusts on the robot's flight dynamics. We hypothesize that this irregularity is due to the modulation of wing kinematics in flight. To this end, we further conducted a benchtop-flapping experiment to obtain more conclusive evidence to support our hypothesis. Section 7 contains our concluding remarks.

## 2. Modelling the robot's flight dynamics and the effects of wind gusts

The flight-capable 80 mg flapping-wing robot in figure 1 was designed after its underactuated predecessor in [43]. The current prototype, with a wingspan of 3.5 cm, is fabricated using laser micro-machining and the *smart composite microstructures* (SCM) process as detailed in [10]. Two piezoelectric bending actuators serve as flight muscles. When a voltage is applied across the piezoelectric plates, it induces approximately linear motion at the tip of the actuator, which is transformed into an angular wing motion by a spherical four-bar mechanism. Each actuator is capable of independently driving a single wing. In operation, the robot is nominally driven with sinusoidal signals near the system's resonant frequency at 120 Hz to maximize the flapping stroke amplitude and minimize the reactive power expended by wing inertia. Lift is generated as elastic flexure hinges allow the wings to passively rotate and interact with the surrounding air. Lift modulation is achieved by altering the amplitudes of the driving signals. Using different flapping kinematics outlined in [10,15], the robot is capable of generating three torques (roll, pitch and yaw) and thrust, all approximately independent of one another.

The inherent instability of the system [13,35,44] necessitates real-time feedback control for the vehicle to achieve stable flight. Owing to the lack of power source and microcontroller on the current robot prototype, a motion capture system and an offboard computer are required for feedback and control computation. The motion capture system provides real-time feedback on the position and orientation of the robot by triangulating the positions of four retro-reflective markers placed on the robot. Driving signals are delivered to the robot via a tether as illustrated in figure 2.

### 2.1. Steady wind on flapping wings

During flight, our flapping-wing robot has a nominal flapping frequency of *f* = 120 Hz and a peak-to-peak flapping amplitude around 90°. Given the wingspan of 3.5 cm, typical hovering conditions lead to a Reynolds number of approximately 10^{3}, indicating that fluid inertial forces dominate viscosity. In this regime, a quasi-steady model states that instantaneous lift (*f*_{L}) and drag (*f*_{D}) varies quadratically with air speed according to *f*_{L,D} = *ρSC*_{L,D}*v*^{2}, where *ρ* is the ambient air density, *S* is the surface or aerofoil area, *C* is a lift or drag coefficient and *v* is the relative air speed [45].

To date, little is understood regarding the effect of steady flow on lift and drag of flapping wings. Based on empirical studies of a simultaneously translating and revolving dynamically scaled wing in [30], a revised quasi-steady model was proposed. This can be summarized as
2.1where *R* is the effective wingspan, the *k _{ij}* terms are dimensionless experimentally determined constants on the order of unity,

*α*is the angle of attack,

*Φ*denotes the wing stroke angle and

*v*is the steady flow air speed as illustrated in figure 3

*a*. The dependence of the angle of attack is qualitatively similar to [45]. In the absence of external flow (

*v*= 0) or if

*k*

_{i}_{0}=

*k*

_{i}_{1}=

*k*

_{i}_{2}, the lift and drag in equation (2.1) become quadratic functions of the flow velocity at point

*R*from the wing root For a flapping wing (instead of a revolving wing) facing frontal wind as shown in figure 3

*a*, this can be re-written as vectors to take into account the direction with respect to the body frame as 2.2To understand the lift and drag characteristics of a flapping wing encountering frontal wind, we numerically evaluate the stroke-averaged lift and drag. This is achieved by neglecting the

*k*

_{20}term (which implies that the wing experiences no drag when the angle of attack is zero) and further assuming a nominal wing trajectory of a single wing to follow with the corresponding angle of attack and

*R*= 1 cm. A similar methodology can be applied for the case of lateral incoming wind (along the axis in figure 3

*a*). For a flapping-wing robot, we compute the stroke-averaged forces contributed by both wings as a function of the air speed of the incoming frontal and lateral flow. The results shown in figure 3

*b*,

*c*suggest that incoming air contributes to increased lift, whereas drag is found to be linearly proportional to the air speed. The force profiles are similar for both frontal and lateral wind, but the effects of frontal wind are more pronounced. The linear relationship between the air speed and frontal drag was previously suggested in [35,36] and investigated more thoroughly in [41], particularly at low advance ratios. Similarly, limited empirical evidence from [35] hints the possibility of linear lateral drag, but it was not explored further.

### 2.2. The effect of wind disturbance on flight dynamics

To describe the dynamics of the robot in the presence of wind disturbance, we first define the inertial frame of reference and the body-attached frame as depicted in figure 2. A rotation matrix *R* (where ) describes the relation between the two coordinate frames. In flight, the flapping motion creates thrust (*Γ*) that nominally passes through the centre of mass and aligns with the of the robot. Additional aerodynamic damping forces caused by the external wind disturbance are denoted by a vector *f*_{w}. The resultant equation of motion describing the translational dynamics of the robot is given as
2.3

To change position, the robot has to rotate its body, so that the thrust vector along takes on a lateral component, as seen in files [46]. The rotational dynamics of the robot is given by Euler's equation
2.4where *τ*_{c} is a 3 × 1 vector of torques generated by the flapping motion as commanded by the controller, and *τ*_{w} represents the torque contributed by the effects of wind disturbances on the robot.

It has been shown that for a flapping-wing system at this scale, body oscillation at wingbeat dynamics can be neglected [38]. For control purposes, we are interested in stroke-averaged force and torques. The relatively short timescale of the actuator dynamics (*τ* ∼ 0.01 s [47]), compared with the robot's dynamics (*τ* ∼ 0.1 s [36]) means that the actuator dynamics can be ignored and we may treat the thrust and torques as four inputs to the system as described by equations (2.3) and (2.4).

Owing to experimental constraints, in this work, we limit the study to the case of horizontal wind only. Previous studies indicate that parasitic drag on the body is two orders of magnitude smaller than the contribution from flapping wings, rendering them negligible in flight [40,48]. Based on previous findings [30,35,41] and the analysis in §2.1, we predict that when horizontally moving air encounters the robot, it gives rise to a drag force linearly proportional to relative air speed. This assumption was shown to be approximately valid for a flapping-wing robot at this scale for both frontal wind and lateral wind at low wind speeds (less than 1.5 m s^{−1}) [35]. Horizontally moving air may also result in additional lift in the direction along the of the robot. As a result, we propose that the effects of a wind disturbance on the translational dynamics of a flapping-wing robot takes the form
2.5where is a vector of relative air velocity expressed in the inertial frame. For a hovering or near-hovering condition, the robot's velocity can be neglected and **v** may be taken as the disturbance wind velocity. Equation (2.5) assumes different force coefficients (*b _{x}*,

*b*,

_{y}*b*) for drag along different body axes, consistent with the numerical results shown in figure 3

_{z}*b*and the wind tunnel test in [35]. The addition of the last term stems from our attempt to model additional lift caused by frontal or lateral wind. This term is always positive and assumed to be a linear function of the air speed on the plane for simplicity. This might not be accurate according to the numerical results above, but the resultant model significantly simplifies the control strategy presented later in §3.

Because *f*_{w} does not necessarily pass through the centre of mass of the robot, it is anticipated to perturb the rotational dynamics of the robot as well. We model the torque contribution from the wind disturbance on the robot as the following
2.6where *a _{x}* and

*a*are corresponding rotational damping coefficients which are composed of

_{y}*b*and

_{y}*b*terms from equation (2.5), respectively. Equation (2.6) can be interpreted as the product of

_{x}*f*

_{w}with appropriate effective moment arms. The absence of torque along the is based on an assumption regarding the symmetry of the robot and its nominal wing kinematics. It is important to mention that equations (2.5) and (2.6) are modelled based on limited evidence. While their validity is debatable, the primary goal is to capture the dominant effects observed in the experiments and are sufficiently accurate for controller design purposes as will be demonstrated.

## 3. Flight control strategies

### 3.1. Adaptive tracking flight controller

In this paper, we present two wind disturbance rejection schemes and demonstrate how they can be implemented on the existing adaptive tracking flight controller presented in [17]. The controller is comprised of two primary components: an attitude controller and an altitude controller, operating in parallel. The attitude controller determines the required torques the robot has to generate to stabilize its orientation, track the predefined trajectory and minimize the unnecessary body rotation, whereas the altitude controller calculates a suitable thrust for the robot to follow the desired height. More detail on the flight controller can be found in [17] and the electronic supplementary material.

In the presence of wind disturbances, the disturbance terms from equations (2.5) to (2.6) appear (in the dynamics equations (2.3)–(2.4)). The exact value of these terms are difficult to determine as they depend on the air speed, wing geometry, wing trajectories, etc. Because these terms are unknown, the disturbance cannot be immediately compensated for by the flight controller and the stability of the system is no longer guaranteed unless a proper estimation or rejection scheme is implemented.

Fortunately, in the current implementation of the flight controller, force disturbances are automatically compensated for without the need for a rejection scheme. This is because the existing control law takes into consideration the acceleration error. To elaborate, the drag force from the wind would result in a non-zero acceleration position error. This prompts the controller to apply a corrective torque command to rotate the robot until the thrust vector cancels out the corresponding drag, attempting to bring the robot into the force equilibrium condition (see also the electronic supplementary material). Nevertheless, the torque equilibrium state is not sufficiently met, and a disturbance rejection scheme is required to maintain a stable flight.

### 3.2. Wind disturbance rejection schemes

Two disturbance rejection schemes for estimating the disturbance torque in this work were initially presented in [42]. The first scheme, adaptive estimation, is derived, using a standard Lyapunov-based adaptive control method [42,49]. The second scheme attempts to model the temporal structure of the disturbance using a finite impulse response filter. Their mathematical descriptions are briefly given in the following sections.

#### 3.2.1. Adaptive estimation

One strategy to counteract wind disturbances is to fully embrace the constant or slowly time-varying wind assumption, and adaptively estimate this effect based on feedback. This is achieved by re-writing the torque disturbance from equation (2.6) as
3.1In this form, the torque disturbance in the body-fixed frame is expressed as a 6 × 3 matrix made of components from a rotation matrix *R*, multiplied by a vector with six unknown elements. In this work, we restrict to the case where *v _{z}* = 0, leaving only four unknown parameters in equation (3.1). Because previous work and the numerical analysis in §2.1 suggests that the effects from frontal wind and lateral wind may be in the same order of magnitude [35], we further simplify the estimation by hypothesizing that

*a*and

_{x}*a*are also comparable in magnitude such that

_{y}*a*≈

_{x}*a*≈

_{y}*a*. This reduces equation (3.1) to 3.2where we have defined a 2 × 2 matrix

*Y*and a 2 × 1 vector of unknowns

*Θ*accordingly.

*R*represents a component in the rotation matrix. The adaptive estimation scheme assumes that the wind velocity is constant or slowly time-varying in the inertial frame, resulting in a constant unknown

_{ij}*Θ*. Equation (3.2) implies that, under constant wind disturbances, the effects on the robot's dynamics do not necessarily appear constant from the robot's perspective as they are coupled to the orientation. However, it can be decoupled for estimation purposes as shown below.

Let the notations and denote an estimated quantity and the estimation error (i.e. ). The adaptive scheme derived in the electronic supplementary material uses a Lyapunov analysis to provide an update law of the estimation, or to ensure that the estimation error asymptotically goes to zero over time. In the meantime, the flight stability is retained given that the controller compensates for the torque disturbance, using the current estimate by augmenting the term to the existing control law.

A similar strategy could be implemented to cope with the disturbance term in the altitude dynamics. Because the altitude controller, taken from [17], employed in this work already possesses the ability to adaptively estimate and compensate for constant thrust offsets, it benefits from the fact that the effect of the constant wind disturbance is identical to an unknown thrust offset. Essentially, we do not need to explicitly amend the altitude controller to deal with horizontal wind, except for changing the adaptive gain to be sufficiently large in order for the estimates to converge quickly.

#### 3.2.2. Least-squares estimation

While the adaptive estimation method is simple and easy to implement, it relies heavily on the assumption that the wind disturbance is constant or slowly time-varying.^{1} Hence, we propose an alternative approach that relaxes the constant assumption slightly, allowing the adaptive algorithm to capture, to some extent, the temporal structure of the wind disturbance. The proposed strategy resembles the adaptive estimation method in that it first estimates the disturbance and applies a counteracting input to eliminate the effects from the disturbance.

To derive the least-squares (LS) estimator for the attitude controller, we obtain an explicit expression of the disturbance torque *τ*_{w} in the rotational dynamics in equation (2.4) by re-arranging the term to obtain
3.3If the quantities on the right-hand side of equation (3.3) are measurable, then we theoretically obtain the current value of *τ*_{w}. However, a memoryless estimate obtained this way is highly susceptible to disturbances and measurement noises in practice. Moreover, to evaluate some quantities, such as on the right-hand side of equation (3.3) requires taking time derivative of the measurements. The result is prone to be highly noisy (and also non-causal).

To overcome the issues, we instead consider the estimate of *τ*_{w} after passing through a simple low-pass filter. First, we define a discrete time transfer function of a low-pass filter as , where *z*^{−1} is a delay operator and is a filter weight. Let *Ψ* denote a quantity on the right-hand side of equation (3.3) after applying the low-pass filter. Equation (3.3) becomes
3.4where [*t _{i}*] is an index for the discrete time domain (i.e.

*t*−

_{i}*t*

_{i}_{−1}is a sample time). That is, using the value of

*Ψ*calculated at the current and previous time step, we recover the estimate of

*τ*

_{w}at the previous time step.

As explained earlier, while the wind disturbance may be constant or slowly time-varying in the inertial frame, the instantaneous disturbance torque perceived by the robot varies at the same timescale as the robot's rotational dynamics. Therefore, it is sensible to consider *τ*_{w}, using the expression in equation (3.2) or which leads to
3.5To obtain the estimate of we assume that may possess some temporal structure. This inspires us to estimate via a *finite impulse response* (FIR) filter as [50]. Here, *N* is the filter length, *k* can be regarded as an arbitrary step size, and *σ _{k}*s are corresponding coefficients to be found. Then, an algorithm for a standard least-squares estimator with an exponential forgetting factor (

*ξ*) is used to determine the

*σ*s that minimize [51]. The current estimate depends entirely on these

*σ*s and

*Θ*s in the past, effectively is a weighted average of past measurements. This approach enables to possess some temporal structure. Because, this is an iterative algorithm,

*σ*s are updated every time step. The controller then compensates for the torque disturbances by projecting back to the appropriate direction as The full derivation of the scheme can be found in the electronic supplementary material and in [42].

The same strategy can be applied to the altitude controller in an attempt to correct for the time-varying disturbance force on the altitude dynamics. To elaborate, the unknown offset is treated as a single parameter using the FIR filter. That is, it may be written as a combination of its past values in the same way that is estimated. In this case, the least-square method is implemented without the need for projection.

## 4. Flight tests

To verify and assess the performance of two proposed disturbance rejection schemes, we performed indoor flight experiments on the insect-scale robot shown in figure 1. A custom low-speed wind generator was used to provide controlled gusts in the confined laboratory setting.

### 4.1. Set-up of the flight arena

Unconstrained flight experiments were performed in a flight arena equipped with between four and eight motion capture *VICON* cameras. The cameras provided position and orientation feedback of the robot at the rate of 500 Hz, covering a tracking volume of 0.3 × 0.3 × 0.3 m. Control algorithms were implemented on a computer running on an xPC Target (*MathWorks*) environment and executed at the rate of 10 kHz for both input sampling and output signal generation. Signals were generated through a digital-to-analogue converter, amplified by a high voltage amplifier (*Trek*) and then delivered to the robot through four 51-gauge copper wires as illustrated in figure 2.

Prior to the disturbance rejection experiments, the robot underwent an extensive characterization and trimming process to identify unknown torque offsets and stable flight. The detailed procedure is outline in [15] and the electronic supplementary material. In case of notable mechanical fatigue or damage to the robot, which occurs regularly in flight experiments (the lifetime of flexure-based wing hinges, for example, varies widely but is generally less than 10 min [52]), manual part replacements or repairs are sometimes possible. In such circumstances, the complete characterizing and trimming process has to be repeated.

### 4.2. Wind disturbance generator

We constructed a low-speed wind disturbance generator from an array of nine 12 V DC fans fitted in a 15 × 15 × 20 cm box, capable of creating wind disturbances in a horizontal plane. In steady state, the wind generator is able to consistently generate wind with the speed ranging from (20–100) ± 2 cm s^{−1} at the hovering flight setpoint 10 cm from the opening as verified by a hot-wire anemometer. The system is capable of generating a time-varying wind profile with a bandwidth of 0.15 Hz. The standard deviation of wind speed in the 14 × 12 cm horizontal plane and 6 × 12 cm vertical plane around the setpoint was found to be 2.5 cm s^{−1} (see the electronic supplementary material for detail).

### 4.3. Flight experiments with wind disturbances

To investigate the effects of horizontal wind disturbances on the flapping-wing robotic insect, we initially carried out a series of experiments with the disturbance generator programmed to produce a constant 60 cm s^{−1} wind in steady state (the chosen speed is comparable to the flight speed of free-flying *Drosophila* [53,54], in our case, this equates to approximately 20 wingspans per second). The wind is oriented along in the positive direction as shown figure 2. For a robot flapping at 120 Hz with 90° peak-to-peak flapping amplitude, the corresponding advance ratio is 0.1 [30]. In the first set of experiments (four flights), no disturbance rejection scheme was implemented. Then, seven flights (three and four) were performed with the adaptive estimation scheme and the least-squares estimator, separately. Lastly, both schemes were tested together simultaneously in three flights. In summary, four sets of experiments were performed.

The time-trajectory plots of all 14 flights are shown in figure 4*a*, with the averaged wind profile at the setpoint from six trials. Figure 4*a* reveals that without the disturbance rejection schemes (black lines), the robot gained significant altitude before it was blown and stayed at approximately 10 cm away from the setpoint. At this position, the robot was able to manoeuvre towards to the desired altitude, possibly owing to the reduced wind speed farther away from the setpoint.

The plots show that when the disturbance rejection schemes were in place, the initial part of the trajectory did not change significantly, nevertheless, the robot managed to manoeuvre back to within a few centimetres of the setpoint in the and a few millimetres in the shortly after. Quantitatively, we compare the performance of the proposed disturbance rejection schemes by computing the *root-mean-square* (RMS) values of the position error of the robot (in all and directions) in the last second of flight (*t* = 6.0 s to *t* = 7.0 s in this case) as shown in figure 5. The plot verifies that both proposed schemes significantly reduced the RMS error from the no compensation case. Specifically, when the least-squares method was implemented, the error decreased by more than 50%.

Next, we repeated the flight experiments using gusts with constant speed of 80 cm s^{−1} (the corresponding advance ratio is ≈0.14). In this condition, the robot generally could not sustain flight in the experimental volume without a compensation scheme and oftentimes resulted in crashes. We executed four flights with the combined adaptive and LS compensation method. To prevent the robot from crashing early in the flight while the estimates from the adaptive and LS schemes had yet to converge, we included the initial estimate of the torque contribution from the wind disturbance (*av _{x}* from equation (3.2)) to be the approximate value corresponding to 60 cm s

^{−1}(≈4.5 × 10

^{−7}Nm) as found from the earlier flights. Snapshots of an example flight taken from a high-speed video camera are compiled and shown in figure 6. The trajectory plots of all four flights performed under this condition are displayed in figure 4

*b*. The plots show that the robot could reach the setpoint within the first 3 s of flight and steadily maintain the position thereafter. Figure 5 reveals that the RMS position error is only 0.5 cm, significantly lower than the previously obtained values from flights with 60 cm s

^{−1}gusts, as the initial estimates allowed the rejection schemes to converge more quickly.

To examine further how effective the proposed strategies are in a simulated environment with time-varying wind gusts, we designed two time-varying disturbance profiles. The first one is sinusoidal, and the second one is based on a well-known turbulence model—the Dryden model [55]. More information on these wind profiles is provided in the electronic supplementary material. Selected flight footage is available as electronic supplementary material, video.

In each case, we performed flight experiments without the disturbance rejection scheme and with the combined adaptive and LS scheme. The trajectory plots are shown in figure 4*c*,*d*, with the summary of the RMS position error given in figure 5. In the absence of the disturbance rejection algorithms, the robot occasionally crashed mid-flight (this can be identified as constant position in figure 4*c*,*d*). With the proposed strategy, flight performance radically improved, resulting in no crashes in any attempted flights. Figure 5 confirms that the robot tracked the position setpoint with comparable error to the constant 60 cm s^{−1} wind case. This provides evidence that the proposed rejection schemes are capable of stabilizing the robot in the presence of slowly varying disturbances.

## 5. Identification of flight dynamics

Thus far, we have demonstrated stable flight in the presence of wind disturbances. With the adaptive technique, stable flight was achieved without precise knowledge regarding damping force and torque coefficients (*b _{x}*,

*b*,

_{y}*a*,

_{x}*a*). It is important to clarify that the fact that stable flight was obtained merely indicates that the proposed dynamic model could capture the essential features of the disturbance, allowing the robot to compensate for their effects. It does not, for instance, prove that our model is accurate and representative of all flight conditions. To gain better insights into the effects of wind gusts on flapping flight dynamics, the collection of acquired trajectory data, spatial distribution of wind speed over the flight arena and control commands can be post-processed using identification techniques similar to [56,57].

_{y}### 5.1. Identification models

#### 5.1.1. Translational dynamics

To examine the influence of time-averaged aerodynamic drag on the flapping-wing robot, we focus on the frontal and lateral direction of the translational dynamics with respect to the robot's body-fixed coordinates. In the limit near the hovering condition (see the electronic supplementary material for the definition and justification), that is the instantaneous speed, acceleration and angular velocity are relatively insignificant, we anticipate the aerodynamic drag to be balanced out by gravity such that and According to our proposed model in equation (2.5), we expect
3.6where **v**_{w} is the vector of wind velocity with respect to the inertial frame. Note that the quantities and represent the frontal and lateral air speed perceived by the robot. That is, the ground speed of the robot which was previously assumed negligible in the controller design, has been taken into account. If the assumption of linear drag along the frontal and lateral direction in equations (2.5) holds, the plots of the gravity components (*gR*_{31}, *gR*_{32}) against their respective air speeds should be linear, with the slopes corresponding to the drag coefficients (*b _{x}*,

*b*).

_{y}#### 5.1.2. Rotational dynamics

Similar to the translational dynamics, the rotational dynamics about the pitch and roll axes of the robot from equation (2.4) can be combined with the aerodynamic drag model from equation (2.6). Again, if we limit our consideration to the near-hovering case ( and see also the electronic supplementary material), it becomes 3.7With the knowledge of the torque commands and the trajectory of the robot from the motion capture system, the quantities from the right-hand side can be determined. They are expected to be a linear function of the frontal and lateral air speeds according to the proposed model.

### 5.2. Identification results

To perform the identification of flight dynamics, we post-processed 17 flight trajectories. Only portions of stable flight after taking off and before landing were selected, totalling more than 62 s or 7500 flapping periods. The data were applied with an acausal low-pass filter with the cut-off frequency of 50 Hz (to get rid of the body oscillation at the flapping frequency of 120 Hz, but the dynamics at approx. 10 Hz or lower is preserved) to eliminate measurement noise and downsampled from 10 to 1 kHz, resulting in over 62 000 data points. Numerical derivatives of position and orientation were calculated to provide translational and angular velocities, and translational and angular accelerations.

To begin, we analyse the frontal drag data by plotting the respective gravity component, *gR*_{32}, against the incoming air speed as suggested by equation (3.6) in figure 7*a*. Each grey point belongs to one data point. To ensure that equation (3.6) is approximately valid, data points near the hovering condition are highlighted in colour. For the rotational dynamics about the pitch axis, which are affected by the frontal air, the torque components on the left-hand side of equation (3.7) are plotted against the frontal air speed in figure 7*b*. Again, coloured points belong to the near-hovering condition. In both cases, the drag force and torque appear somewhat nonlinear with respect to the perceived air speed, The observed nonlinearity, however, does not conform with the quadratic drag profile as one may anticipate for the case of a non-flapping wing. Instead, the drag appears to increase less dramatically at air speeds greater than 0.5 m s^{−1} (*J _{a}* ≈ 0.08). Nevertheless, to correlate the data to the proposed linear model, we constructed best-fit lines via the least-squares method, using only data points with the air speed lower than 0.4 m s

^{−1}(shown in green). The best-fit lines produce the

*R*

^{2}-values of 0.93 and 0.88 for the translational and rotational directions, indicating a reasonable goodness of fit for the model in the range of interest.

The resultant drag force and torque contributed by the air speed in the lateral direction of the robot plotted in figure 7*c*,*d*, on the other hand, do not immediately suggest nonlinear behaviour. Particularly for the translational drag, the data suggest that it could be highly linear with respect to the air speed up to 0.8 m s^{−1}. Nevertheless, the result could be regarded as inconclusive owing to the lack of incoming wind from the opposite direction. The *R*^{2}-test for the linear fit yields the value 0.75 and 0.59, respectively.

Linear fits from the identification results here provide the estimates of translational and rotational drag coefficients as and The differences in magnitude between *b _{x}* and

*b*and

_{y}*a*and

_{x}*a*are around 20% and 35%, respectively.

_{y}### 5.3. Interpretation of the identification results

It can be seen that the estimates of translational drag coefficients *b _{x}* and

*b*from the experimental data are comparable in magnitude. This agrees with the rudimentary findings from the wind tunnel experiments from [35]. The modified quasi-steady model in equation (2.2) and the subsequent numerical simulation in figure 3, nonetheless, predict a wider difference between

_{y}*b*and

_{x}*b*. For example, if we assume that the simplified quasi-steady model used to generate figure 3 is valid, would predict

_{y}*k*

_{10}= 2.2, whereas would predict

*k*

_{10}to be 0.70. In terms of the normalized drag coefficient used in conventional quasi-steady models for translational or rotating wings, these

*k*

_{10}s correspond to and 1.4 sin

^{2}

*Φ*, comparable to the drag coefficient of a similar (but not identical) robotic wing in [54] or a

*Drosophila*wing [58]. Apart from the validity of the simplified model, the discrepancy could come from the assumption on nominal wing kinematics, the neglect of aerodynamic–wing interactions, etc.

More interestingly, the evident nonlinearity of the observed frontal drag at low advance ratios (*J _{a}* < 1) was not theoretically predicted by the quasi-steady model nor observed in related experiments [30,35,41]. We hypothesize that this nonlinear behaviour is caused by the modification of wing kinematics in flight as the robot attempted to hover while subject to the disturbance torque. To elaborate, when a flapping-wing robot with the centre of mass situated below the centre of pressure encounters gusts incoming from the frontal direction, it has to shift the mean flapping stroke dorsally to produce a pitch torque in order to counteract the torque induced by the gusts. This shift in the mean stroke position reduces the effective wing area that is perpendicular to the wind, which is likely to have subsequent influence on the aerodynamics. This phenomenon has not been thoroughly scrutinized in previous studies [30,35,41].

The identification results in figure 7*c*,*d* do not immediately suggest any deviation from the anticipated linear trend. In this situation, the robot counteracts the induced drag torque by differentially varying the flapping amplitudes of both wings to produce a roll torque. It is possible that the change in wing kinematics may also affect the overall drag in a similar fashion to translational drag, but the limited amount of flight data collected does not provide sufficient range and resolution for us to observe the change. To further investigate the possibility of nonlinear drag profile owing to the adjustment of wing trajectories, we turned to a static platform equipped with force sensors for more accurate measurements in controlled experiments.

## 6. Benchtop-flapping experiments

To investigate the influence of constant airflow on drag force production of a flapping robot, we mounted a single-winged robot on a dual-axis force sensor and measured its performance against the following parameters: mounting orientations, flapping kinematics and airflow conditions.

### 6.1. Bench top experimental set-up

Owing to the short-lived nature of the robot, we were unable to use the robot from the flight experiment in this set-up. We fabricated a similar single-winged robot half with a revised wing geometry that offers an enhanced aerodynamic performance according to [59] for static tests and airflow tests for ease of mounting to the force sensor. We can estimate force production for a complete robot, because the robot halves are symmetrical. Figure 8*a* shows the robot mounted on the force sensor for lateral and vertical force measurements [60]. Detail on the experimental procedure is given in electronic supplementary material.

Figure 8*b*–*d* illustrates three mounting configurations, wind directions and flapping kinematics of the frontal and lateral wind tests. For the frontal wind test, shown in figure 8*b*, we mounted the robot half such that the wing was perpendicular to the drag axis. The incoming frontal wind was directed in the drag axis. The bottom diagram shows the corresponding variation of mean stroke offset for the frontal wind experiments. Figure 8*c*,*d* shows the set-up for lateral wind tests. The wing driver was mounted such that the mean stroke was parallel to the drag axis. The incoming lateral wind was directed in the direction parallel to the drag axis. We positioned the wind generator to direct flow from both directions. The bottom diagrams illustrate how the wing stroke amplitude was varied for the lateral wind experiments.

### 6.2. Flapping experiments

To assess the effects from the mean stroke offset and amplitude variation from flapping wings in the presence of 60 cm s^{−1} frontal and lateral wind on drag, we performed five sets of experiments. In the first two sets, one with frontal wind and one without wind, we measured the drag force while varying the mean stroke offset (11 values) as depicted in figure 8*b*. The difference in drag measurements between these experiments represents the effect of wind on the flapping wing alone (i.e. excluding the drag from the body or structural components of the sensor). The other three sets were performed with the robot mounted in the orientation shown in figure 8*c*,*d*: one without airflow and two with the airflow from the wing root/tip towards the wing tip/root. In this case, the flapping amplitude was varied (19 values). Measurements from the three sets can be processed and combined to simulate the scenario where the robot adjusts the flapping amplitude differentially to generate roll torque to stabilize against the lateral wind in flight (see top diagrams in figure 9).

At each data point, the robot half was commanded to flap for 5 s at 130 Hz. Middle portions of the data lasting 1.3 s (approx. 160 cycles) were processed with a low-pass filter with the cut-off frequency of 600 Hz to eliminate measurement noise. Measurements from five experimental sets are integrated to simulate the in-flight wing kinematics in response to the airflow as we described earlier under the assumption that the aerodynamic interaction between the two wings and the body is negligible. The results are shown in figure 9.

According to figure 9*a*, when a pair of flapping wings in the neutral position (no offset angle, peak-to-peak amplitude of 60°) encounters 60 cm s^{−1} frontal gust, we measured nearly 70 µN of drag force. This diminishes rapidly when the offset angle was introduced to mimic the response of the robot in flight. The measured drag almost disappeared as the offset angle reached 5° in both positive and negative directions.

Figure 9*b* reveals that when a pair of flapping wings (no offset angle, mean peak-to-peak amplitude of 55°) faces the incoming airflow from the lateral direction, the measured drag force varied as a function of the peak-to-peak differential amplitude. At small differential amplitude values (less than or equal to 10°), the variation in measured drag is relatively small, of the order of 10 µN. It was found that drag force reduces drastically at large values of differential amplitude.

### 6.3. Quasi-steady model predictions and analysis

To supplement the benchtop-flapping experiments, we employed the simplified quasi-steady model introduced in §2.1 to numerically predict the resultant drag on a pair of flapping wings subject to frontal and lateral wind. Similar to the flapping experiments, offset angle was varied in the presence of frontal wind. For the lateral wind, flapping amplitudes were differentially altered. Nominal wing kinematics (such as flapping amplitude, angle of attack, frequency, etc.) were identical to those in §2.1. While these parameters might not match our experimental conditions exactly, we believe they should qualitatively capture the trend. The numerical results are shown at the bottom of figure 9.

In the frontal wind case, the numerical simulation generally agrees with the flapping experiment. That is, a notable reduction in drag force was seen at non-zero offset angle. The experimental data demonstrate substantially higher sensitivity to the offset angle than the prediction from the quasi-steady model. However, it is qualitatively consistent with the identification result in figure 7 (where the impact of frontal wind on drag was seen decreased at greater wind speed) in regards to the hypothesis that the robot adjusted the wing trajectories to stabilize its attitude against the drag torque.

The numerical simulation also predicts very little variation in the drag force when the flapping amplitude was varied differentially in the presence of lateral wind (figure 9*b*, bottom). This is because of the antagonistic effect from two wings. On the other hand, the experimental data hints that, when the difference in amplitude is small (less than or equal to 10°), the lateral drag may not be symmetric about positive and negative differential amplitude values as suggested by the quasi-steady model. We believe this is likely owing to the disruption of airflow from the wing closer to the wind source preventing the flow reaching the other wind with the same speed. As a consequence, the observed drag is greater when the nearer wing has a larger stroke amplitude. At a large amplitude difference (approx. ±25°), an appreciable drop in drag force is seen. We are unable to offer a concrete explanation for this phenomenon, but factors related to passive wing rotation, complex interactions between the dynamics of actuators, wings and aerodynamic forces could be responsible for the observed behaviour.

## 7. Conclusion and discussion

In this paper, we study the effects of horizontal wind disturbances on the flight of an insect-scale flapping-wing robot. Based on limited previous studies and experimental evidence, we proposed simple models to capture the effects of wind disturbances on the translational and rotational dynamics of the robot for control purposes. With a few simplifying assumptions, two disturbance rejection schemes compatible with the adaptive tracking flight controller previously developed in [17] were presented. The strategies were implemented and verified in a series of flight control experiments, including gusts with constant and time-varying wind speeds with air speed up to 80 cm s^{−1} (advance ratio *J*_{a} = 0.14). Predictably, we found that the initial part of the trajectories, when the disturbance rejection schemes were implemented, does not differ from when the schemes were absent. This is due to the fact that the proposed schemes were not aware of the presence of the disturbance at the beginning and they require some time for the estimates to converge. With partial prior knowledge of the wind speed and direction, the robot rapidly stabilized and gave rise to markedly smaller position error as highlighted in figure 5. Overall, the proposed schemes could prevent the robot from crashing and significantly reduced the position error to the level comparable to non-disturbed flight.

A careful analysis into collected flight data revealed an unexpected behaviour concerning the interaction between flapping wings and steady airflow. We observed that the time-averaged profile of drag force or torque on the robot, particularly in the frontal or pitch direction may not be a linear function of surrounding air speed as previously suggested in the literature [30,35,41]. We speculated that the discrepancy results from the adjustment of wing kinematics during flight in order for the robot to counteract the drag torque. To this end, we modelled this phenomenon, using a simplified quasi-steady model and verified the hypothesis experimentally, using a single flapping wing mounted on a force sensor. The results appeared to strengthen our hypothesis. That is, the overall damping force is decreased when the wing kinematics are modified to reflect how the robot stabilizes its attitude in flight when there is a wind disturbance. The effect is more pronounced in the case of the robot changing its wing offset position in response to frontal wind.

From the flight control perspective, this finding has an interesting implication. In some aspects, the effort of our flapping-wing robot, or flying insects alike, to stabilize flight in the presence of a disturbance itself indirectly alleviates the effect of wind gusts. Potentially, this allows flying machines to withstand disturbances of greater magnitudes than we originally predicted.

The experiments in this work were carried out on an insect-scale robot with passive wing hinge. There has been evidence suggesting that flying insects primarily rely on passive wing rotation, with occasional changes in the rest angles to generate torques for enhanced manoeuvrability [61,62]. In this study, we assumed a prescribed time-varying angle of attack in the quasi-steady model and have not experimentally considered the precise role of passive rotation or wing flexibility on the dynamics of the wing motion and robots or how this affects the average drag upon encountering the wind disturbance. This is one of the avenues yet to be explored in future work.

## Authors' contributions

P.C. conceived the study. P.C., Y.C., E.F.H., K.Y.M. and R.C. set up the experiments. P.C., Y.C. and E.F.H. performed the experiments. P.C. analysed the data and prepared the manuscript. R.J.W. supervised the study. R.J.W. and K.Y.M. edited the manuscript. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This work was partially supported by the National Science Foundation (award no. CMMI-1251729), and the Wyss Institute for Biologically Inspired Engineering.

## Disclaimer

Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## Acknowledgements

We wish to thank N. Gravish for comments, discussions, and his assistance with videography during the experiments.

## Footnotes

One contribution of 19 to a theme issue ‘Coevolving advances in animal flight and aerial robotics’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3575855.

↵1 Slowly time-varying, in this case, is with respect to the dynamics of the robot. This is so that the algorithm could distinguish the disturbance dynamics from the robot's dynamics. Fortunately, dynamics of natural gusts (approx. 0.1 Hz [50]) are typically slower than our robot's dynamics (approx. 10 Hz).

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.