**1. Introduction**

Nature provides investigators with abundant insights and inspirations that aid in solving engineering problems and creating new inventions [1,2]. With the advancement of science, engineering, and technology, a growing number of robots have been designed in different shapes and sizes to mimic biological systems [3,4]. These biomimetic robots mimic animals, such as four-legged animals (quadrupeds) [5], monkeys [6], snakes [7], fish [8,9], bats [10,11], and various birds and insects [12–14]. These robots have diverse potential

**Citation:** Aurecianus, S.; Ha, G.-H.; Park, H.-C.; Kang, T.-S. Longitudinal Mode System Identification of an Insect-like Tailless Flapping-Wing Micro Air Vehicle Using Onboard Sensors. *Appl. Sci.* **2022**, *12*, 2486. https://doi.org/10.3390/ app12052486

Academic Editors: Luis Gracia and Carlos Perez-Vidal

Received: 26 January 2022 Accepted: 24 February 2022 Published: 27 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

applications, ranging from simple toys to exploration robots in dangerous and hazardous environments. They have the potential to perform tasks that conventional robots cannot accomplish, such as climbing steep surfaces [15,16], maneuvering through rubble [17], and walking on uneven terrains [5]. These unique capabilities that only animals can perform have motivated many researchers to develop biomimetic robots.

The flapping-wing micro air vehicle (FWMAV) is one such type of robot. It has been popular for more than a decade [12,18–20]. With birds and insects as the main source of inspiration for designing robots, FWMAVs are typically categorized into bird- or insectinspired FWMAVs. Alternatively, they are differentiated into tailed or tailless FWMAVs based on the existence of the control surface on their tails, by which the differentiation agrees with the physical configurations of birds or insects, respectively [21–25]. In general, birds and insects have unique flight characteristics, which can be observed from the manner in which their wings are manipulated and modulated. Compared to birds, insects normally flap with higher flapping frequency and larger stroke amplitude to produce sufficient force to fly. The flapping stroke plane of insects is nearly horizontal, which causes similar lift force to be produced during a downstroke and an upstroke, whereas the flapping stroke plane of birds is practically vertical, which causes most of the lift force to be produced during downstroke [23,24]. These flight characteristics are mimicked by FWMAVs using various approaches and engineering designs, especially for producing wing flapping motion and other kinematics to generate force and moment [19,20,26].

Unlike tailed FWMAVs, tailless FWMAVs can exhibit steady hover and high maneuverability during a flight [27], indicating their potential applications in restricted spaces. To date, several tailless FWMAV designs have successfully performed stable flight, including the Nano Hummingbird [27], RoboBee [28], Delfly Nimble [29], KUBeetle [30–35], robotic hummingbird [36], Colibri [37], NUS-Robobird [38], and Purdue Hummingbird [39]. There are some similarities between each design. Nano Hummingbird, KUBeetle, robotic hummingbird, and Colibri use only one DC motor to generate thrust through one pair of wings and a few servos to manipulate the wing kinematics for control moment generation. They are mainly different in system integration and various transmission mechanisms, such as combination of string and crankshaft [27], combination of 4-bar linkage and pulleystring [30], and combination of slider-crank and 4-bar linkage [37]. In comparison, Delfly Nimble and NUS-Robobird use two DC motors to manipulate the flapping of two pairs of wings to generate thrust and roll moments. Servos are used to adjust the relative angle of two flapping mechanisms for pitch and yaw moment generation. A different control strategy was found in Purdue Hummingbird where each wing, one left and one right, is controlled by a DC motor to realize fully decoupled and independent wing kinematics manipulation. In this way, no servos are needed because control moments are generated by changing the wing flapping range and speed. While most tailless FWMAVs with successful flights use DC motors and servos, RoboBee uses piezoelectric actuators driven with high voltage. Stabilizing tailless FWMAVs is a significantly important objective in the development of tailless FWMAVs for flying. Owing to the inherent instability of tailless FWMAVs, feedback controllers are required to stabilize them during flight. In comparison, tailed FWMAVs, such as the Microbat [40], Delfly II [41], TL-Flowerfly [42], Konkuk University Biomimetic Ornithopter [43], H2Bird Ornithopter [44], and Robird [45], have an inherent stability [46] so that they can maintain equilibrium without any active manipulation of wing or tail kinematics. Recent findings on Hawkmoth [47] have indicated that the vibration of flapping wings potentially contributes to the improved stability of large insects and FWMAVs through a passive vibrational stabilization mechanism.

Owing to their feasibility and simplicity, proportional–integral–derivative (PID) feedback controllers and their simple variants are often implemented, particularly during the early development of tailless FWMAVs [19,30,37,38]. More complex variants of PID feedback controllers incorporating different strategies (such as neural network [48], particle swarm optimization [49], and dominant pole placement [50]) and autotuning algorithms (such as optimal approach [51,52] and fuzzy adaptive [53]) are being developed for various

