Adv Syst Sci Appl 2023; 03:1–15 Published online at https://ijassa.ipu.ru. Frequency Domain Identification of the Quadcopter Attitude Dynamics Vadim Alexandrov*, Ilya Rezkov, Dmitrii Shatov, Yuri Morozov V.A. Trapeznikov Institute of Control Sciences, Russian Academy of Sciences, Moscow, Russia Abstract: The paper presents the identification of a linearized model of the quadcopter attitude dynamics. The attitude is described by the Euler angles of roll, pitch and yaw. Closed loop identification is performed, when the quadcopter control system operates providing flight stability. The experimental flight data, when the sequence of the sine waves is fed as required value for each angle separately, are used. The transfer functions from the controls, which provide the torques related to the body frame axes, to the Euler angles are obtained via the finite frequency identification procedure. Moreover, components of the rotational dynamics are considered in detail. The unknown parameters of the transfer functions are found by optimization procedure using the experimentally obtained values of the frequency response for the set of test frequencies. The difference in the parameters values of the fixed linearized model structure, identified for two operating points, shows that the tilt angle dynamics is sufficiently nonlinear. Keywords: system identification, frequency domain identification, quadcopter attitude, Euler angles, transfer function 1. INTRODUCTION The complex flight dynamics of a quadcopter is a widely studied benchmark for control design methods. Although there are several well-functioning control systems for commercial quadcopters and there are open source projects for experimental multirotor setups, research is ongoing on the application of classical and original methods of control design for the quadcopter flight control system. Identification is an essential step for model based control design. The quadcopter model identification was studied by various methods [14]. In [7, 10, 13] coefficients of the linearized quadcopter models were identified using the software CIFER. The test signal called frequency sweep is applied during the test flight of the quadcopter, then the corresponding frequency response is calculated and model coefficients are found. In [13] the full linearized quadcopter model in form of the several transfer functions was identified, in particular, for lateral and longitudinal motions there were obtained the third order systems with pure differentiation plus time delay and the poles were stable, one was real and the rest were complex-conjugate. The authors of [7, 10] identified quadcopter models in state-space form basing on the results of transfer function identification. In [7] the third order functions similar to [13] were used but augmented with a real zero. The article [8] develops a method of linear model identification in state-space form for a multirotor vehicle based on the equation-error maximum likelihood estimator which reduces to least-square approach applied to experiment data in the time domain. The real-time identification method for a quadcopter model was proposed in [1]. Here, the recursive Fourier transform coupled with least-square estimation was used to find the state-space model coefficients estimates. The proposed method uses multisine excitation signal added to the full motor command. ∗Corresponding author: va.alexandrov@yandex.ru 2 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV The finite frequency identification approach [2] is used in our works to find the transfer function of the control plant model. The problem for altitude control loop identification has been studied in [5]. The transfer function from the horizontal projection of a full thrust to the horizontal velocity of the quadcopter was experimentally identified, and aerodynamic drag was detected in [4]. It should be noted, that the horizontal part of thrust is determined by the combination of pitch and roll angles. Identification of the tilt angle dynamics is considered in [6], where, in addition to frequency domain identification, the time domain nonlinear gray- box parameters estimation is used. This article develops the results of [6], considering also the rotational dynamics for the yaw angle. The problem of transfer functions identification for separated loops of roll, pitch and yaw angles that describe the quadcopter attitude is stated in Section 2 more detailed. Section 3 briefly describes the finite frequency identification method and results of the experimental flight data processing. The structures of the transfer functions for roll, pitch and yaw loops are introduced and evaluated via optimization procedure in Section 4. Moreover, the test flight data for another operating point of the pitch angle is analyzed. Conclusion is in Section 5. The experimental quadcopter [4–6] is the frame S500 with 10-inch propellers, the on- board flight controller Pixhawk 4 and the standard ArduPilot software suite. 2. PROBLEM STATEMENT The quadcopter rotational motion description based on Newton-Euler formalism [14] is ϕ̇ θ̇ ψ̇  = [ 1 sinϕ tan θ cosϕ tan θ 0 cosϕ − sinϕ 0 sinϕ/ cos θ cosϕ/ cos θ ][ p q r ] , (2.1) [ ṗ q̇ ṙ ] =  Jy − Jz Jx qr Jz − Jx Jy pr Jx − Jy Jz pq +  1 Jx τϕ 1 Jy τθ 1 Jz τψ  , (2.2) where ϕ, θ, ψ are roll, pitch, yaw Euler angles, p, q, r are angular velocities in body frame coordinates x, y, z respectively, Jx, Jy, Jz are inertia values assuming that the quadcopter is symmetric about the body axes, τϕ, τθ, τψ are torques around the respective axes, which are shown on Fig. 2.1. The torques are determined by the difference in the rotational speeds Ω1...4 of each propeller. Formulas for the case of symmetric X-type quadcopter frame are τϕ = Kτ (Ω2 2 + Ω3 2 − Ω1 2 − Ω4 2), τθ = Kτ (Ω2 2 + Ω4 2 − Ω1 2 − Ω3 2), τψ = Kψ(Ω1 2 + Ω2 2 − Ω3 2 − Ω4 2), (2.3) where Kτ , Kψ are proportionality coefficients. Thus, the control signals u1...4 for each motor are formed as u1 = uF − 0.5uϕ − 0.5uθ + uψ, u2 = uF + 0.5uϕ + 0.5uθ + uψ, u3 = uF + 0.5uϕ − 0.5uθ − uψ, u4 = uF − 0.5uϕ + 0.5uθ − uψ, Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 3 Fig. 2.1. Body frame axes where uF is the control value for thrust providing the altitude hold, uϕ, uθ, uψ are controls for torques related to body frame axes x, y and z. It is assumed that u1...4 ∈ [0, 1], uϕ, uθ, uψ ∈ [−0.5, 0.5]. Some dynamics of motors is present from ui to Ωi. Usually installed electronic speed controllers (ESC) for each motor provide a constant gain from control ui to the rotational speed Ωi. A simplified description for the rotor speed dynamics is a first order transfer function [9]: PΩ(s) = KΩ TΩs+ 1 (2.4) with the time constant TΩ and the gainKΩ. A time delay Tdelay should be added to this transfer function to account the ESC discrete operation. If a separate control for a one channel is considered, when, for example, uϕ = const ̸= 0 and uθ = uψ = 0, then the rotational speed can be divided into parts ΩF and Ωϕ dependent from the controls uF and uϕ Ωi = ΩF ± Ωϕ and then Ωi 2 = ΩF 2 ± 2ΩFΩϕ + Ωϕ 2. Therefore, squares are canceled in (2.3), and the torque is linear with respect to the additional rotational speed generated from the angular control uϕ, when ΩF is assumed as a constant: τϕ = 4KτΩFΩϕ. Moreover, when uθ = uψ = 0 for roll control consideration separately, then one can assumed that angular velocities q = r = 0 in (2.2). Thus, a basic model for roll angular velocity p from the angular control uϕ can be described by transfer function Pp(s) = Kϕ s(TΩs+ 1) e−sTdelay , (2.5) Kϕ = 4KτΩFKΩ Jx , where values of gain Kϕ, time constant TΩ and time delay Tdelay should be determined. The transfer function for roll angle ϕ can be obtained from (2.5) by adding one more integrator if Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 4 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV pitch and yaw angles values are negligible in (2.1): Pϕ(s) = Kϕ s2(TΩs+ 1) e−sTdelay . (2.6) The basic models for pitch and yaw angles loops are similar to (2.6) on the strength of (2.2). Nevertheless, the experimentally identified models in a number of papers have a different structure. Transfer functions for roll and yaw angular velocities obtained in [13] are Pp(s) = K1s (T1s+ 1)(s2 + 2ξωps+ ωp2) e−sTdelay , (2.7) Pr(s) = K2 T2s+ 1 e−sTdelay , (2.8) and the transfer function for the pitch angular velocity is similar to the roll one. Another model structure for roll and pitch angular velocities is considered in [12], where the state- space model for the roll loop is[ v̇ ṗ ϕ̇ ] = [ a1 0 g a2 0 0 0 1 0 ][ v p ϕ ] + [ b1 b2 0 ] uϕ(t− Tdelay), (2.9) where a1, a2, b1, b2 are identified parameters, g is the gravitational acceleration, v is the translational velocity projection onto the body frame y axis. The transfer function from control uϕ to angular velocity p for the model (2.9) is Pp(s) = K̃s(s+ b̃1) s3 + ã2s2 + ã0 e−sTdelay , where K̃, b̃1, ã2, ã0, are some parameters, calculated from parameters of the system (2.9). The pitch loop model in the paper [12] is similar to (2.9) and the yaw loop model is the same as (2.8). The first order rotor speed dynamics is added to pitch and roll model of the form (2.9) in [7]. Moreover, the yaw loop model is complemented by a lead time constant Tlead, i.e. the transfer function for yaw angular velocity is Pr(s) = K2(Tleads+ 1) (T2s+ 1)(TΩs+ 1) e−sTdelay . (2.10) One more example of the pitch and roll loop model is described in [10]. For brevity, the roll model only is separated:[ v̇ ṗ ϕ̇ ] = [ a11 a12 g a21 a22 0 0 1 0 ][ v p ϕ ] + [ 0 b2 0 ] uϕ, (2.11) where a11, a12, a21, a22, b2 are identified parameters. The corresponding transfer function is Pp(s) = K̃s(s+ b̃1) s3 + ã2s2 + ã1s+ ã0 , (2.12) where K̃, b̃1, ã2, ã1, ã0 are coefficients calculated from the parameters of (2.11). Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 5 Thus, the model structure (2.6) is not confirmed by the identification results. However, there is no generally accepted model, parameters of which need to be identified. Therefore, the considered problem is to find the transfer functions from the controls uϕ, uθ, uψ to roll, pitch and yaw angles in the form P (s) = b̂ms m + . . .+ b̂0 sn + ân−1sn−1 + . . .+ â0 (2.13) with orders m,n under the condition m < n and parameters b̂0, . . . , b̂m, â0, . . . , ân−1. 3. FINITE FREQUENCY IDENTIFICATION The finite-frequency identification approach is based on analysis of the experimental data when the sine wave test signals on several different frequencies from wide enough range are fed to the system input. The quadcopter cannot be tested in open loop mode, since horizontal stabilization is provided by default attitude controllers. In the case of closed loop identification [3,5] the sine wave test signals are added to the controller setpoint. The required values of the pitch, roll and yaw angles θ∗, ϕ∗, ψ∗ for attitude controllers in identification test for the pitch angle loop should be the follows ones: θ∗ = ηi sin(ωi(t− t0)), ϕ∗ = 0, ψ∗ = 0, (3.14) where ηi and ωi are amplitude and frequency of the test signal (i = 0, .., nω, where nω is a number of the test frequencies from a given test set that are fed consequently to the system input), t is time and t0 is the start time of the experiment on the current frequency. Then, uθ, uϕ and uψ are outputs of corresponding controllers to stabilize the roll and yaw angles values and provide the sinusoidal variation of the pitch angle. Identification tests for the roll and yaw angles loops are similar. The ArduPilot onboard control system is used for the test flight, where the author’s new software procedures are added to form the test signals (3.14) and collect the experimental data. The finite-frequency identification procedure consists of two steps: 1) obtaining estimates of frequency response for several test frequencies ωi and 2) finding estimates of coefficients of the identified transfer function using the frequency response estimates. The frequency response estimates are obtained according to P (jωi) = αi + jβi = αyi + jβyi αui + jβui , i = 1, . . . , nω, (3.15) where αi and βi are called the frequency parameters, which are found using output αyi, jβyi and input αui, jβui values calculated from the experimental data by Fourier filter formulae [2, 5]: αyi = 2 ηiTi ∫ t0+Ti t0 θ(t) sin(ωi(t− t0))dt, βyi = 2 ηiTi ∫ t0+Ti t0 θ(t) cos(ωi(t− t0))dt, αui = 2 ηiTi ∫ t0+Ti t0 uθ(t) sin(ωi(t− t0))dt, βui = 2 ηiTi ∫ t0+Ti t0 uθ(t) cos(ωi(t− t0))dt, i = 1, . . . , nω, (3.16) Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 6 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV where Ti = 2πq/ωi is the filtration time (experiment duration), where q ∈ N is the number of periods of the test frequency. Some examples of the experimental data are shown on Fig. 3.2. It is seen that the system output value follows the setpoint sine wave very accurately for low frequency but with lag and reduced amplitude for high frequency, which is in line with expectations for the closed loop system. In other hand, the control signal absorbs the noises and nonlinearities. The frequency parameters estimates obtained from the experimental flight data, where inputs are the controls uϕ,θ,ψ ∈ [−0.5, 0.5] and outputs are angles ϕ, θ, ψ in radians, are in table 3.1. Fig. 3.2. Sine wave test for the pitch closed loop with frequencies 16 and 2 rad/s Table 3.1. The frequency parameters estimates Roll Pitch Yaw i Frequency, ωi, rad/s αi βi αi βi αi βi 1 30.28 -0.044 -0.008 2 20.11 -0.397 0.344 -0.385 0.287 -0.094 -0.055 3 16.01 -0.736 0.552 -0.724 0.438 -0.139 -0.0864 4 12.03 -1.532 0.722 -1.409 0.628 -0.195 -0.122 5 8.00 -3.642 1.517 -3.571 0.972 -0.316 -0.243 6 5.34 -8.642 3.662 -6.815 2.142 -0.630 -0.422 7 4.00 -9.436 9.005 -9.040 5.771 -1.039 -0.552 8 3.50 -8.479 12.253 -7.689 6.591 -1.317 -0.731 9 3.00 -4.108 11.115 -6.914 7.417 -1.657 -0.838 10 2.00 0.559 8.897 -2.384 7.335 -3.217 -2.439 11 1.00 0.843 3.979 0.007 3.643 These estimates are used to write a system of linear equations, where they are equated for each test frequency to (2.13) with unknown parameters, but with chosen order values n,m : b̂m(jωi) m + . . .+ b̂0 (jωi)n + ân−1(jωi)n−1 + . . .+ â0 = αi + jβi, (3.17) where ωi, αi, βi are values from experiments and b̂0, . . . , b̂m, â0, . . . , ân−1 are unknowns. When condition 2nω > n+m+ 1 is satisfied, the overdetermined system of 2nω linear algebraic equations with n+m+ 1 unknowns is obtained by separating the real and imaginary parts in (3.17). Then, the solution of this system provides estimates of the transfer function coefficients by the least squares criterion. The values n,m are chosen by decreasing until the identified transfer function contains no obviously cancellable stable poles and zeros. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 7 The identification results for the data from table 3.1 are Pϕ(s) = 58.65(s+ 0.304) (0.050s+ 1)(0.193s+ 1)(s2 − 3.44s+ 13.8) , (3.18) Pθ(s) = 49.39(s− 0.156) (0.046s+ 1)(0.181s+ 1)(s2 − 3.79s+ 12.9) , (3.19) Pψ(s) = 1294(s+ 10.93) s(0.834s+ 1)(s2 + 28.9s+ 1036) . (3.20) Comparison of these functions with the data of table 3.1 is shown on Fig. 3.3 – 3.5. Fig. 3.3. Nyquist diagram of identified transfer function for roll Fig. 3.4. Nyquist diagram of identified transfer function for pitch Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 8 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV Fig. 3.5. Nyquist diagram of identified transfer function for yaw The identified transfer functions differ from those given in Section 2. The function structures for roll and pitch are expectedly the same. This structure is similar to the form of (2.12), where the differentiating element is added in the numerator, since (2.12) is the transfer function for angular velocity. The smallest time constant in the denominators of (3.18) and (3.19) can be considered a first-order rotor speed model as proposed in [7] added to (2.12). More detailed study of the experimental data is presented in the next section. 4. IDENTIFICATION IN DETAILS 4.1. Rotors dynamics Electronic speed controllers (ESC) type BLHeli32 are installed on the experimental quadcopter, so the rotors angular velocities Ω1...4 are measured. Then, the rotor speed dynamics can be identified from the collected experimental data. The frequency parameters obtained via (3.15)–(3.16) during the yaw tests from the controls u1...4 ∈ [0, 1] to the rotors angular velocities Ω1...4 in rad/s are in the table 4.2 for each rotor. Table 4.2. The frequency parameters estimates for the rotors dynamics rotor 1 rotor 2 rotor 3 rotor 4 i Frequency, ωi, rad/s αi βi αi βi αi βi αi βi 1 30.28 251 -570 285 -443 273 -508 221 -589 2 20.11 475 -568 538 -521 518 -560 390 -591 3 16.01 627 -551 665 -517 575 -557 620 -516 4 12.03 751 -503 772 -457 754 -490 737 -478 5 8.00 899 -378 908 -341 870 -334 907 -386 6 5.34 991 -281 936 -254 995 -248 991 -257 7 4.00 997 -211 990 -202 996 -203 980 -199 8 3.50 1012 -189 995 -147 1003 -178 1008 -175 9 3.00 1013 -145 1015 -125 -1002 -158 1009 -160 The structure of the main dynamics for rotors can be set as the first order (2.4) with delay. The delay is due to the discreteness of the flight controller and ESC. Thus, the parameters Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 9 KΩ, TΩ, Tdelay of the transfer function should be identified: PΩ(s) = KΩ TΩs+ 1 e−sTdelay . (4.21) The estimates of (4.21) parameters can be found via optimization procedure with variables vector x = [KΩ, TΩ, Tdelay]. For the convenience of notation, hereinafter the vector of variables for the optimization procedure will be denoted as a row vector. The objective function is min x nω∑ i=1 4∑ k=1 |PΩ(x, jωi)− (αki + jβki)| |αki + jβki| , where PΩ(x, jωi) is the value for s = jωi of the transfer function (4.21) with parameters from the variables vector x, k is the rotor number and αki, βki are corresponding estimates of the frequency parameters from the table 4.2, i.e., an average function for four rotors is found. The identification result is PΩ(s) = 1033 0.048s+ 1 e−0.0042s. (4.22) The delay value correlates with the fast loop of the flight controller, which operates with step of 0.0025 s. The rotors time constant is expectedly small, but it cannot be considered negligible when the controller sampling time is essentially less. Its Nyquist diagram and experimental points from the table 4.2 are shown on Fig. 4.6. Thus, one part of all considered loops is identified. Fig. 4.6. Nyquist diagram of identified transfer function for rotors 4.2. Roll and pitch model Further, the pitch model will be considered mainly, since the roll model is similar. Translation velocity is included into the differential equation describing the angular velocity dynamics in (2.9) and (2.11) [10, 12], which is different from the base model Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 10 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV (2.2). For a detailed consideration, it is proposed here to use a system of three differential equations for the angular velocity dynamics with the state variables: angle, angular velocity and translation velocity. The nonlinear model for translation velocity is obtained in [4]. Then, the pitch model is: v̇ = g tan θ − Cx|v|v, q̇ = −a1θ − a2q − a3v + a4τθ, θ̇ = q, (4.23) where v is horizontal translation velocity under assumptions that the vertical velocity is zero, i.e. the vertical projection of the thrust mg compensates to the quadcopter mass m, the pitch loop is considered separately, i.e. roll and yaw angles are zero and v corresponds to the pitch angle, and no wind, since the airspeed is implied in (4.23) but the ground speed is measured in the flight controller; g is acceleration of gravity; Cx is the aerodynamic parameter, which is not a constant, but can be considered as some value for the operating point; a1, a2, a3, a4 are other parameters of the model. The torque τθ is determined by the rotors dynamics (4.22) with some gain. The aerodynamic parameter Cx is the product of 1) the aerodynamic coefficient cx, which in turn depends on the velocity, 2) air density, i.e. the value ofCx may vary slightly depending on the weather conditions and decreases noticeably in the mountains, and 3) surface area, which depends on the tilt angle, which in the case of a quadcopter can be a very complex dependence, since the rotating propellers are not a rigid body. The experimental flight data contains the horizontal translation velocity values. An analysis of these data shows that for low velocity up to 10 m/s the nonlinear first formula in (4.23) is not confirmed, and the dynamics of the translation velocity depending on the angle can be described by the transfer function Pv(s) = g s+ a5 , (4.24) where g = 9.8065 and a5 is the parameter, which value is identified from the experimental flight data: a5 = 0.16. Then, the transfer function from the control uθ to the pitch angle formed from the system (4.23) taking into account (4.21) and (4.24) is Pθ(s) = ã4(s+ a5) (TΩs+ 1)((s+ a5)(s2 + a2s+ a1) + ã3) e−sTdelay , (4.25) where ã4 = a4 · K̃τ and ã3 = ga3 are new unknowns as products of the corresponding parameters of (4.23) and the gains of torque and translation velocity; TΩ and Tdelay are evaluated in (4.22), a5 is already identified too. Thus, the estimates of unknowns x = [a1, a2, ã3, ã4] of the model (4.25) should be found. The GlobalSearch optimization procedure from the MATLAB Global optimization Toolbox can be used to find these values from the experimentally obtained frequency parameters given in table 3.1. This procedure is based on the multistart approach [11] with nonlinear programming solver fmincon, where interior-point algorithm is chosen. This procedure usually provides a reliable result, while the most effective optimization procedure for the problem being solved is not discussed here. The objective function is min x nω∑ i=1 |P (x, jωi)− (αi + jβi)| |αi + jβi| . (4.26) Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 11 The estimates are found for roll and pitch: xϕ = [−1.28, 1.86, 73.07, 286.6], (4.27) xθ = [−14.61, 2.97, 83.06, 295.5]. (4.28) The identified transfer functions are Pϕ(s) = 56.67 · (s+ 0.16) (0.048s+ 1)(0.198s+ 1)(s2 − 3.04s+ 14.4) e−0.0042s, (4.29) Pθ(s) = 42.91 · (s+ 0.16) (0.048s+ 1)(0.145s+ 1)(s2 − 3.75s+ 11.7) e−0.0042s. (4.30) These functions match enough to the experimental data that is shown on Fig. 4.7, 4.8 and differ little from (3.18), (3.19). Fig. 4.7. Nyquist diagram of identified transfer function for roll The problem of the validity of the identified models under other flight conditions remains. The following flight test was carried out, where the required value of the pitch angle is formed as θ∗ = ηi sin(ωi(t− t0)) + θ0, where θ0 is an operating point. The values for the test are θ0 = 5°, ηi = 3°. The obtained frequency parameters estimates are in table 4.3. The estimates of the model (4.25) parameters other than (4.28) are found for the pitch loop: xθ = [−34.01, 3.41, 26.38, 241.9]. (4.31) The identified transfer function is Pθ(s) = 30.05 · (s+ 0.16) (0.048s+ 1)(0.124s+ 1)(s− 3.79)(s− 0.685) e−0.0042s. (4.32) The correspondence of the identified transfer function to the experimental data is shown on Fig. 4.9 Thus, the identified transfer functions (4.29), (4.30) represent the roll and pitch dynamics only for small values of the tilt angle and the translation velocity. The parameters a1 and a3 of (4.23) change significantly when the tilt angle is not near zero. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 12 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV Fig. 4.8. Nyquist diagram of identified transfer function for pitch Table 4.3. The frequency parameters estimates for the pitch loop (5°operating point) i Frequency, ωi, rad/s αi βi 1 20.11 -0.328 0.223 2 16.01 -0.562 0.304 3 12.03 -1.073 0.404 4 8.00 -1.971 0.390 5 5.34 -3.68 0.296 6 4.00 -5.003 0.545 7 3.50 -5.885 0.225 8 3.00 -5.735 0.672 9 2.00 -6.186 1.727 10 1.00 -4.554 3.292 Fig. 4.9. Nyquist diagram of identified transfer function for pitch (5°operating point) Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 13 4.3. Yaw model It should be noted that the structure of the identified transfer functions (4.29), (4.30) for roll and pitch angles does not contain an integration element, in contrast to the preliminary structure (2.6), while the integrator is present in (3.20) for the yaw loop. Thus, the transfer function for angular velocity r can be considered, and then, the transfer function for yaw angle is obtained by adding the integrator. This approach is further approved by the fact that in the case of yaw loop the measured data of the angular velocity during the sine wave test have less fluctuations than the angle data that makes it possible to obtain more reliable values of the frequency parameters. It is obvious that the rotors dynamics (4.22) should be used for the model from control uψ to torque τψ. The next idea follows from the presence of numerator in identified transfer function (3.20). This suggests the existence of some state variable, such that the dynamics of the angular velocity can be described as a system ν̇ = −a4ν +Kνr, ṙ = −a1ν − a2r + a3τψ, (4.33) where ν is some state variable, a1..4, Kν are model parameters. Then, the transfer function from the control uψ to the yaw angular velocity r has the form Pr(s) = ã3(s+ a4) (TΩs+ 1)((s+ a4)(s+ a2) + ã1) e−sTdelay , (4.34) where ã3 = a3 · K̃τ and ã1 = a1 ·Kν . The estimates of unknowns xψ = [ã1, a2, ã3, a4] are found via the GlobalSearch optimization procedure from the experimentally obtained frequency parameters for the yaw angular velocity: xψ = [−1875, 352.6, 992.0, 5.62]. (4.35) The identified transfer function is Pr(s) = 992 · (s+ 5.62) (0.048s+ 1)(s+ 358)(s+ 0.296) e−0.0042s. (4.36) The identified function (4.36) match enough to the experimental data that is shown on Fig. 4.10. The transfer function for yaw angle can be presented: Pψ(s) = 992 · (s+ 5.62) s(0.048s+ 1)(s+ 358)(s+ 0.296) e−0.0042s. (4.37) This function differs from (3.20), since the frequency parameters for yaw angle are less accurate than for yaw angular velocity one, and the result obtained with criterion (4.26) and structure (4.34) is more reasonable. 5. CONCLUSION Frequency domain identification allows to identify and analysis complex non-obvious dynamics. The transfer functions (4.29), (4.30), (4.37) for roll, pitch and yaw angles of the quadcopter attitude are obtained in this article. The advantage of separated identification for simple loops is shown, such as the rotors dynamics, to reduce the complicity of the non- obvious part of the system. Solution of the optimization problem can provide better results than the least squares due to the correct cost function and also allows to find the delay Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 14 V. ALEXANDROV, I. REZKOV, D. SHATOV, YU. MOROZOV Fig. 4.10. Nyquist diagram of identified transfer function for yaw angular velocity estimation. The unknown parameters for optimization procedure can be more sensible, than the coefficients of the transfer function. Identification of the pitch loop for operating points of 0° and 5° proves that the system is nonlinear, since the obtained transfer functions differ significantly. Further research is required to find limits and dependencies of the system parameters. ACKNOWLEDGEMENTS Research partially supported by RSF, project No. 23-29-00588, https://rscf.ru/en/project/23-29-00588/. REFERENCES 1. Alabsi, M.I. & Fields, T.D. (2019) Real-Time Closed-Loop System Identification of a Quadcopter. J. Aircraft, 56(1), 1–12. 2. Alexandrov, A.G. (1999) Finite-frequency identification and model validation of stable plant. IFAC Proceedings Volumes, 32(2), 4064–4069. 3. Alexandrov, A. & Palenov, M. (2011) Self-Tuning PID-I Controller. IFAC Proceedings Volumes, 44(1), 3635–3640. 4. Alexandrov, V., Rezkov, I. & Shatov, D. (2020) Finite-frequency Identification of the Quadcopter Translation Dynamics. Proc. of the 24th International Conference on System Theory, Control and Computing (ICSTCC), Sinaia, Romania, 471–476. 5. Alexandrov, V., Rezkov, I. & Shatov, D. (2021) Linearized model identification for quadcopter vertical translation dynamics. Proc. of the 25th International Conference on System Theory, Control and Computing (ICSTCC), Iasi, Romania, 278–283. 6. Alexandrov, V., Rezkov, I., Shatov, D. & Morozov, Y. (2023). Identification of the Quadcopter Rotational Dynamics for the Tilt Angle. In: Ronzhin, A., Sadigov, A., Meshcheryakov, R. (eds) Interactive Collaborative Robotics. ICR 2023. Lecture Notes in Computer Science (pp. 334–344), vol 14214. Springer, Cham. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) FREQUENCY DOMAIN IDENTIFICATION OF THE QUADCOPTER ATTITUDE DYNAMICS 15 7. Cho, S.H., Bhandari, S., Sanders, F.C., Tischler, M. & Cheung, K.K. (2019) System Identification and Controller Optimization of Coaxial Quadrotor UAV in Hover. Proc. of the AIAA Scitech 2019 Forum, San Diego, CA, 1–14. 8. Cunningham, M.A. & Hubbard Jr., J.E. (2020) Open-Loop Linear Model Identification of a Multi-rotor Vehicle with Active Feedback Control. J. Aircraft, 57(3), 1–18. 9. Murphy, M.D. (2016) A modular virtual laboratory for quadrotor control simulation. IFAC-PapersOnLine, 49(6), 93–98. 10. Niermeyer, P., Raffler, T. & Holzapfel, F. (2015) Open-Loop Quadrotor Flight Dynamics Identification in Frequency Domain via Closed-Loop Flight Testing. Proc. of the AIAA Guidance, Navigation, and Control Conference, Kissimmee, FL, 1–14. 11. Ugray, Z., Lasdon, L., Plummer, J., Glover, F., Kelly, J., et. al. (2007) Scatter search and local NLP solvers: a multistart framework for global optimization. INFORMS Journal on Computing, 19(3), 328–340. 12. Wei, W., Cohen, K., & Tischler, M.B. (2015). System identification and controller optimization of a quadrotor UAV. Proc. of the AHS International’s 71st Annual Forum and Technology Display, Virginia Beach, Virginia. 13. Wei, W., Tischler, M.B., Schwartz, N. & Cohen, K. (2014) Frequency-Domain System Identification and Simulation of a Quadrotor Controller. Proc. of the AIAA Modeling and Simulation Technologies Conference, National Harbor, MD, 1–9. 14. Zhang, X., Li, X., Wang, K. & Lu, Y. (2014) A survey of modelling and identification of quadrotor robot. Abstract and Applied Analysis, 2014, 320526. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) Introduction Problem statement FINITE FREQUENCY IDENTIFICATION Identification in details Rotors dynamics Roll and pitch model Yaw model Conclusion