System identification of wind effects on multirotor aircraft

Multirotor airplanes are widely used in many outdoor applications, e.g., agriculture, transportation, and public safety, where winds might be strong and prevalent. However, the effects of wind on multirotor aircraft are still not fully understood yet. The objective of this paper is to investigate and model wind effects on a real hovering octocopter. The wind is directly measured and considered as one of the inputs to the bare-airframe model. Then a series of models, each corresponding to a different wind condition, are identified from real flight data through a system identification approach. The time-domain validation results show that an average of 15% error reduction can be achieved by considering wind effects, captured by a wind correction term. The identified models will play an important role for the future development of model-based controllers for outdoor multirotor aircraft.


Introduction
Multirotor unmanned aerial vehicles (UAVs) are becoming increasingly popular over the past few years due to their high maneuverability and vertical takeoff and landing capability. They have been widely used in many outdoor applications, such as the aerial photography, infrastructure inspection, agriculture monitoring, transportation, public safety, and delivery services (Valavanis and Vachtsevanos 2015;Iost Filho et al. 2020). When flying in outdoor environments, such multirotor airplanes may be subject to the influences of various environmental disturbances, particularly winds. A strong wind may cause dramatic changes in the forces and moments acting on a multirotor UAV, which, if not taken care of properly, can significantly degrade the overall UAV flight performance. Therefore, the effects of wind on multirotor UAVs need to be formally investigated and modeled, which is the objective of this paper.

Related work
Existing approaches to analyze and model aerodynamic effects, including wind effects, on multirotor aircraft can be roughly classified into the following four categories: 1. First-principle-based methods: The blade element theory and the momentum theory are widely used to derive the forces and the moments generated at the rotor plane of a multirotor airplane (Martin and Salaun 2010;Bristeau et al. 2009;Hoffmann et al. 2007). The blade flapping effect, which is significant in helicopters, has also be studied and adapted for multirotor aircraft (Hoffmann et al. 2011;Pounds et al. 2010;Mahony et al. 2012). These first-principle approaches all make assumptions that are mostly derived from single-rotor aircraft. Moreover, these approaches regularly neglect the interferences between different rotors, as well as between the rotors and the vehicle airframe. Flight experiments are required to validate such assumptions and subsequently the derived models. Recent work in Jain et al. (2019) presents a way to model the downwash effect during proximity flight of two quadrotors by using the velocity field model and the result is validated in flight experiments.

Computational fluid dynamics (CFD) simulations:
High-fidelity CFD simulations have been carried out to investigate aerodynamic interactions on a small quadrotor in Yoon et al. (2017). However, such simulations are computationally expensive and cannot directly produce an analytic model, necessary for the purpose of either dynamic characteristics analysis or controller design. 3. Wind tunnel experiments: NASA Ames has conducted several wind tunnel experiments to determine the forces and the moments of a range of commercial multirotor airplanes under different wind speeds and attitudes (Russell et al. 2016). However, the cost of such wind tunnel experiments is very high. Furthermore, at least with the wind tunnel used in Russell et al. (2016), only static forces and moments can be measured at fixed conditions. Dynamic characteristics still remain indeterminable from such wind tunnel experiments and free flight tests are required to obtain the actual dynamic aircraft model. 4. System identification: System identification has been used successfully in the modeling of full-scale fixedwing and rotary-wing airplanes (Tischler and Remple 2012). The to-be-identified airplane needs to execute designated maneuvers to enable effective excitation on its dynamics. Then the airplane's model parameters can be estimated by analyzing the gathered flight data in either a time domain or a frequency domain (Valasek and Chen 2003;Grauer and Morelli 2014;Morelli 2012Morelli , 2000. This approach has also been applied to the modeling of small UAVs (Mettler et al. 1999;Wei et al. 2015;Gong et al. 2019), which usually have unique configurations and the aerodynamic interactions between their various components can be very complicated.
When evaluating the turbulence effects on fixed-wing airplanes, the Dryden spectral model, which assumes a "frozen-field", is normally used (Hale et al. 2017). However, the Dryden model is not suitable for a lowspeed rotorcraft as the mean wind speed becomes dominant (Dahl and Faulkner 1978). The authors in Lusardi et al. (2004) proposed a method called control equivalent turbulence inputs (CETI) to model the turbulence on a low-speed helicopter and treated the turbulence as an equivalent control input. This approach has been proven to capture the effects caused by the turbulence on a UH-60 helicopter in Lusardi (2004) and a 3DR Iris+ quadrotor in Juhasz et al. (2017). Nevertheless, the wind speeds in these cases were not measured directly but inferred. More importantly, they assumed the mean wind is constant temporally and therefore didn't consider the effects of wind changes on the rotorcraft dynamics.