applications. As with many biomimetic robots, tailless FWMAVs have less complicated body structures and control mechanisms than their biological counterparts. Thus, they can be considered as simplified versions of insects. However, owing to the nonlinear, time-varying, and periodic nature of these systems, describing the mechanics of tailless FWMAVs is still a considerable challenge. In addition, as tailless FWMAVs have different wing configurations [54–57], actuators, and control mechanisms [19,27–39], they cannot be fully generalized. Many efforts have been made to model the flapping kinematics, aerodynamics, and flight dynamics of flyers and FWMAVs. They are derived computationally or experimentally with various assumptions, ranging from time-invariant linear to nonlinear time-periodic models [58–60]. One of the main purposes of obtaining a dynamic model is to design a flight controller with various control techniques and schemes, such as PID [34,61], linear–quadratic–Gaussian [62], state feedback [63,64], nonlinear [65–67], sliding mode [68,69], robust [70,71], adaptive [72], adaptive neural network [73,74], and model-free [75]. While many dynamic models and controllers have been successfully developed, only a limited number have been implemented on tailless FWMAVs because some tailless FWMAVs cannot generate sufficient lift to perform free flight for implementation and verification.

In this study, the system identification of a tailless FWMAV called KUBeetle was conducted to improve the previously derived longitudinal mode extended dynamic model [34] and obtain a more accurate dynamic model to design a controller. The fabrication of tailless FWMAVs with successful flight makes it possible for system identification based on measured input and output data to be conducted [65,76–82]. This approach can potentially identify real dynamics and address the design changes of tailless FWMAVs to obtain a dynamic model with a high degree of similarity to actual systems. In a previous study [65], three single-axis closed-loop system identifications of the robotic hummingbird were conducted using a gimbal system and a Vicon motion tracking system with step reference command as the input. A second-order system was used to describe the pitch, roll, and yaw dynamic. A theoretical linear model for the piezo actuator positioning of Robobee wings' [76] was developed using an off-board displacement laser sensor. The identified model had an error of around 5%, but it was restricted to the wing driving system. A longitudinal gray box model system identification using open-loop system identification procedures was conducted on Delfly Nimble [77]. Data were obtained using the OptiTrack motion tracking system and the onboard system with doublet input as a reference command. The obtained linear time-averaged longitudinal dynamic models for hover and forward flight up to 1 m/s showed more than 80% fit.

The majority of system identification research uses high-speed camera systems with high data rates to obtain experimental flight data. High-speed camera systems are very expensive. In addition, the camera marker balls are relatively heavy compared to the small mass of the KUBeetle (only 17.7 g), which could affect the dynamics. Moreover, there exists little difference between the angle measured by a high-speed camera system and that obtained by an onboard sensor, due to misalignment, scale factor errors, and added measurement noises. So, the model obtained using the high-speed camera system should be modified to be used in a real feedback control system, while the dynamic model obtained using onboard sensors can be used directly without modification.

Many system identification methods have been developed. They are available for openloop and closed-loop system identification [78–82]. Because KUBeetle flight is unstable, a closed-loop system identification was chosen. Considering prior knowledge of the dynamics of the KUBeetle and the PD feedback controller, a gray box model approach with indirect methods [79,82] was utilized. The longitudinal mode extended dynamic model of KUBeetle, categorized as a linear, time-invariant, and single-input multipleoutput system with state space model structure, was refined with the help of system identification tools in MATLAB®. Several variations of input and output data pairing and different identification methods on model refinement were investigated and compared. By investigating dynamic model parameters and time-domain responses, the most appropriate

refined dynamic model was obtained. The main contributions of this paper are: (1) dynamic model parameters for the new KUBeetle using only low-cost onboard sensors outputs; and (2) stability and frequency domain analyses, which could be used for designing a closed-loop feedback controller.

The remainder of this paper is organized as follows. Section 2 describes the design and longitudinal mode dynamic model of the KUBeetle. Section 3 presents experimental flight, system identification setups, and the results. Finally, Section 4 provides conclusions and indicates towards future research topics.

#### **2. KUBeetle**

#### *2.1. KUBeetle Design*

Figure 1 shows the KUBeetle [30–35,83,84], a tailless FWMAV used for parameter identification in this study. The design of the KUBeetle in this study slightly differed from that of a previously developed KUBeetle [34]. They were mainly different in their flapping mechanisms and the location of the control board. One pair of wings of the KUBeetle were modulated by a DC motor using a rack and pinion mechanism to perform a high-amplitude flapping motion for thrust generation [85], whereas previous models used a pulley-string mechanism [30,31]. In the current KUBeetle, the control board was placed between the yaw servo motor and battery, which could lower the center of gravity (CG) position compared with the older version of the KUBeetle. The control board was designed to receive the reference command to control the FWMAV and transmit flight data in real time through a wireless communication RF module extension. A previous study [34] explains the detailed design and specifications of the control board. Figure 1b shows that the flight stability and control are maintained by manipulating the wing stroke plane using three servo motors to generate control moments for pitch, roll, and yaw.

**Figure 1.** A flight model of the KUBeetle: (**a**) detailed construction components used in the KUBeetle; (**b**) principle of the stroke-plane-change mechanism for control moment generation.

#### *2.2. Dynamic Model*

System identification using a gray box model requires an initial dynamic model with starting parameters. The longitudinal mode dynamic model of the KUBeetle used in this study was based on Newton–Euler equations of motion at near-hover conditions [86–89]. The dynamic model was extended to include filters and servo motor dynamics [34]. Parameters (specifically, the moments of inertia and filter parameters) were updated according to the latest KUBeetle design. Table 1 shows the stability and control derivatives obtained from the computational fluid dynamics (CFD) analysis and force and moment measurements. Figure 1 shows the body frame of the KUBeetle. The extended longitudinal mode dynamic model is described below [34,86–89]:

$$
\dot{\mathbf{x}}\_{\rm ext} = A\_{\rm ext} \,\mathbf{x}\_{\rm ext} + B\_{\rm ext} \,\boldsymbol{u}\_{\rm ext} \,\tag{1}
$$

$$y\_{\text{ext}} = \mathbb{C}\_{\text{ext}} \,\, \mathbf{x}\_{\text{ext}} \,\, \tag{2}$$

*xext* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *u w q ϑ q*ˆ*lpf ϑ*ˆ *lpf ϑ*ˆ *k f* ,*angle ϑ*ˆ *k f* ,*bias uservo* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (3)

where *u* and *w* are velocities along the *x*- and *z*-axis, respectively; *q* is the pitch rate (rad/s); *ϑ* is the pitch angle (rad); *uext* is the control input in the form of pitch servo angle; *yext* is the output or response of the KUBeetle; *q*ˆ*lpf* is the low-pass filtered pitch rate; *ϑ*ˆ *lpf* is the low-pass filtered pitch angle; *ϑ*ˆ *k f* ,*angle* and *ϑ*ˆ *k f* ,*bias* are the pitch angle and bias estimation of the Kalman filter, respectively, and *uservo* is the servo output. The system and input matrices of the extended dynamic model are given below [34,86–89]:

$$A\_{\rm ext} = \begin{bmatrix} \frac{X\_x}{m} & \frac{X\_x}{m} & \frac{X\_y}{m} & \frac{X\_y}{m} & 0 & 0 & 0 & 0 & 0 & \frac{X\_{\rm Tx}}{m} \\ \frac{Z\_x}{m} & \frac{Z\_y}{m} & \frac{Z\_y}{m} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{Z\_{\rm Tx}}{m} \\ \frac{M\_x}{I\_{\rm TF}} & \frac{M\_y}{I\_{\rm TF}} & \frac{M\_y}{I\_{\rm TF}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 204.3 & 0 & -204.3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 10.13 & 0 & -10.13 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0.329 & -0.329 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -0.05 & 0.05 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -30.67 \end{bmatrix}, \tag{4}$$

$$B\_{\rm ext} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}, \tag{5}$$

where *X* and *Z* are aerodynamic forces along the *x*- and *z*-axes, respectively, and *M* is the aerodynamic pitching moment about the CG. Subscripts *u, w*, and *q* indicate their aerodynamic derivatives with respect to *u*, *w*, and *q*. The subscript *γp* implies control input by the pitch servo. *<sup>m</sup>* is the mass of the KUBeetle (17.7 × <sup>10</sup>−<sup>3</sup> kg), *Iyy* is the pitching moment of inertia (12.43 <sup>×</sup>10−<sup>6</sup> kg m<sup>2</sup> ), and *g* is the gravitational acceleration 9.81 m s−<sup>2</sup> . The first four rows of the system matrix consisted of parameters of the basic dynamic model [34]. In the fifth and sixth rows, parameters of low-pass filters for pitch rate and pitch angle were included, respectively. The seventh and eighth rows included Kalman filter parameters for pitch angle and bias estimation. Finally, the last row contained servo motor dynamic parameters. The output matrix of the extended dynamic model, *Cext* (openloop) or *Cext*,*cl* (closed-loop), can be adjusted according to the number of available outputs

for system identification. For pitch angle and pitch rate outputs, the output matrix is given as follows:

$$\mathbb{C}\_{\text{ext}1} = \mathbb{C}\_{\text{ext},d1} = \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix} . \tag{6}$$

For pitch angle output only, the output matrix was given as:

$$\mathbb{C}\_{\text{ext}\mathfrak{L}} = \mathbb{C}\_{\text{ext}\mathfrak{L}} \mathfrak{L} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix} . \tag{7}$$

**Table 1.** Stability and control derivatives.


In general, it is not feasible or practical to identify an unstable system. Therefore, a simple PD feedback controller was added to the unstable system to obtain a stable closedloop control system. The same PD gains in a previous study [34] were used as the starting reference. They were then slightly tuned through experimental flight. Figure 2 shows the data acquisition structure of system identification with a PD feedback controller. Pitch rate estimation was obtained after the pitch rate output of the dynamic model was filtered by a low-pass filter [34]. Considering real time implementation, a first order low-pass filter (cutoff frequency 32 Hz) was used to reduce the high-frequency noise because of its simplicity which took a small calculation time and memory of the microprocessor. Pitch angle estimation was obtained by combining low-pass filtered pitch rate output and pitch angle estimated from acceleration outputs using the Kalman filter [34]. The KUBeetle showed stable flight performance with the current filters, especially on hovering, although there were some oscillating motions, indicating that the current setup was sufficient for an experimental flight.

**Figure 2.** Data acquisition structure of system identification with a PD feedback controller.

A closed-loop system identification is more challenging than an open-loop system identification due to existing feedback. The feedback can affect system identification performance by making the closed-loop system less sensitive to parameter changes in an open-loop system [79]. In addition, the natural dynamic of the KUBeetle toward the reference command might be regarded as a disturbance that needs to be damped and minimized by the controller [77,78]. However, the feedback also allows for more input in a specific frequency range without increasing the output magnitude [79]. Two different methods [79,82], direct and indirect, were tested to identify the longitudinal mode dynamic model of KUBeetle. According to Figure 2, the direct method used sets of control input and output data to obtain the open-loop dynamic model, while the indirect method used sets of reference command and output data to obtain the closed-loop dynamic model. It is possible to use the indirect approach, because the structure of the state space of the closed-loop extended dynamic model is fully known.

A simple PD controller with reference command is given bellow [34]:

$$
\hbar \mu\_{pd} = K\_P \left( \theta\_{ref} - \hat{\theta}\_{kf,angle} \right) - K\_D \theta\_{lpf,} \tag{8}
$$

where *upd* is the control input, *ϑref* is the angle reference command, *KP* is the proportional controller gain, and *KD* is the derivative controller gain. The system matrix and input matrix of the closed-loop extended dynamic model can be expressed as follows:

$$A\_{\rm ext,cl} = \begin{bmatrix} A\_{\rm ext} (1:8, 1:9) \\ 0 & 0 & 0 & -30.91 K\_D & 0 & -30.91 K\_P & 0 & -30.67 \end{bmatrix} \tag{9}$$

$$B\_{ext,cl} = \begin{bmatrix} B\_{ext}(1:8,1) \\ 30.91 K\_P \end{bmatrix}.\tag{10}$$

The addition of a PD controller to the extended dynamic model only changed the last row of the state space matrix where the control input was given to the servo. Either an openloop or a closed-loop extended dynamic model could be used as the starting reference in system identification with the gray box model approach. The open-loop extended dynamic model could be retrieved simply by setting controller gain parameters on the last row of the system matrix to zero and changing the last parameter of the input matrix to 30.91, according to Equations (4) and (5). Moreover, it should be noted that velocities along the xand z-axes were not measured. They could not be used for system identification, which makes the system identification more challenging.

Parameters of the filters, servo motor, and PD controller were fixed, because their values were already known and verified [34]. The gravity, integrator, and zero constants were also fixed. Thus, only nine stability derivatives and three control derivates were refinable in system identification. Each parameter had different physical characteristics and an effect on the dynamic of the KUBeetle [87–89]. Based on CFD analysis [89], *Xu*, *Zw*, and *Mu* should be negative and have relatively large magnitudes; *Xw*, *Zu*, *Zq*, and *Mw* should be close to zero; and *Xq* and *Mq* should have relatively small magnitudes. Overall, *Xu*, *Xq*, *Zw*, *Mu*, and *Mq*, the five parameters of stability derivatives, and *Xγp*, *Zγp*, and *Mγp*, the three parameters of control derivatives, were mainly focused on for system identification. During system identification, these eight parameters were refined with a relatively wide range of restriction. Conversely, non-focused parameters of stability derivatives were kept close to zero with a narrow range of restriction. Restrictions were adjusted considering computational fluid dynamics (CFD) results [87–89] and current KUBeetle design specifications to maintain the integrity of the identified dynamic model in relation to the physical system. While relaxing parameter restrictions might improve the result of system identification, some refined parameters might not match the real physical system.

#### **3. Experimental Flight and System Identification**