Contributions
In this paper, we use the system identification approach to investigate and model the wind effects on an octocopter platform based on data collected from real flight experiments. Wind speeds and orientations are directly measured from a high-precision wind sensor and treated as an additional input entering the bare-airframe model. A series of linearized state-space models, each corresponding to a different wind condition, are identified in the frequency domain with comprehensive identification from frequency response (CIFER®) software. In addition, the eigenvalues and modes of the identified models are calculated and compared to evaluate the wind effects on aircraft dynamics. Finally, the identified models are validated with real flight doublet tests. The major contributions of this paper are as follows: -To the best of our knowledge, this paper is the first instance where the system identification approach has been used to model the wind effects on multirotor airplanes. The wind is directly measured and considered as one of the inputs to the model. -Similarly, we believe that this paper is also the first instance where the effects of wind on multirotor airplanes, particularly with respect to different wind speeds, have been studied and validated from outdoor experiments with real wind data. We also discuss how wind speeds affect the damping ratio and the natural frequency of the pitch mode of the bare-airframe model. Such kind of insights will play an important role for the future development of model-based controllers for outdoor multirotor aircraft.

Hardware and data collection
In Sect. 2.1, we will introduce the UAV platform as well as the wind sensor used in our flight experiments. Then, in Sect. 2.2, we will introduce a working hypothesis regarding the wind distribution, which provides the necessary rationale to our experiment design. Finally, in Sect. 2.3, we will provide the details on how our flight experiment is conducted and how the data on both the UAV and the wind are collected.

Hardware description
The multirotor UAV used in our experiments is a retrofitted octocopter, of which the main body is built from a DJI Spreading Wings S1000+. The Pixhawk Cube is selected as the flight controller due to its fully programmable capability. A FrySky X8R receiver and a Taranis X9D Plus transmitter are used as the RC pair. A Here+ RTK GPS is used to provide precise position and velocity estimation of the UAV and a telemetry radio is installed to monitor the real-time flight status. Each propeller of the octocopter has a diameter of 38.1 cm and a maximum thrust of 2.8 kg. The octocopter has a motor-to-motor distance of 104.5 cm and a takeoff weight of 7.0 kg. A 16,000 mAh 6S lithium battery is used, which can provide a flight duration of around 20 minutes. A photo of our multirotor UAV is provided in Fig. 1. The wind sensor used in our experiments is a Model 91000 ResponseONE ultrasonic anemometer, which can measure horizontal winds. Winds in the vertical direction are less significant than the horizontal winds and thus neglected in our study. The sensor is mounted on top of a telescopic antenna mast and has a distance of 4.6 m (15 ft) above the ground. Figure 2 shows the setup of our outdoor flight experiments, which will be explained in detail in the next sub-section. Since the goal of this paper is to study the wind effects on multirotor UAVs, ideally we would like to measure the wind at the location of the UAV, which unfortunately cannot be realized in practice. Therefore, one working hypothesis that we are making in this paper is that after considering transport delay, the wind at the location of the UAV is roughly the same as the wind at the location of the wind sensor, as long as the UAV and the wind sensor are at the same height and not far away from each other. Such a hypothesis is in accordance with the Taylor's hypothesis, which states that an eddy, satisfying M < 0.5M with M being the wind speed and M being the standard deviation of the wind speed, has negligible change as it advects past a sensor Stull (2012).