In this study, the system identification process was performed in the following steps: (1) experimental flight setup; (2) data acquisition for system identification via experimental flight; (3) linear model refinement in system identification using the gray box model approach; and (4) dynamic model analysis and verification.

#### *3.1. Experimental Flight Setup*

Figure 3 shows the experimental flight setup to obtain flight data for system identification. This setup consisted of a remote control, the KUBeetle, and a ground station. The remote control was used for throttle adjustment for taking off, hovering, and landing the KUBeetle. It was also used for adjusting and maintaining trim conditions during the experimental flight. The flight was conducted indoors without wind disturbance. All experimental flight data were obtained and processed with an onboard control board equipped with an inertial measurement unit (MPU-9250, InvenSense, San Jose, CA, USA). Flight data were then transmitted via wireless communication (NRF24L01+, Nordic Semiconductor, Trondheim, Norway) to the ground station. Considering that the bandwidth of the longitudinal mode dynamics with the PD controller was less than 20 rad/s, the sampling rate was set to be 100 Hz. Recorded data consisted of time, pitch reference command, pitch angle, and pitch rate. In addition, the control input was calculated from the reference command and output data with known control gains.

**Figure 3.** Experimental flight setup for data acquisition.

#### *3.2. Data Acqusition for System Identification via Experimental Flight*

Initially, the KUBeetle was manually kept in a hovering condition with minimum drift via a remote control. In hovering condition, the preprogrammed command in the control board was triggered through the remote control to automatically apply reference commands to the KUBeetle. Preprogrammed command was used to ensure the repeatability and consistency during each experimental flight. The reference command is a very important component of the system identification because it affects entire aspects of the system identification process. The reference command should be able to excite as many dynamic modes of interest and natural dynamics as possible to provide comprehensive data for system identification [77–80]. With the dynamic of the KUBeetle in consideration, two different reference commands were used in the experimental flight. The first reference command was in the form of a regular square signal with a relatively long and constant period (amplitude: ±17◦ (0.3 rad); period: 6 s) to analyze transient and steady state responses of the KUBeetle (Supplementary Materials Video S1). The second reference command consisted of an irregular square signal with varying periods (amplitude: ±17◦ (0.3 rad); period: 0.2 to 1.4 s) to cover the frequency bandwidth of interest (Supplementary Materials Video S2). The second reference command was generated using a pseudo-random binary signal, which was a periodic and deterministic signal with properties such as white noise [79]. The amplitude and period of the reference command were chosen considering the operating range of the KUBeetle. Assuming the same amplitude and frequency range, the square signal covered a similar range to that of a chirp signal, as indicated by the frequency spectrum of reference commands shown in Figure 4:

For clarity, the extended dynamic model with parameters obtained from CFD and force and moment measurements, later used as the initial dynamic model in system identification, is called the CFD dynamic model. Figure 5 shows experimental and simulation results of the KUBeetle with a regular square command for pitch. Together with the preprogrammed command, the constant trim control input from the remote control was applied to the KUBeetle. The constant trim command was subtracted to obtain the reference command from the actual command input applied. Small noises were added while the trim control input was being read from the pulse-width-modulated signal. The reference command was the summation of the preprogrammed command and the added trim command noise. Figure 5a shows that the closed-loop KUBeetle has a relatively slow response, high overshoot, and large steady state error. After the reference command was provided, the KUBeetle attempted to follow the reference command in the beginning by responding quickly according to the stroke plane change via the servo motor. However, it was gradually stabilized to a low pitching angle, which maintained the large steady state error. This phenomenon was mainly caused by the position of the CG that affected the force and the moment generated by the gravity force during flight. The CG was located at a relatively large distance below the aerodynamic center (AC). Figure 5b presents the pitch rate responses to a regular square reference command, showing that the measured pitch rate still had a large periodic and high-frequency noise due to wing flapping.

Figure 6 demonstrates experimental and simulation results of the KUBeetle with an irregular square reference command. Figure 6a shows the pitch angle responses. In some intervals, the pitch angle response misleadingly looks as if it is following the command signal well, although some phase delays exist. This is due to overshooting responses at the beginning of step inputs. However, as the frequency of the reference command increased, it could not follow the reference command. Figure 6b illustrates that the pitch rate was contaminated with large noise, as in the regular square reference command case shown in Figure 5b. The high-frequency noise was still very large on the angular rate signal, although it was low-pass filtered (cutoff frequency 32 Hz). However, it was not easy to remove it without invoking much phase delay, which degraded the closed-loop stability. Thus, the current setup was kept for the experimental flight.

**Figure 4.** Frequency spectrum of reference commands.

**Figure 5.** Experimental flight data and simulation results for system identification with regular square reference inputs: (**a**) pitch angle responses, and (**b**) pitch rate responses.

**Figure 6.** Experimental flight data and simulation results for system identification with irregular square reference inputs: (**a**) pitch angle responses, and (**b**) pitch rate responses.

#### *3.3. Linear Model Refinement in System Identification*

System identification processes were performed with the help of system identification tools in MATLAB®. The dynamic model generated from the system identification using the black box model approach was found to be very complex. It could not be matched with real physical parameters of the KUBeetle. Thus, this research only focused on the system identification using the gray box model approach. The gray box model approach with direct and indirect methods was tested to refine open- and closed-loop dynamic models. Two different pairings were considered for the direct and indirect methods, respectively. Firstly, the control input was given as input, while the pitch rate and pitch angle were given as outputs. Secondly, the reference command was given as input, while the pitch rate and pitch angle were given as outputs. Equations (4) and (5) present the structure of the system and input matrices of the dynamic model as fixed for the open-loop dynamic model, while Equations (9) and (10) present those for the closed-loop dynamic model. In addition, given values of parameters in Equations (4)–(10), obtained from computations and measurements, were used as the initial values for system identification. The pitch angle output was given higher weighting, because it had smaller noise compared to the pitch rate output which had larger high-frequency noise that could affect the accuracy of the identified dynamic model. According to the experimental flight setup and data, the obtained dynamic model mainly represented the dynamic of KUBeetle during hover and forward flights with pitch angle under 0.3 radian. Adaptive Gauss–Newton (GNA), Levenberg–Marquardt (LM), and gradient search (GRAD) methods were utilized for model refinements. Table 2 summarizes the general characteristics of each identification method [90–92]:

The total iteration for each search method was limited to 2000 iterations, as further iterations provided negligible improvement. Results from each search method were compared and evaluated. One of the indicators used was the fit percentage of the responses of the identified dynamic model to the experimental flight data. It was defined as:

$$Fit\,Percentage = 100 \left( 1 - \frac{\sqrt{\sum\_{i=1}^{N} \left( y\_{measured} - y\_{model} \right)^2}}{\sqrt{\sum\_{i=1}^{N} \left( y\_{measured} - \overline{y\_{measured}} \right)^2}} \right), \tag{11}$$

where, *N* was the number of samples, *ymeasured* was the experimental flight data measured, *ymodel* was the simulated response of the dynamic model identified, and *ymeasured* was the mean of experimental flight data measured.


**Table 2.** General characteristics of each identification method.

Firstly, system identification using the gray box model approach with the direct method was conducted. Closed-loop system identification with the direct method is often used, because it works regardless of the complexity of the regulator. In addition, it does not require knowledge of the characteristic of the feedback [79]. However, the system identification result depended on the quality of the disturbance model, which should also be identified during model refinement. The open-loop CFD dynamic model in Equations (4) and (5) was used as the initial dynamic model. However, the direct method was found to be unsuitable for the system identification of the KUBeetle using the current experimental flight data as the disturbance model could not be identified well, which caused responses of the identified model not to match experimental flight data. Note that some working identification methods for open-loop system identification may fail when they are applied to a closed-loop system identification [79].

Secondly, considering that the open-loop dynamic model could easily be retrieved from the closed-loop dynamic model, system identification using the gray box model approach with the indirect method was conducted. In the indirect method, the system matrix of the closed-loop CFD dynamic model in Equations (9) and (10) was used as the initial dynamic model. Through system identification with several different methods and various restriction settings, several identified dynamic models were obtained. The criteria to select the identified dynamic model were: (1) the fit percentage should be high enough; (2) time-domain responses of the closed-loop dynamic model should have similar characteristics to the experimental flight data; and (3) the open-loop dynamic model should have similar stable and unstable poles and zeros to those of the computationally calculated dynamic model, whereas the closed-loop dynamic model with the PD feedback controller should be stable.

#### *3.4. Dynamic Model Analysis and Verification*

Single or a combination of multiple sets of experimental flight data were tested to obtain the best identified dynamic model. This was completed to evaluate the consistency of the identification results and sensitivity of the identification methods. It was found that identified parameters related to velocities along the x- and z-axes had relatively larger deviations because they were not available using onboard sensors. In addition, high-frequency noises and nonlinear body dynamics not included in the model also affected identification results. By following all the selection criteria, six identified dynamic models were considered. Table 3 shows parameters of the identified dynamic model with the indirect method. Each of the selected dynamic models was named based on the identification method used to refine it and numbers represent datasets used for the identification. Fit percentages of