Working hypothesis on wind distribution
In order to test our hypothesis, we use two identical wind sensors and compare their wind measurements in two different scenarios (as shown in Fig. 3): (i) the sensors are positioned along the wind direction and (ii) they are positioned perpendicular to the wind direction. The sensors are placed at the same height (4.6 m) above the ground. The wind speed and direction measurements from the two sensors are plotted in Fig. 4. A transport delay has been considered in scenario i) to account for the wind traveling time. The root mean square errors (RMSEs) between the two wind sensor readings are calculated for both scenarios and the results are provided in Table 1. It can be observed that the RMSEs are very small compared to the mean wind speed and the mean direction. Moreover, the difference along the wind direction is smaller than the one perpendicular to the wind direction.
Based on these observations, in our flight experiments, we decide to command the UAV to hover at a position a distance of 4.6 m away from a single wind sensor and along the wind direction. The measured wind speed and direction by Fig. 1 A photo of our multirotor UAV, a DJI Spreading Wings S1000+ retrofitted with (1) a Here+ RTK GPS, (2) a FrySky X8R receiver, (3) a telemetry radio, (4) a Pixhawk Cube flight controller, and (5) a battery eliminator circuit the wind sensor will be regarded as wind data at the location of the UAV after taking a transport delay into consideration. The wind sensor is placed upstream of the UAV so that the downwash from UAV propellers does not affect the wind measurements. Moreover, the UAV is commanded to fly at the same height as the wind sensor so that we can ignore the vertical distribution of the wind.

Comment
We would like to point out that, even though our hypothesis is still not ideal (ideally we would hope to directly measure the wind at the exact location of the UAV), it does reflect the fact that winds change more rapidly temporally than spatially (Stull 2012). Furthermore, our treatment of winds, whereas we only assume winds to be spatially uniform, is far more sophisticated than state of the art (Lusardi 2004;Juhasz et al. 2017), where winds are assumed to be uniform both spatially and temporally.