identified models obtained using irregular and high-frequency reference signals were much lower than those obtained using regular and low-frequency reference signals. Thus, they were discarded. The fit percentage of the dynamic model from CFD was only about 22%. This low fit percentage was mainly because only wing parts were considered to obtain the stability derivative-related parameters, instead of using the whole body of the KUBeetle. The whole body of the KUBeetle, including the wings, was considered to obtain mass and moments of inertia. The average fit percentage of the identified dynamic models was approximately 49%, which was more than double that of the CFD dynamic model. Most parameter values obtained from the three different identification methods and CFD were quite similar to each other. The CG of the new KUBeetle design, being lower than the CG of the previous KUBeetle design, increased the magnitude of *Mu*. This was in accordance with the shift of the CG position due to the relocation of the control board and RF module from the top of the flapping mechanism to under the motor. The increase in *Mq* showed that the system had higher damping than that initially given from CFD. This indicated that the body of the KUBeetle greatly affected system stability. According to system identification results, aerodynamic forces along the x- and z-axes, *Xγ<sup>p</sup>* and *Zγp*, had larger magnitudes than that previously measured. When measuring the control derivative, only the wings and frames near the wing were included. Because of a mounting problem, the lower body of the KUBeetle could not be included, which could be one of the possible reasons for the difference. Other parameters maintained similar values to CFD values.


**Table 3.** Refined parameters of identified dynamic models using an indirect method.

Using the current flight data and dynamic model structure, it was not easy to obtain dynamic model parameters with high fit percentage, even using the black box model approach where there were no parameter restrictions. Moreover, it was found that the parameters obtained were a little sensitive to the datasets used for system identification. The possible reasons for this were: (1) velocities along x- and z-axes were not measured, which degraded the accuracy of the related parameters, *Xu*, *Xw*, *Xq*, *Zu*, *Zw*, and *Zq*; (2) highfrequency noise also contributed to accuracy degradation; and (3) there were nonlinear coupling dynamic terms not included in the model. It has been reported that accuracy can be improved if velocities along the x- and z-axes are available and nonlinear coupling terms are included [77].

Although all dynamic model candidates showed similar fit percentages, GRAD1 was chosen for analysis and evaluation since it had the highest fit percentage. Identified parameters were evaluated in the basic form of the dynamic model [34], where filters and servo dynamics were not included. The basic dynamic model is described as follows:

$$
\dot{\mathbf{x}}\_{basic} = A\_{basic\ \mathbf{x}} \mathbf{x}\_{basic} + B\_{basic\ \mathbf{u}} \mathbf{u}\_{basic\ \mathbf{w}} \tag{12}
$$

$$y\_{\text{basic}} = \mathbb{C}\_{\text{basic}} \, \mathbf{x}\_{\text{basic}}.\tag{13}$$

$$\begin{aligned} \left. \begin{array}{c} \boldsymbol{w} \\ \boldsymbol{w}\_{\text{basic}} = \begin{bmatrix} \boldsymbol{u} \\ \boldsymbol{w} \\ \boldsymbol{q} \\ \boldsymbol{\theta} \end{bmatrix} \right\} \end{aligned} \tag{14} $$

$$A\_{\text{basic}} = \begin{bmatrix} \frac{\chi\_{\text{u}}}{m} & \frac{\chi\_{\text{u}}}{m} & \frac{\chi\_{\text{d}}}{m} & \mathcal{S} \\ \frac{Z\_{\text{u}}}{m} & \frac{Z\_{\text{w}}}{m} & \frac{Z\_{\text{q}}}{m} & 0 \\ \frac{M\_{\text{w}}}{T\_{\text{yy}}} & \frac{M\_{\text{w}}}{T\_{\text{yy}}} & \frac{M\_{\text{q}}}{T\_{\text{yy}}} & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix},\tag{15}$$

$$B\_{\rm brsic} = \begin{bmatrix} \frac{\chi\_{\gamma\_F}}{\bar{Z}\_{\gamma\_F}^{\rm yn}}\\ \frac{\bar{Z}\_{\gamma\_F}^{\rm yn}}{\bar{I}\_{\rm yy}}\\ 0 \end{bmatrix},\tag{16}$$

and *Cbasic* was set as <sup>0</sup> <sup>001</sup> to describe the pitch angle measurement.

Figure 7 shows the open-loop frequency responses and a Nyquist diagram of the GRAD1 and CFD dynamic models. When the frequency was under 0.01 rad/s, the lowfrequency gain of the identified GRAD1 dynamic model was similar to that of the CFD dynamic model. As the frequency increased to higher than 0.3 rad/s, the gain of the GRAD1 dynamic model was smaller than that of the CFD dynamic model. Moreover, at mid-frequency, the GRAD1 dynamic model had a local minimum point. As the frequency further increased to higher than 5 rad/s, both dynamic models had similar characteristics. Poles of the GRAD1 were (2.86 ± 8.46i, −9.35, −1.64), while those of the CFD model were (3.69 ± 7.58i, −9.22, −0.76). The unstable pole locations did not move much. However, the magnitude of the stable dominant pole of the GRAD1 dynamic model was larger than that of the CFD dynamic model. Thus, it is expected that response of the GRAD1 dynamic model is faster than that of the CFD dynamic model. Both basic dynamic models encircle the (−1, 0) point two times in the Nyquist diagram, meaning that if the loop is closed, the two systems will be stable. However, Figure 7b shows that both models have very small phase margins, with less than 8 degrees.

The simulation was used to analyze and compare extended dynamic models as shown in Figure 8. Gain and phase margins of the GRAD1 dynamic model were 17.2 dB and 39.0◦, respectively, while the gain and phase margins of the CFD dynamic model were 16.5 dB and 30.3◦, respectively. In the Nyquist diagram, both extended dynamic models encircled the (−1, 0) point two times. Both dynamic models had small gains in the low-frequency range, which resulted in a large steady state error. To reduce the steady state error, the loop gain should be increased. However, this may decrease the stability margins.

Figure 9 shows experimental flight data and simulation results of the GRAD1 dynamic model using the PD feedback controller. As shown in Figure 9a, the pitch angle response of the GRAD1 dynamic model was closer to the pitch angle from the experimental flight data compared to that of the dynamic model from CFD. The peak and valley points predicted using the GRAD1 dynamic model were more accurate than those of the CFD dynamic model. Figure 9b also shows that the GRAD1 dynamic model was better than the CFD dynamic model in following the pitch rate from the experimental flight data, neglecting the high-frequency noise. The high-frequency noise was due to wing flapping, which degraded the performances of system identification and control.

**Figure 7.** Frequency responses (basic dynamic models only): (**a**) Bode plot; (**b**) Nyquist diagram.

**Figure 8.** Frequency responses of extended dynamic models with a PD controller: (**a**) Bode plot; (**b**) Nyquist diagram.

**Figure 9.** Experimental flight data and simulation results of the GRAD1 dynamic model using a PD feedback controller: (**a**) pitch angle, and (**b**) pitch rate for the provided reference command of the square signal with long period and constant frequency.

#### **4. Conclusions**

In this study, we identified a dynamic model of the KUBeetle with several identification approaches and methodologies using angle and angular rate data from onboard sensors only. The gray box model approach was preferred to the black box model approach, because the structure and integrity of parameters of the dynamic model could be maintained. Higher weighting was given to the angle output because it had much lower noises compared to angular rate output. It was found that dynamic models with high fit percentage were obtained using datasets which included the data obtained using regular and low-frequency reference commands. In addition, obtained parameters were slightly sensitive to the quality of the dataset. GNA, LM, and GRAD identification methods showed similar results to each other for two different datasets, while GRAD1 showed the highest fit percentage. These consistent results indicate that the dynamic model was quite reliable and could represent the main dynamic of the KUBeetle, although nonlinear dynamics were not fully identified. Compared to that of the CFD dynamic model, the fit percentage was improved from 22% to 49%. Increases in *Mu* complied with the lower CG of the KUBeetle. In addition, the increases in *Mq* indicated higher damping caused by the body dynamics of the KUBeetle that were not considered previously. The obtained dynamic model mainly works for predicting hover and forward flight responses with pitch angle reference command under 0.3 radian, which is used in an experimental flight. The obtained dynamic model described the pitch angular motion quite well, while the linear velocity-related parameters needed more improvement. This is because there was no onboard speed sensor.

From frequency domain analysis, it was found that the open-loop GRAD1 model had a small phase margin indicating that KUBeetle stability could be easily affected by the angular rate feedback gain changes. Moreover, the GRAD1 dynamic model had a faster response because its stable dominant pole had a larger magnitude than the CFD dynamic model. In addition, it showed a similar response to the experimental flight data, neglecting high-frequency noises. It is expected that if the high-frequency noise due to flapping could be reduced and the velocity could be estimated, the fit percentage could be further improved. The time-domain and frequency responses of the identified dynamic model were considerably closer to experimental flight data than those obtained by numerical computation. Further studies are needed to reduce the high-frequency noise due to flapping, to better estimate forward and vertical velocities, and to improve the accuracy of identification.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/app12052486/s1, Video S1: Experimental Flight PD controller with regular frequency command (speed 0.5×); Video S2: Experimental Flight PD controller with irregular frequency command (speed 0.5×).

**Author Contributions:** Conceptualization, T.-S.K.; Data curation, S.A. and G.-H.H.; Formal analysis, S.A.; Investigation, S.A. and G.-H.H.; Resources, H.-C.P.; Supervision, T.-S.K.; Validation, S.A. and H.-C.P.; Writing—original draft, S.A.; Writing—review & editing, H.-C.P. and T.-S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by a grant (No. 2020R1F1A1077026) of the National Research Foundation (NRF) funded by MSIT and a grant (No. 2020R1A6A1A03046811) from the Ministry of Education, Korea.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The authors also wish to express their gratitude to Hoang Vu Phan, Thi Kim Loan Au, and Khanh Nguyen for their support in this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