Flight Experiments
The goal of our modeling (system identification) effort is to identify the bare-airframe model under the effects of wind. What we meant by "bare-airframe model" is simply a model consisting of the mixer, the actuators, and the vehicle, as shown in Fig. 5. The mixer generates the command for each motor based on the controller output. The commands are then executed by actuators which include the ESCs, motors, and rotors, resulting in different forces and torques and eventually change the state of the vehicle. The bare-airframe is known to be unstable, therefore we must have the control system (see Fig. 5) active during our flight experiments. The pilot command pilot and the computer-generated command auto (to be described later) are injected into the control system, resulting in four inputs to the mixer: lat , the roll control, lon , the pitch control, yaw , the yaw control, and thr , the throttle control. Our control system is customized in order to provide a position feedback to keep the UAV roughly a 4.6 m clearance above the ground (to ensure that the ground effects do not contribute to the identified model (Wei et al. 2019). The measured wind is used as described in Sect. 2.2 and treated as an additional input, wind , into the bare-airframe model.
The onboard inertial measurement unit (IMU) measures the translational accelerations and the angular rates of the UAV, while the RTK GPS and the magnetic compass measure the global position and the heading angle of the UAV, respectively. All these measurements pass through an extended Kalman filter, resulting the follow set of estimated states: (1) the roll, pitch, yaw Euler angles { , , } , (2) the roll, pitch, yaw angular velocities {p, q, r} in the body frame, (3) the linear velocities {v n , v e , v d } in the inertial frame where the subscripts n, e, d represent for the north, east, downward directions, and (4) the linear accelerations {a x , a y , a z } in the  Table 1 The RMSEs of two sensor readings in the two scenarios shown in Fig UAV body frame where the subscripts x, y, z align with the vehicle's body axes. All these states, together with the mixer inputs { lat , lon , yaw , thr } which are the magnitudes of latitude, longitude, yaw, and throttle commands normalized to 0~100%, are logged at 100 Hz. On the other hand, the horizontal wind speeds ( v w n , v w e ) in north and east directions are recorded at 20 Hz by the wind sensor, where the superscript w represents for the variable related to wind. They are resampled to 100 Hz through the fourth-order spline interpolation to be consistent with the estimated states. Moreover, the airspeed in the body frame can be calculated by using T are the inertial velocity and the inertial wind speed vectors, respectively, and is the rotation matrix from the inertial frame to the body frame.
We apply two different types of computer-generated commands (signals) in our flight experiments: (1) a frequency sweep signal for the frequency-domain identification purpose, and (2) a doublet signal for the time-domain model validation purpose. A frequency sweep signal spans a frequency range of 0.1-10 Hz with a fade-in added in the low frequency portion to avoid sharp inputs to the motors (Tischler and Remple 2012), while a doublet signal is a two-sided pulse with the same magnitude in both the positive and the negative directions. Two sets of sample flight data, one for a longitudinal sweep and the other one for a longitudinal doublet, are shown in Figs. 6 and 13, respectively.
We call one of our data collection sessions a "trial". Each trial corresponds to one particular signal (frequency-sweep or doublet), one particular axis (longitudinal, lateral, directional, or throttle), and one particular wind condition (calm, light, and strong, more details will be provided later). Each trial starts and ends at the same trim condition and lasts 34 s [the duration is selected to capture enough spectral contents and ensure the vehicle's attitude controls do not degrade as the battery weakens (Michel et al. 2019)]. Before each trial, the magnitudes of the to-be-injected signal, whether frequency-sweep or doublet, is determined to sufficiently and effectively excite the multirotor UAV dynamics, e.g., ensuring a good signal-to-noise ratio. We have conducted 5 trials for each signal-axis-wind combination, therefore, the total trial number is 5 × 2 × 4 × 3 = 120 . To allow us to capture the effects of wind on UAV dynamics, we select three days with quite different wind conditions: a calm wind day (with an average wind speed of 0.9 m/s), a light wind day (with an average wind speed of 2.7 m/s) and a strong wind day (with an average wind speed of 5.4 m/s) and conduct a total of 5 × 2 × 4 = 40 trials on each day. Use the Beaufort Scale to categorize the wind conditions, the three wind speeds fall into force 0, force 2, and force 4, respectively. It is noted that even though we only selected three wind conditions here, the proposed method can be extended to different wind conditions. The three wind conditions we surveyed in this manuscript are the wind speeds common in the area we lived in. Higher wind speeds are possible but can rarely be seen during our data collection season, so we temporarily don't have the data for them.

System identification method
In this section, we follow the procedures prescribed in book (Tischler and Remple 2012) (for the system identification of fixed-wing and rotary-wing aircraft) to identify the wind effects on multirotor UAV, based on the data collected from our trials. Roughly speaking, we take the following two steps:

Step I: model structure determination
We use commonly used six-DOF flight-dynamics equations of motion (Powers et al. 2015) as the starting point to derive the state-space model of the bare-airframe. A first-order actuator model is explicitly considered to capture the highfrequency roll-off in angular rate frequency responses, as recommended in Gong et al. (2019): where is the original control signal, ′ is the lagged control signal after considering the first-order dynamics, and T a is the time constant. To reflect such effect of actuator dynamics, we add three additional states, ′ lat , ′ lon , and ′ thr , to the original rigid body state vector, resulting the final state vector of the state-space model being u, v, w, p, q, r, , , , � lat , � lon , � thr ] T . As for the controls of the state-space model (see Fig. 5), the control input downstream the controller is = [ lat , lon , yaw , thr ] T , while the wind is modeled as an additional input, = [ w lat , w lon , w thr ] T . Moreover, the time delay in each control axis is included into the state-space model to account for unmodeled high-frequency dynamics. Finally, we get the following state-space model: where is an identity matrix in our case, is the state derivative matrix as follows: with [U 0 , V 0 , W 0 , 0 , 0 ] T being the trim conditions (the way to obtain these conditions will be presented in Sect. 4.1), and are the control derivative matrices with respect to the controller and the wind, respectively, as follows: = [ lat , lon , yaw , thr ] T and = [ w lat , w lon , w thr ] T are the time delays corresponding to the control input and the wind, respectively, = [u, v, w, p, q, r, a x , a y , a z ] T is the observed output, and and are the measurement matrices, allowing the measurement vector to be expressed in terms of the state vector and its derivative ̇ . We can derive and using the following two observations: (1) the body translational and angular velocities can be directly found in the state vector, and (2) the accelerometer measurements are of the form: The state-space model structure (Eq. (3)) can be further reduced by analyzing the coherence of each input and output pair from the flight experiment data. We implement such a reduction procedure by following the guideline specified in Tischler and Remple (2012), i.e., where ̂2 xy is the coherence value between an input x (an element of ) and an output y (an element of ), and min and max are the minimum and maximum frequencies of the fitting range, respectively. Only those responses that meet both requirements will be used in the model parameter identification step, i.e., Step II (Sect. 3.2).
Comment A second working hypothesis that we are making in this paper is that the parameters of the and the matrices are decoupled from the parameters of the matrix. In other words, the wind effects can be identified separately from those of control input (see Fig. 5). Such a hypothesis is pragmatical in nature. We will provide the practical reason for such a hypothesis in Sect. 4.3. But the main idea is it allows us to identify and separately in Step II (Sect. 3.2). Even though the hypothesis is not absolutely correct we will show in Sect. 4 that it offers an approximation which in turn leads to a model that not only captures the wind effects mathematically but also fits flight data reasonably well.

Step II: model parameter identification
The parameters of the determined model (Eq. (3)) are estimated by minimizing the following cost function: where n TF is the number of different frequency-response pairs, n is the number of frequency points in each frequency range ( 1 , n ), W , W g , and W p are the weighting functions, and T c and T are the frequency responses obtained from the flight data and the identification solution, respectively. The cost function is derived from minimizing the summed cost for the n TF transfer functions with high coherency values and for each transfer function, we minimize the weighted sum of the square errors of both magnitude and phase at different frequency points. The optimization problem in Eq. (8) is solved by the secant method (McDaniel and Burrows 1968) which has been proven robust for the model structures involving large numbers of identification parameters and frequency-response pairs (Tischler and Remple 2012). A CR ≤ 20% (Cramer-Rao bound) and I ≤ 10% (Insensitivity) are considered as acceptable for each identified parameter. An overall average cost function that achieves J ave ≤ 150 is considered as achieving an acceptable level of accuracy.

Trim conditions
The bare-airframe model is identified at three different wind conditions, calm, light, and strong, as mentioned in Sect. 2.3. The trim states are obtained by forcing the aircraft to hover in the wind for 30 seconds without any pilot input. These trim states are averaged and presented in Table 2. Since in all our trials, the pilot commands the UAV to fly in a direction roughly facing the wind, meaning the wind is in the longitudinal axis, we have both 0 and V 0 to be zeros. Notice that a higher wind speed in the longitudinal direction increases both the trim pitch angle 0 and the trim forward speed U 0 . It can also be observed that the calm wind speed is very small, therefore we can assume all the three trim values are zeros under the calm wind condition. These trim values can be used in Eq. (4).

Results under calm wind condition
We implement the system identification process described in Sect. 3 to the data collected on the calm wind day. We find that all off-axis responses have low coherence values and thus are not used in the identification process. This also indicates that the UAV dynamic is well-decoupled. The onaxis angularg rate and accelerometer frequency responses along the longitudinal axis are provided in Fig. 7. The high coherence value indicates a good linearity of the model. Moreover, a first-order actuator model (see Eq. (1)) can successfully capture the roll-off at the high-frequency range for the angular rate response. The lateral response has a similar performance thus is not shown here. The finalized lateral and longitudinal models are and The frequency responses of the yaw and heave axes are provided in Fig. 8, which show that a first-order model is again adequate to capture the dynamics along each of the two axes. However, the yaw damping N r and the heave damping Z w both have high CR% values and thus are dropped in the finalized models. This indicates that the heave and the yaw rate dampings are not identifiable from the sweep signals we injected into the bare-airframe. Other excitations are needed to capture the value of these two terms. The equivalent time delay thr is found to be unrelated and thus taken as zero, which means that our first-order actuator model (see Eq. (1)) is sufficient to capture the delay in the vertical axis. The finalized yaw and heave models are and Equations (9), (10), (11), and (12) together complete the bare-airframe model under the calm wind condition. The identified parameters of these equations are provided in Table 3. As we can see, the CR% s and I% s for all parameters are all within the guideline prescribed in (Tischler and (9)

Results under light and strong wind conditions and identification of wind effects
In natural outdoor environments, the wind is uncontrollable, meaning, different from the computer-generated command auto that we can design and pre-determine to sufficiently excite our UAV dynamics, we cannot just "design" a wind profile, e.g., we cannot frequency sweep a wind! This is the main reason underneath our second working hypothesis (see Sect. 3.1). By assuming the parameters of the and the matrices of the state-space model (Eq. (3)) are decoupled from the parameters of the matrix, in practice, we can then decouple the identification of the and matrices and that of the matrix as follows: we can complete the first step, the identification of the and matrices (Sect. 4.3.1), by following the same steps as those used in the calm wind condition, and then complete the second step, the identification of the matrix (Sect. 4.3.2), by considering the residual after removing the effects due to the and matrices.

Identification of and
The longitudinal frequency responses under the light and strong wind conditions are provided in Figs. 9 and 10, respectively. It can be seen that the coherence in the lowfrequency range is affected by the wind. Based on these data, we follow the same steps as those used in the calm wind condition to identify a model with the and the matrices determined. The identified parameter values are provided in Table 3. It can be observed that such a model, i.e., a model with only the and the matrices, has a higher cost value than the one under the calm wind condition, indicating that the wind affects the longitudinal dynamic response. Moreover, it can be observed that the wind has a limited effect on the yaw and heave dynamics. So their responses are not plotted here. Finally, it is noted that when the wind speed is high, the yaw angular rate damping N r becomes significant. Therefore, we add N r back to the model.

Identification of
Based on our second working hypothesis, the residual, defined as the difference between the actual UAV response and the predicted response by only considering the and terms: accounts for the wind effects, i.e., the term in Eq. (3). Specifically, we perform an one-step prediction by using the partial model we have already identified in Sect. 4.3.1 (i.e., the model with the and matrices but not the matrix) and the state vector we recorded from our flight experiments. Since the body translational velocities and the angular rates are directly measured, only the accelerometer measurements are left during the residual calculations. Sample residuals of a x and a z under the light and the strong wind conditions are plotted in Figs. 11 and 12, respectively. We can then use these residuals to identify the parameters in the matrix. Due to the symmetrical configuration of our UAV platform, we assume that the wind has the same effect along the longitudinal and the lateral axes, thus we only need to identify four parameters, X w,lon , Z w,lon , w lon , and w thr . Their identified values are provided in Table 5.

Validation results
The doublet response (the data collected for such response is not used in the system identification process but only for the validation purpose) for the model identified under the calm wind condition (obtained in Sect. 4.2) is shown in Fig. 13. It shows excellent agreement between the model output and the flight data. The doublet responses for the two models identified under the light wind condition, one without the wind correction term (i.e., the term, obtained in Sect. 4.3.1) and the other with the wind correction term (obtained in Sect. 4.3.2) are shown in Fig. 14. The similar doublet responses for the two models identified under the strong wind condition are shown in Fig.15. Furthermore, the RMSE values for different wind conditions and with different models are given in Table 6. The table shows that when we consider the wind correction term, , the RMSE of a x can be reduced by 13% under the light wind and 16% under the strong wind condition. Meanwhile, the RMSE of a z is reduced by 15% under the light wind and 7% under the strong wind. All these results show that the predictive power of our models is significantly improved by considering the wind effects.

Eigenvalues and dynamic modes under different wind conditions
Based on the bare-airframe models identified for different wind conditions, we calculate their eigenvalues. These eigenvalues, together with their corresponding dynamic modes, are shown in Table 7. We can see that there are two unstable poles in both the lateral and the longitudinal dynamics, indicating that the bare-airframe multirotor airplane itself is inherently unstable, therefore demonstrating the necessity of a control system in the loop. Furthermore, it can be observed that the eigenvalues in the pitch and the roll modes are very close to each other, which indicates the high symmetry of our UAV platform (Fig. 16).
The results regarding eigenvalues and dynamic modes show that when a multirotor UAV, e.g., the octocopter used in this paper, is hovering under a windy condition (either with a light or strong wind), its dynamic characteristics are significantly affected by the wind. The wind effects can be summarized as follows: -Both the damping ratio and the natural frequency of the longitudinal dynamics increase as the wind speed increases. For instance, the natural frequency is increased  . This means that the bare-airframe becomes more unstable as the wind speed increases and a more robust control system is required to stabilize the vehicle.
-The eigenvalues in the roll mode remain unchanged with respect to the wind speed as long as the wind direction is aligned with the vehicle's longitudinal axis, indicating that the wind has a less significant effect on the dynamics in the lateral axis. -The wind has a limited effect on the yaw response of a multirotor UAV under a hovering condition. A damped yaw mode is only observed when the UAV is under the strong wind condition.
All these findings are helpful in the design of control system for outdoor multirotor UAVs.

Current limitations
Though we successfully identified the wind effects on multirotor aircraft under different wind conditions, there are currently some limitations of our proposed method. First, the additional aerodynamic forces caused by the wind are lumped together. A comprehensive understanding and modeling of different sources and components can be difficult. For instance, the change of interaction effects can not be separated from the drag acting on the rotors induced by the wind. Further work need to be done to study and model different components of the wind effects which we discuss it in the next section. Second, our method requires accurate wind measurements as the input. However, such accuracy and even the wind data may not always be available in the actual flight. One potential solution is to equip a second UAV with the high-precision wind sensor and command it to fly around the target UAV. The second UAV is able to provide real-time and accurate wind measurements at the desired location.

Possible aerodynamic explanation and future work
Our results show that some aerodynamic forces and moments are directly affected by the wind and the wind effects can be partially accounted for by considering the wind as one of control inputs entering the bare-airframe model. Such additional aerodynamic forces can arise from (1) the drag force acting on the rotor plane, (2) the thrust changes generated by the rotors, (3) the interactions between the rotors, and (4) the interactions between the rotors and the body.
Our study provides a natural starting point for the systematic understanding of wind effects on multirotor UAVs. In order to obtain a comprehensive understanding, a multipronged methodology should be adopted. First, further analysis can be conducted by combining first-principle models and our system identification results. For example, by explicitly calculating the blade flapping and the induced drag from existing literature (Johnson 2013) and comparing such first-principle-based results to flight experiment data, we can find discrepancies, which can in turn be used to improve existing models. The discovery of such discrepancies can also improve our understanding of the complex aerodynamic effects of wind on multirotor aircraft. Second, forced oscillation wind tunnel experiments (Yeo and Romander 2013) can be conducted in the future to allow systematically injecting winds with proper spectral content to sufficiently excite multirotor UAV dynamics. For instance, a frequency sweep wind can be generated by such type of wind tunnels to excite a multirotor UAV and then a similar approach as presented in this paper can be used to identify the wind effects on multirotor UAV dynamics. However, such experiments can be costly.

Conclusions
In this paper, we used a system identification approach to systematically identify wind effects on multirotor UAVs under a hovering condition. The wind was directly measured by a wind sensor and considered as one of control inputs into the state-space model. Our time-domain validation showed that such a treatment can reduce the prediction error by 15% in average. Our identified models can not only offer us unique insights into the effects of wind on multirotor dynamics but also be used directly for model-based flight control system design in the future.