- Email: [email protected]
Recursive Identification of Motion Model Parameters for ultralight UAV∗
Preprints, 1st IFAC Conference on Modelling, Identification and Preprints, 1st IFAC Conference on Modelling, Preprints, IFAC Conference Modelling, Identification Identification and and Control of 1st Nonlinear Systems on Preprints, 1st IFAC on Preprints, 1st IFAC Conference Conference on Modelling, Modelling, Identification Identification and and Control of Nonlinear Systems Control of Nonlinear Systems Available online at www.sciencedirect.com June 24-26, 2015. Saint Petersburg, Russia Control of Nonlinear Systems Control of Nonlinear Systems June 24-26, 2015. Saint Petersburg, Russia June 24-26, 24-26, 2015. Saint Saint Petersburg, Russia Russia June June 24-26, 2015. 2015. Saint Petersburg, Petersburg, Russia
ScienceDirect
IFAC-PapersOnLine 48-11 (2015) 233–237
Recursive Identification of Motion Model Recursive Identification of Motion Model Recursive Identification of Motion Model Recursive Identification of Motion Model Parameters for ultralight UAV. Parameters for ultralight UAV. Parameters Parameters for for ultralight ultralight UAV. UAV.
∗,∗∗∗ ∗∗,∗∗∗ Amelin ∗,∗∗∗ Stanislav Tomashevich ∗∗,∗∗∗ ∗,∗∗∗ ∗∗,∗∗∗ Amelin Tomashevich ∗,∗∗∗ Stanislav ∗∗,∗∗∗ ∗,∗∗∗ Amelin Stanislav Tomashevich ∗,∗∗∗ Amelin Stanislav Tomashevich Boris Andrievsky Amelin Stanislav∗,∗∗∗ Tomashevich ∗∗,∗∗∗ ∗,∗∗∗ Boris Andrievsky ∗,∗∗∗ Boris Andrievsky Andrievsky ∗,∗∗∗ Boris Boris Andrievsky ∗ ∗ Saint-Petersburg State University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg State State University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg State University (SPbSU), Russia (e-mail: [email protected]). Saint-Petersburg State University (SPbSU), Russia (e-mail: [email protected]). ∗∗ [email protected]). [email protected]). University, Russia (e-mail: ∗∗ ITMO [email protected]). ∗∗ ITMO University, Russia (e-mail: ∗∗ ITMO University, University, Russia (e-mail: (e-mail: ∗∗ ITMO [email protected]) ITMO University, Russia Russia (e-mail: [email protected]) ∗∗∗ [email protected]) [email protected]) Institute for Problems of Mechanical Engineering, the Russian ∗∗∗ [email protected]) ∗∗∗ for Problems of Mechanical Engineering, the Russian ∗∗∗ Institute Institute for Problems of Mechanical Engineering, the Russian ∗∗∗ Academy Institute for Problems of Mechanical Engineering, the of Sciences (e-mail: [email protected]) Institute for Problems(e-mail: of Mechanical Engineering, the Russian Russian Academy of Sciences [email protected]) Academy of Sciences (e-mail: [email protected]) Academy Academy of of Sciences Sciences (e-mail: (e-mail: [email protected]) [email protected]) Abstract: In the paper the possibility of using parameter estimations is examined. Identification Abstract: In In the the paper paper the the possibility possibility of of using using parameter parameter estimations estimations is is examined. examined. Identification Identification Abstract: Abstract: In the paper the possibility of using parameter estimations is examined. Identification algorithm is used to get quadrocopter model parameters. It is shown that this method is Abstract: In the paper the possibility of using parameter estimations is examined. Identification algorithm is used to get quadrocopter model parameters. It is shown that this method is is algorithm is used to getofquadrocopter quadrocopter model parameters. parameters. It configurations is shown shown that that were this method method algorithm is used to get model It is this is applicable for this type UAVs. Experiments with different conducted algorithm is used to get quadrocopter model parameters. It is shown that this method is applicable for for this this type type of of UAVs. UAVs. Experiments Experiments with with different different configurations configurations were were conducted conducted applicable applicable for type of UAVs. with configurations were conducted and described. The results show the efficiency of the applied identification applicable for this this of obtained UAVs. Experiments Experiments with different different configurations were procedure. conducted and described. described. The type results obtained show the the efficiency efficiency of the the applied applied identification procedure. and The results obtained show of identification procedure. and described. The results obtained show the efficiency of the applied identification procedure. Comparison with the results of the non-recursive least-squares identification scheme is presented. and described. The results obtained show the efficiency of theidentification applied identification procedure. Comparison with the results of the non-recursive least-squares scheme is presented. Comparison with with the the results results of of the the non-recursive non-recursive least-squares least-squares identification identification scheme scheme is is presented. presented. Comparison Comparison with the results of the non-recursive least-squares identification scheme is presented. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: UAV, quadrotor, identification, estimation, altitude hold. Keywords: UAV, UAV, quadrotor, quadrotor, identification, identification, estimation, estimation, altitude altitude hold. hold. Keywords: Keywords: UAV, quadrotor, identification, estimation, altitude hold. Keywords: UAV, quadrotor, identification, estimation, altitude hold. 1. INTRODUCTION The rest of the paper is organized as follows. The onboard 1. INTRODUCTION The rest paper is as follows. The onboard 1. The rest of of the thealgorithm paper is organized organized as follows. The onboard 1. INTRODUCTION INTRODUCTION The paper as The identification is described in Section 2. The 1. INTRODUCTION The rest rest of of the thealgorithm paper is is organized organized as follows. follows. The onboard onboard identification is described in Section 2. The identification algorithm is described in Section 2. The identification algorithm is described in Section 2. model which is used for identification, and procedure In recent times many ideas of using ultralight UAVs beidentification algorithm is identification, described in Section 2. The The model which is used for and procedure In recent times many ideas of using ultralight UAVs bemodel which is used for identification, and procedure In recent times many ideas of using ultralight UAVs bemodel which is used for identification, and procedure the quadrotor model identification based In recent times many ideas of of using ultralight UAVs bebe- of comes more and more global. These UAVs (quadrocopters) model which is used forparameters identification, and procedure In recent times many ideas using ultralight UAVs of the quadrotor model parameters identification based comes more and more global. These UAVs (quadrocopters) of the the quadrotor model parameters identification based comes more and These UAVs (quadrocopters) of quadrotor model parameters identification based on the flight test data (by LSE method see (Ljung, comes more advantages and more more global. global. These UAVs (quadrocopters) have many like mobility, low weight, speed of the quadrotor model parameters identification based comes more and more global. UAVs (quadrocopters) the flight test data (by LSE method see (Ljung, have many many advantages like aaaThese mobility, low weight, speed speed on on the flight test data (by LSE method see (Ljung, have advantages like mobility, low weight, on the flight test data (by LSE method see (Ljung, 1999)) is described in Section 3. Section 4 is devoted have many advantages like a mobility, low weight, speed and relative cheapness. Generally, the setting of control on the flight test data (by LSE method see (Ljung, have many advantages like a mobility, low weight, speed 1999)) is described in Section 3. Section 4 is devoted and relative relative cheapness. cheapness. Generally, Generally, the the setting setting of of control control 1999)) is described in Section to3. Section 44 is devoted and 1999)) is in Section devoted application of identification existing and relative cheapness. Generally, the setting in of manual control to system parameter is provided by operators 1999)) is described described in Section Section to3. 3. the Section 4 is isquadrotor devoted and relative cheapness. Generally, setting of control to application of identification the existing quadrotor system parameter is provided provided by the operators in manual to application of identification to the existing quadrotor system parameter is by operators in manual to application of identification to the existing quadrotor and preliminary estimation of the parameters. Results of system parameter is provided by operators in manual mode. Identification algorithm will allow to automate the to application of identification to the existing quadrotor system parameter isalgorithm providedwill by allow operators in manual and preliminary estimation of the parameters. Results mode. Identification to automate the and parameter preliminaryidentification estimation of ofare thedescribed parameters. Results of of mode. Identification algorithm and preliminary estimation the parameters. Results of in Section 5. mode. Identification algorithm will will allow allow to to automate automate the the the process of setting up. and preliminary estimation of the parameters. Results of mode. Identification algorithm will allow to automate the the parameter identification are described in Section 5. process of of setting setting up. up. the parameter identification are described in Section 5. process the parameter identification are described in Section 5. Concluding remarks and the future work intentions are process of setting up. the parameter identification are described in Section 5. process of setting up. Concluding remarks and the future work intentions are The knowledge of models, their parameters and dynamics Concluding remarks and the future work intentions are Concluding remarks and the future work intentions are given in Section 6. The knowledge of models, their parameters and dynamics Concluding remarks and the future work intentions are The knowledge of models, their parameters and dynamics given in Section 6. The knowledge of models, their parameters and dynamics should help to designers make control laws, given in in Section 6. 6. The knowledge of models,to andcontrollers dynamics given should help to to designers designers totheir makeparameters control laws, laws, controllers given in Section Section 6. should help to make should help toaccurate designersand to efficient. make control control laws, controllers controllers and etc more Moreover, when it is should help to designers to make control laws, controllers and etc etc more more accurate and and efficient. Moreover, Moreover, when when it is is and and etconboard more accurate accurate and efficient. efficient. Moreover, when it it to is doing and in real time, it gives aa possibility 2. IDENTIFICATION ALGORITHM and more accurate and efficient. Moreover, when it is doingetconboard onboard and in in real real time, it it gives gives possibility to 2. IDENTIFICATION ALGORITHM doing and time, a possibility to 2. IDENTIFICATION ALGORITHM doing onboard and inin inthe realparameters time, it it gives gives possibility to respond to changes online. For example, 2. IDENTIFICATION ALGORITHM doing onboard and real time, aa possibility to 2. IDENTIFICATION ALGORITHM respond to changes in the parameters online. For example, respond to changes in the parameters online. For example, respond to changes in the parameters online. For example, if UAVs mass is changing it is needed to build a new control respond to changes in the parameters online. For example, There are many algorithms to get model parameters, see, if UAVs mass is changing it is needed to build a new control if UAVs mass it to There are to parameters, see, if UAVs mass is is changing changing it is is needed needed to build build aaa new new control control e.g. law or reconfigure controllers parameters. There are many many algorithms algorithms to get get model model parameters, see, if UAVs mass is changing it is needed to build new control There are many algorithms to get model parameters, see, (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure reconfigure controllers parameters. There are many algorithms to get model parameters, see, law or controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the interest to use algorithm, which measures few information There are many applications when quadrocopters are e.g. (Kalouptsidis and Theodoridis, 1993). There is the interest to use algorithm, which measures few information There are many applications when quadrocopters are interest to use algorithm, which measures few information interest to use algorithm, which measures few information There are many applications when quadrocopters are as possible, because we can not know, for example, full There are many applications when quadrocopters are needed: ultralight UAVs can perform many different tasks interest to use algorithm, which measures few information possible, because we can not know, for example, full There are many UAVs applications whenmany quadrocopters are as needed: ultralight can perform different tasks as possible, possible, because we range can not not know, for example, example, full as because we can know, for full needed: ultralight can many different tasks state of plant in aa wide of tasks. Therefore we use needed: ultralight UAVs UAVs can perform perform many different tasks from photography the terrain to carry out flight in the as possible, because we range can not know, for example, full state of plant in wide of tasks. Therefore we use needed: ultralight UAVs can perform many from photography photography the terrain terrain to carry carry outdifferent flight in intasks the identification state of plant in a wide range of tasks. Therefore we use state of plant in a wide range of tasks. Therefore we use from the to out flight the algorithm described in (Gawthrop, 1987) befrom photography the terrain to carry out flight in the mountains. The main condition of all this tasks is stable state of plant in a wide range of tasks. Therefore we identification algorithm described in (Gawthrop, 1987) befrom photography thecondition terrain toofcarry outtasks flightis in the identification algorithm described in (Gawthrop, 1987) use mountains. The main all this stable beidentification algorithm described in (Gawthrop, 1987) bemountains. The main condition of all this tasks is stable cause it is based on knowledge about input-output signals mountains. The main condition of all this tasks is stable and precise altitude holding. When aircraft is flying at an identification algorithm described in (Gawthrop, 1987) because it is based on knowledge about input-output signals mountains. The main condition of all this tasks is stable and precise altitude holding. When aircraft is flying at an cause it is based on knowledge about input-output signals cause it is based on knowledge about input-output signals and precise altitude holding. When aircraft is flying at an and about orders of transfer function only. Let us abbreand precise altitude holding. When aircraft is flying at an altitude of several thousand meters above terrain without cause it is based on knowledge about input-output signals about orders of transfer function only. Let us abbreand precise altitudethousand holding. meters When aircraft is flying at an and altitude of several above terrain without and about about orders of of transfer function only. only. Let us abbreabbreand orders transfer function us altitude of thousand meters terrain viate this method as CIA (Continuous-time Identification altitude of several several thousand meters above abovebecause terrain without without mountains there is no such requirement error of and about orders of transfer function only. Let Let us abbreviate this method as CIA (Continuous-time Identification altitude of several meters above terrain without mountains there is isthousand no such such requirement requirement because error of of Algorithm). viate this method as CIA (Continuous-time Identification viate this method as CIA (Continuous-time Identification mountains there no because error It is easy to see that it had a good approach mountains thereis is is no no such such requirement because error of Algorithm). several meters for this But when viate this method as CIA (Continuous-time Identification It is easy to see that it had aa good approach mountains there because error of several meters meters is not not critical criticalrequirement for this this altitude. altitude. But when Algorithm). It is easy to see that it had approach Algorithm). It is easy to see that it had aa good good approach several is for But when to identification for “LAAS Helicopter Benchmark” with several meters is not not critical critical for over this altitude. altitude. But when aircrafts flight trajectory lies low the ground error in Algorithm). It is easy to see that it had good approach to identification for “LAAS Helicopter Benchmark” with several meters is not critical for this altitude. But when aircrafts flight trajectory lies low over the ground error in to identification for “LAAS Helicopter Benchmark” with to identification for “LAAS Helicopter Benchmark” with aircrafts flight trajectory lies low over the ground error in a good quality, see, e.g. (Fradkov et al., 2010). This alaircrafts flight trajectory lies low over the ground error in several meters may cause a crash. to identification for “LAAS Helicopter Benchmark” with quality, see, e.g. (Fradkov et al., 2010). This alaircrafts flight trajectory lies low over the ground error in aa good several meters may cause a crash. good quality, see, e.g. (Fradkov et al., 2010). This ala good quality, see, e.g. (Fradkov et al., 2010). This alseveral meters may cause a crash. gorithm is described in Section 2. It would be helpful to several meters may cause a crash. a good quality, see, e.g. (Fradkov et al., 2010). This algorithm is described in Section 2. It would be helpful to several meters may cause a crash. gorithm is described in Section 2. It would be helpful to Identification problem can be solved by two ways: ongorithm is described in Section 2. It would be helpful to compare results with another way to identify plant model Identification problem can be solved by two ways: ongorithm is described in Section 2. It would be helpful to compare results with another way to identify plant model Identification problem can be solved by two ways: onIdentification problem can be solved by two ways: oncompare results with another way to identify plant model board in real time and after receiving all data results. compare results with another way to identify plant model Identification problem can be receiving solved byalltwo ways: on- parameters with using non-recursive least squares identiboard in real time and after data results. compare results with another way to identify plant model parameters with using non-recursive least squares identiboard in real time and after receiving all data results. board in real real by timefirst andway after receiving all much data results. results. parameters with using using non-recursive least squares identiidentiIdentification is assumed in papers parameters with non-recursive least squares board in time and after all data algorithm which is described in subsection 3.2 Identification by first way is receiving assumed in in much papers fication parameters with using non-recursive least squares identialgorithm which is described in subsection 3.2 Identification by way is papers Identification bytofirst first waydata is assumed assumed in much much papers fication fication algorithm which is described in subsection 3.2 but it is easier collect after flying and simulate fication algorithm which is described in subsection 3.2 Identification by first way is assumed in much papers but it is easier to collect data after flying and simulate fication algorithm which is described in subsection 3.2 but it is easier to collect data after flying and simulate Consider the following single-input–single-output (SISO) but it is easier to collect data after flying and simulate this identification process separately (how we make in this Consider the following single-input–single-output (SISO) but it is easier to collect data after flying and simulate this identification identification process process separately separately (how (how we we make make in in this this plant Consider the following single-input–single-output (SISO) Consider the following single-input–single-output (SISO) this model: this identification process separately (howStanculeanu we make make in in and this plant paper). Papers (Abas et al., 2011, 2012; Consider the following single-input–single-output (SISO) this identification process separately (how we this paper). Papers (Abas (Abas et al., al., 2011, 2012; 2012; Stanculeanu and plant model: model: plant model: paper). Papers et 2011, Stanculeanu and paper). Papers (Abas et al., 2011, 2012; Stanculeanu and Borangiu, 2011) are devoted to second way. plant model: paper). Papers et al., 2011, 2012;way. Stanculeanu and Borangiu, 2011)(Abas are devoted to second Borangiu, (n−1) Borangiu, 2011) 2011) are are devoted devoted to to second second way. way. yy (n) Borangiu, are devoted to second way. (n) (t) + a1 y (n−1) (t) + . . . + an y(t) = The work 2011) (n) (n−1) (t) + . . . + an y(t) = + a 1y (n) (t) was performed in IPME RAS, supported by RSF (grant y (t) + a y (t) + . .. .. + a 1 y (n−1) (n−1) yy (n)(m) (t) + a (t) + ann y(t) y(t) = = The work was performed in IPME RAS, supported by RSF (grant 1 (m−1) (t) + a y (t) (t) + ..+ . ..+ + = The work was performed in IPME RAS, supported by RSF (grant 1 n y(t) u (t) + b .. .. a+ bm u(t) (1) b 14-29-00142). The work was performed in IPME RAS, supported by RSF (grant (m) (m−1) The work was performed in IPME RAS, supported by RSF (grant 0 u(m) (t) + b1 u (m−1) u (t) + . + (1) b 14-29-00142). 0 u(m) (t) + b1 u(m−1) (t) + . . . + b m u(t) b u(t) (1) b 14-29-00142). 0 u(m) (t) + b1 1 u(m−1) (t) + . . . + bm m u(t) (1) b 14-29-00142). 0 (t) + . . . + b u(t) (1) b u (t) + b u 14-29-00142).
Konstantin Konstantin Konstantin Konstantin Konstantin
0
1
Copyright © 2015, IFAC IFAC 2015 (International Federation of Automatic Control) 237Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © Copyright © IFAC 2015 237 Copyright ©under IFAC responsibility 2015 237Control. Copyright 237 Peer review© of International Federation of Automatic Copyright © IFAC IFAC 2015 2015 237 10.1016/j.ifacol.2015.09.189
m
MICNON 2015 234 Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237 June 24-26, 2015. Saint Petersburg, Russia
where ai , i = 1, 2, . . . , n and bj , j = 1, 2, . . . , m are unknown parameters which are needed to be identified. y(t), u(t) is output and input respectively. Let us introduce two state filters ζ˙ = Ad ζ + bd y(t), σ˙ = Ad σ + bd u(t)
(2) where ζ(t), σ(t) ∈ R ; pair (Ad , bd ) has a regular canonical form and det(sIn −Dd ) = D(s) is an arbitrary Hurwitz polynomial D(s) = sn + d1 sn−1 + . . . + dn . N
Identification algorithm has the following form ˜ ˙ ρ(t)ζ(t) Ω(t) = −Γ(t)˜ ρ(t)˜ ρ(t)T Ω(t) + Γ(t)˜ ˙ Γ(t) = −Γ(t)˜ ρ(t)˜ ρ(t)T Γ(t) + αΓ(t)
(3) (4)
n+m+1
where Ω(t) ∈ R – vector of estimates for ai , i = 1, 2, . . . , n and bj , j = 1, 2, . . . , m; Γ(t) is gain matrix; ρ˜ = [ζn , . . . , ζ1 , σm+1 , . . . , σ1 ] is vector of filters (2) n outputs; ζ˜ = y(t) − dn−i+1 ζi and α > 0 is algorithm i=1
parameter. It is needed to take initial value Γ(0) = Γ(0)T and this parameter is responsible for identification speed. 3. QUADROTOR MODEL 3.1 Design modeling dynamics
For quadrocopters the problem of altitude hold is simple but requires precision and reliability. The quadrotor helicopter dynamics model may be derived using EulerLagrange formalism under the assumptions that the frame structure is rigid, the structure is symmetrical, the centre of gravity (CoG) and the body fixed frame origin are assumed to coincide and the propellers are rigid. There is significant amount of literature dedicated to quadrotors, where the model of platform dynamics is derived, see papers by Fantoni and Lozano (2002); Altug et al. (2002); Castillo et al. (2005); Das et al. (2008); Janusz et al. (2013) for mentioning a few. This model has the following form: Ix p˙ = (Iy − Iz )qr + Tϕ (t), Iy q˙ = (Iz − Ix )pr + Tθ (t), Iz r˙ = (Ix − Iy )pq + Tψ (t), (5) m¨ x = (sin ψ · sin ϕ + cos ψ · sin θ · cos ϕ)F (t), m¨ y = (sin ψ · sin θ cos ϕ − cos ψ · sin ϕ)F (t), m¨ z = cos θ · cos ϕ · F (t) − mg,
where ϕ, θ, ψ are the Euler angles (roll, pitch, yaw, respectively); p, q, r – angular velocities; (x, y, z) gives the quadrotor CoG position in the normal Earth reference frame, Ix , Iy , Iz stand for inertia moments; Tϕ (t), Tθ (t), Tψ (t) are external torques, acting to the quadrotor; m denotes the quadrotor mass; g is the gravity acceleration. F (t) denotes the following total thrust, produced by the rotors: F (t) = f1 (t) + f2 (t) + f3 (t) + f4 (t), (6) where fi (t), i = 1, . . . 4, are forces, applied from each one rotor. It is worth mentioning that fi (t) depend on the blades rotation speed, aerodynamic drag and external disturbances. 238
In the present paper we focus our attention to quadrotor movement in height, assuming that the quadrotor is horizontally positioned. Therefore angles ϕ, θ in (5) should be taken zero, yaw angle ψ may have an arbitrary value. This leads to the following equation for isolated vertical motion dynamics: m¨ z = F (t) − mg, (7) where z denotes the quadrotor altitude. Interaction between rotating blades and surrounding air has a complex nature, cf. (Leishman, 2006; Janusz et al., 2013). Based on (Janusz et al., 2013) we assume that this interaction may be modeled as a viscous friction of the vertical motion. Then the vertical motion dynamics may be described as follows: 4 ˙ = Ku Ui (t) − mg, (8) m¨ z (t) + Kz˙ z(t) i=1
where Ui denotes the rotation speed of ith propeller (i = 1, . . . , 4), Kz˙ is the viscous friction parameter, Ku may be called the propeller effectiveness gain.
In what follows, due to the assumption that the quadrotor is positioned horizontally, it is taken that U1 (t) = . . . U4 (t). Provisionally neglecting the brushless motors dynamics (which are fast in compare with the quadrotor body dynamics), we obtain the following model of the vertical motion dynamics z¨(t) + K z(t) ˙ = Ku Uz (t) − g, (9) 4 where Uz = Ui is the altitude control signal (the sum i=1
of PWM signals Ui , applied to the motors), K, Ku are normalized model parameters. K, Ku depend on various factors and they are difficult to calculate. Therefore we use an identification procedure to find them based on the experimental results
In this way altitude hold is needed in working of servos: their generated traction must compensate gravity. To make this it is available to use PID-controller as example. 3.2 Estimate by LSE Coordinates z, acceleration z¨ and control signal Uz were obtained according to the results of the flight from the autopilot dataflash. With the help of differentiation by Eulers method the velocities and the accelerations vector z¨ can be obtain on the basis of the vector z or z˙ can be obtained by integration z¨. Moreover, the production of the vector z˙ by differentiation practically coincides with the data obtained from the quadrotors accelerators. The equation describing the dynamics (9) can be rewritten in matrix form χ = Y, (10) where Y = row(¨ z (1), z¨(2), . . . , z¨(n)), χ has the form Uz (1) 1 z(1) ˙ ˙ Uz (2) 1 z(2) , (11) χ= .. .. ... . . Uz (n) 1 z(n) ˙ where n – number of measurements. – vector of unknown parameters. From the expression (10) vector of unknown parameters is found as = χ+ Y , where χ+ – pseudoinverse matrix to matrix χ.
MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237
235
Let us identify parameters of quadrotor with model (9). In this subsection quadrotor model was found only for test No. 1 to show a rough estimate of the identification parameters (parameters of quadrotor for this case can be seen in Tab. 1). Result of identification by this method for other experiments is shown in Tab. 2. Following parameters have been found: Ku = 0.0208, g = 9.5841, K = 0.4631, then the dynamic equation (9) takes the form (12) z¨ + 0.4631z˙ = 0.0208Uz − 9.5841. From (12) it can be seen that the obtained acceleration of gravity close to how it should be gr ≈ 9.81, so we can assume it g = 9.81. 4. THE PRACTICAL APPLICATION FOR QUADROCOPTER For experiments four configuration are used, all variants is shown in Table 1 Table 1. Configurations parameters of quadrocopter. Test No. 1 2 3 4
Battery, mAh 4000 2600 2600 4000
Prop, inch 10 × 4.5 10 × 4.5 8×6 8×6
Mass, g 1220 1120 1118 1218
Fig. 3. Velocity z˙ of quadrotor. where Kz Uz = F1 + F2 + F3 + F4 and Kz is conversion factor. This equation is given for moments in time, when z¨ = 0, z˙ = 0 and ϕa = ψa = 0 because we consider the movement in a vertical plane. The equation (13) allows us to estimate conversion factor Kz = mg Uz ≈ 0.018. This estimation coincides with that obtained by the least squares method in the model (12). In the light of the above, equation (9) takes the form (14) m¨ z + K z˙ = 0.018Uz − m9.81. Knowledge of the equation (14) simplify identification, because this is necessary to evaluate one parameter less. But it will be shown below that this parameter Kz can be estimate too, but with some problems. 5. RESULTS OF IDENTIFICATION
Altitude is taken as the output value. These measurements are red from GPS and barometer, then mean value is used. The input signal is sum of PWM signals which are sending on each of motors. This signal is nondimensional quantity: it takes values from 1118 to 1919 for each motor and converts in autopilot with the following law: at the minimum (1118) rotation speed is zero and at the maximum (1919) rotation speed is 920 rmp/V × 12V , linear dependence. Input signal for one case is presented on Figure 1, acceleration z¨ and velocity z˙ are presented on Fig. 2 and on Fig 3 respectively.
Using identification algorithm which is described in Section 2, for four variations of equipment gives us result which is shown in the Table 2. It should be noted too that gravitational constant is known and this knowledge can simplify the identification because we need to estimate only one parameter: coefficient K. Since we use the method which is described in Section 2 it is needed to choose parameters for identification algorithm. Let us use next settings: α = 0.001, Γ(0) = 1000, ρ˜(0) = col{0, 0, 0}. Since there was four different experiments with various configurations and it means that it is necessary to allocate equal time intervals and identify parameters for these intervals to compare setting process later. These time histories is presented on Figs. 4-6.
Fig. 1. Input signal Uz . Fig. 4. Time histories of identification parameter g. Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively.
Fig. 2. Acceleration z¨ of quadrotor. As we can see, when input signal Uz ≈ 5800, then z¨ = 0, z˙ = 0. It means that equation (9) can be represented as follows Ku Uz = mg (13) 239
It can be seen that the parameters identifications for all four cases are the same and almost equal. There are some problems with parameter K: identification has a small variation and does not converge as the two other parameters, but it can be estimated separately as it is stated above. Finally, it is needed to show how the identified model reproduces input signal for one of the cases (test No 2).
MICNON 2015 236 Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237 June 24-26, 2015. Saint Petersburg, Russia
by LSE. Quadrocopters model dynamics for motion in a vertical plane are obtained. The results obtained are compared with the theoretical statements and efficiency of the resulting model is demonstrated. The difficulties encountered in the identification are noted.
Fig. 5. Time histories of identification parameter K. Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively.
The future work intentions are decide difficulties with equipment and make possible to have more accurately onboard identification. Also it is planned to identify other parameters quadrocopter to receive full model. REFERENCES
Fig. 6. Time histories of identification parameter Ku . Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively. Let us use model (15) to show efficiency of the resulting model (15) z¨ + 0.804z˙ = 0.016Uz − 9.744. Result of representing can be seen on Fig. 7. There is a good match with modeling output signal with the resulting model and real output.
Fig. 7. Output signal z, estimate z end error – blue solid line, black solid line and dashed line respectively. In the Figure 7 solid blue line is quadrocopter real altitude which was taken on certain time interval for 40 seconds, solid blue line is estimation of the altitude, obtained using the identified model and dashed line is error e(t) = z(t) − z(t). There are some problems with identification: input signals were very noisy at some time intervals and sensors accumulated errors over time, making it difficult to identify UAVs dynamic at large intervals of time. This is the reason why problems with identification of parameter K were found. 6. CONCLUSION Parameter identification is studied experimentally by two methods: by continuous-time identification algorithm and Table 2. The estimates obtained by two methods Test No. 1 2 3 4
g, m/s2 CIA LSE 9.638 8.516 9.071 7.945 11.213 9.972 11.714 8.313
K, 1/s CIA LSE 0.424 0.758 0.301 0.704 0.406 0.456 0.863 0.848
Ku CIA LSE 0.015 0.014 0.012 0.013 0.016 0.014 0.017 0.013
240
Abas, N., Legowo, A., and Akmeliawati, R. (2011). Parameter identification of an autonomous quadrotor. In 2011 4th International Conference On Mechatronics (ICOM), 1–8. Abas, N., Legowo, A., Ibrahim, Z., Rahim, N., and Kassim, A.M. (2012). Modeling and system identification using extended kalman filter for a quadrotor system. In Proceedings of 2012 International Conference on Electrical Engineering and Applications (ICEEA 2012). Altug, E., Ostrowski, J., and Mahony, R. (2002). Control of a quadrotor helicopter using visual feedback. In Proc. IEEE Int. Conf. on Robotics and Automation, 72–76. IEEE Press, Washington, DC, USA. Amelin, K., Amelina, N., Granichin, O., Granichina, O., and Andrievsky, B. (2013). Randomized algorithm for uavs group flight optimization. In Proc. of 11th IFAC International Workshop on Adaptation and Learning in Control and Signal Processing, 205–208. Castillo, G.P., Lozano, R., and Dzul, A.E. (2005). Modelling and Control of Mini-Flying Machines. Series: Advances in Industrial Control. Springer-Verlag Ltd. Das, A., Subbarao, K., and Lewis, F. (2008). Dynamic inversion of quadrotor with zero-dynamics stabilization. In Proc. 17th IEEE Int. Conf. on Control Applications, Part of 2008 IEEE Multi-conference on Systems and Control (MSC 2008). IEEE Press, San Antonio, USA. Fantoni, I. and Lozano, R. (2002). Non-linear Control for Underactuated Mechanical Systems. Series: Communications and Control Engineering. Springer-Verlag Ltd., London. Fradkov, A.L., Andrievsky, B., and Peaucelle, D. (2010). Estimation and control under information constraints for laas helicopter benchmark. Control Systems Technology, IEEE Transactions on, 18(5), 1180–1187. Gawthrop, P. (1987). Continuous-Time Self-Tuning Control. Letchworth, U.K.: Research Studies Press. Janusz, W.C., Czyba, R., and Szafra´ nski, G. (2013). Model identification and data fusion for the purpose of the altitude control of the VTOL aerial robot. In Proc. 2nd IFAC Workshop on Research, Education and Development of Unmanned Aerial Systems (RED-UAS 2013), 263–269. IFAC, Compiegne, France. Kalouptsidis, N. and Theodoridis, S. (eds.) (1993). Adaptive system identification and signal processing algorithms. Prentice-Hall, Inc., Upper Saddle River, NJ, USA. Leishman, J. (2006). Principles of Helicopter Aerodynamics. 2nd ed. Cambridge University Press, New York. Ljung, L. (1999). System Identification: Theory for the User. Prentice-Hall information and system sciences series. PTR Prentice Hall.
MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237
Stanculeanu, I. and Borangiu, T. (2011). Quadrotor blackbox system identification. In World Academy of Science, Engineering and Technology, volume 5, 276–279. Appendix A. TECHNICAL CHARACTERISTICS We use the ultralight UAV which is described in (Amelin et al., 2013). Autopilot Ardupilot Mega 2.6 with inertial system (gyroscope, accelerometer, magnetometer, barometer) and receiver signal from the satellite navigation system GPS. Autopilot Ardupilot Mega 2.6 is built on chips Atmels ATMEGA 2560, has 8 inputs and outputs for the PWM signal, three UART port, a built-in FTDI, 8 analog inputs, 4 Mb Flash memory for recording log data telemetry. It is equipped with Invensenses 6 DoF Accelerometer/Gyro MPU-6000, a barometer MS5611-01BA03 and 4 sonars MB1260 XL-MaxSonar-EZL0. Corps UAV – quadrocopter frame DJI 450. Driving body – radial-beam with the size of 450 mm in diameter. At the end of each of the four beams is a brushless motor with an external rotor and the turns controller (control mechanisms). The main UAV characteristics are as follows: • • • • •
The length of the beams, mm: 225; Flight weight, gram: varied in experiments; Flight speed, km/h: from 0 to 10; Engines, watts: 250 each of the four engines; Battery, mAh: 2600 or 4000 (depending on the experiment); • Autopilot and additional microcomputer: installed.
Fig. A.1. Photo of quadrocopter, which was used for the tests.
241
237
ScienceDirect
IFAC-PapersOnLine 48-11 (2015) 233–237
Recursive Identification of Motion Model Recursive Identification of Motion Model Recursive Identification of Motion Model Recursive Identification of Motion Model Parameters for ultralight UAV. Parameters for ultralight UAV. Parameters Parameters for for ultralight ultralight UAV. UAV.
∗,∗∗∗ ∗∗,∗∗∗ Amelin ∗,∗∗∗ Stanislav Tomashevich ∗∗,∗∗∗ ∗,∗∗∗ ∗∗,∗∗∗ Amelin Tomashevich ∗,∗∗∗ Stanislav ∗∗,∗∗∗ ∗,∗∗∗ Amelin Stanislav Tomashevich ∗,∗∗∗ Amelin Stanislav Tomashevich Boris Andrievsky Amelin Stanislav∗,∗∗∗ Tomashevich ∗∗,∗∗∗ ∗,∗∗∗ Boris Andrievsky ∗,∗∗∗ Boris Andrievsky Andrievsky ∗,∗∗∗ Boris Boris Andrievsky ∗ ∗ Saint-Petersburg State University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg State State University (SPbSU), Russia (e-mail: ∗ Saint-Petersburg State University (SPbSU), Russia (e-mail: [email protected]). Saint-Petersburg State University (SPbSU), Russia (e-mail: [email protected]). ∗∗ [email protected]). [email protected]). University, Russia (e-mail: ∗∗ ITMO [email protected]). ∗∗ ITMO University, Russia (e-mail: ∗∗ ITMO University, University, Russia (e-mail: (e-mail: ∗∗ ITMO [email protected]) ITMO University, Russia Russia (e-mail: [email protected]) ∗∗∗ [email protected]) [email protected]) Institute for Problems of Mechanical Engineering, the Russian ∗∗∗ [email protected]) ∗∗∗ for Problems of Mechanical Engineering, the Russian ∗∗∗ Institute Institute for Problems of Mechanical Engineering, the Russian ∗∗∗ Academy Institute for Problems of Mechanical Engineering, the of Sciences (e-mail: [email protected]) Institute for Problems(e-mail: of Mechanical Engineering, the Russian Russian Academy of Sciences [email protected]) Academy of Sciences (e-mail: [email protected]) Academy Academy of of Sciences Sciences (e-mail: (e-mail: [email protected]) [email protected]) Abstract: In the paper the possibility of using parameter estimations is examined. Identification Abstract: In In the the paper paper the the possibility possibility of of using using parameter parameter estimations estimations is is examined. examined. Identification Identification Abstract: Abstract: In the paper the possibility of using parameter estimations is examined. Identification algorithm is used to get quadrocopter model parameters. It is shown that this method is Abstract: In the paper the possibility of using parameter estimations is examined. Identification algorithm is used to get quadrocopter model parameters. It is shown that this method is is algorithm is used to getofquadrocopter quadrocopter model parameters. parameters. It configurations is shown shown that that were this method method algorithm is used to get model It is this is applicable for this type UAVs. Experiments with different conducted algorithm is used to get quadrocopter model parameters. It is shown that this method is applicable for for this this type type of of UAVs. UAVs. Experiments Experiments with with different different configurations configurations were were conducted conducted applicable applicable for type of UAVs. with configurations were conducted and described. The results show the efficiency of the applied identification applicable for this this of obtained UAVs. Experiments Experiments with different different configurations were procedure. conducted and described. described. The type results obtained show the the efficiency efficiency of the the applied applied identification procedure. and The results obtained show of identification procedure. and described. The results obtained show the efficiency of the applied identification procedure. Comparison with the results of the non-recursive least-squares identification scheme is presented. and described. The results obtained show the efficiency of theidentification applied identification procedure. Comparison with the results of the non-recursive least-squares scheme is presented. Comparison with with the the results results of of the the non-recursive non-recursive least-squares least-squares identification identification scheme scheme is is presented. presented. Comparison Comparison with the results of the non-recursive least-squares identification scheme is presented. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: UAV, quadrotor, identification, estimation, altitude hold. Keywords: UAV, UAV, quadrotor, quadrotor, identification, identification, estimation, estimation, altitude altitude hold. hold. Keywords: Keywords: UAV, quadrotor, identification, estimation, altitude hold. Keywords: UAV, quadrotor, identification, estimation, altitude hold. 1. INTRODUCTION The rest of the paper is organized as follows. The onboard 1. INTRODUCTION The rest paper is as follows. The onboard 1. The rest of of the thealgorithm paper is organized organized as follows. The onboard 1. INTRODUCTION INTRODUCTION The paper as The identification is described in Section 2. The 1. INTRODUCTION The rest rest of of the thealgorithm paper is is organized organized as follows. follows. The onboard onboard identification is described in Section 2. The identification algorithm is described in Section 2. The identification algorithm is described in Section 2. model which is used for identification, and procedure In recent times many ideas of using ultralight UAVs beidentification algorithm is identification, described in Section 2. The The model which is used for and procedure In recent times many ideas of using ultralight UAVs bemodel which is used for identification, and procedure In recent times many ideas of using ultralight UAVs bemodel which is used for identification, and procedure the quadrotor model identification based In recent times many ideas of of using ultralight UAVs bebe- of comes more and more global. These UAVs (quadrocopters) model which is used forparameters identification, and procedure In recent times many ideas using ultralight UAVs of the quadrotor model parameters identification based comes more and more global. These UAVs (quadrocopters) of the the quadrotor model parameters identification based comes more and These UAVs (quadrocopters) of quadrotor model parameters identification based on the flight test data (by LSE method see (Ljung, comes more advantages and more more global. global. These UAVs (quadrocopters) have many like mobility, low weight, speed of the quadrotor model parameters identification based comes more and more global. UAVs (quadrocopters) the flight test data (by LSE method see (Ljung, have many many advantages like aaaThese mobility, low weight, speed speed on on the flight test data (by LSE method see (Ljung, have advantages like mobility, low weight, on the flight test data (by LSE method see (Ljung, 1999)) is described in Section 3. Section 4 is devoted have many advantages like a mobility, low weight, speed and relative cheapness. Generally, the setting of control on the flight test data (by LSE method see (Ljung, have many advantages like a mobility, low weight, speed 1999)) is described in Section 3. Section 4 is devoted and relative relative cheapness. cheapness. Generally, Generally, the the setting setting of of control control 1999)) is described in Section to3. Section 44 is devoted and 1999)) is in Section devoted application of identification existing and relative cheapness. Generally, the setting in of manual control to system parameter is provided by operators 1999)) is described described in Section Section to3. 3. the Section 4 is isquadrotor devoted and relative cheapness. Generally, setting of control to application of identification the existing quadrotor system parameter is provided provided by the operators in manual to application of identification to the existing quadrotor system parameter is by operators in manual to application of identification to the existing quadrotor and preliminary estimation of the parameters. Results of system parameter is provided by operators in manual mode. Identification algorithm will allow to automate the to application of identification to the existing quadrotor system parameter isalgorithm providedwill by allow operators in manual and preliminary estimation of the parameters. Results mode. Identification to automate the and parameter preliminaryidentification estimation of ofare thedescribed parameters. Results of of mode. Identification algorithm and preliminary estimation the parameters. Results of in Section 5. mode. Identification algorithm will will allow allow to to automate automate the the the process of setting up. and preliminary estimation of the parameters. Results of mode. Identification algorithm will allow to automate the the parameter identification are described in Section 5. process of of setting setting up. up. the parameter identification are described in Section 5. process the parameter identification are described in Section 5. Concluding remarks and the future work intentions are process of setting up. the parameter identification are described in Section 5. process of setting up. Concluding remarks and the future work intentions are The knowledge of models, their parameters and dynamics Concluding remarks and the future work intentions are Concluding remarks and the future work intentions are given in Section 6. The knowledge of models, their parameters and dynamics Concluding remarks and the future work intentions are The knowledge of models, their parameters and dynamics given in Section 6. The knowledge of models, their parameters and dynamics should help to designers make control laws, given in in Section 6. 6. The knowledge of models,to andcontrollers dynamics given should help to to designers designers totheir makeparameters control laws, laws, controllers given in Section Section 6. should help to make should help toaccurate designersand to efficient. make control control laws, controllers controllers and etc more Moreover, when it is should help to designers to make control laws, controllers and etc etc more more accurate and and efficient. Moreover, Moreover, when when it is is and and etconboard more accurate accurate and efficient. efficient. Moreover, when it it to is doing and in real time, it gives aa possibility 2. IDENTIFICATION ALGORITHM and more accurate and efficient. Moreover, when it is doingetconboard onboard and in in real real time, it it gives gives possibility to 2. IDENTIFICATION ALGORITHM doing and time, a possibility to 2. IDENTIFICATION ALGORITHM doing onboard and inin inthe realparameters time, it it gives gives possibility to respond to changes online. For example, 2. IDENTIFICATION ALGORITHM doing onboard and real time, aa possibility to 2. IDENTIFICATION ALGORITHM respond to changes in the parameters online. For example, respond to changes in the parameters online. For example, respond to changes in the parameters online. For example, if UAVs mass is changing it is needed to build a new control respond to changes in the parameters online. For example, There are many algorithms to get model parameters, see, if UAVs mass is changing it is needed to build a new control if UAVs mass it to There are to parameters, see, if UAVs mass is is changing changing it is is needed needed to build build aaa new new control control e.g. law or reconfigure controllers parameters. There are many many algorithms algorithms to get get model model parameters, see, if UAVs mass is changing it is needed to build new control There are many algorithms to get model parameters, see, (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure reconfigure controllers parameters. There are many algorithms to get model parameters, see, law or controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the law or reconfigure controllers parameters. e.g. (Kalouptsidis and Theodoridis, 1993). There is the interest to use algorithm, which measures few information There are many applications when quadrocopters are e.g. (Kalouptsidis and Theodoridis, 1993). There is the interest to use algorithm, which measures few information There are many applications when quadrocopters are interest to use algorithm, which measures few information interest to use algorithm, which measures few information There are many applications when quadrocopters are as possible, because we can not know, for example, full There are many applications when quadrocopters are needed: ultralight UAVs can perform many different tasks interest to use algorithm, which measures few information possible, because we can not know, for example, full There are many UAVs applications whenmany quadrocopters are as needed: ultralight can perform different tasks as possible, possible, because we range can not not know, for example, example, full as because we can know, for full needed: ultralight can many different tasks state of plant in aa wide of tasks. Therefore we use needed: ultralight UAVs UAVs can perform perform many different tasks from photography the terrain to carry out flight in the as possible, because we range can not know, for example, full state of plant in wide of tasks. Therefore we use needed: ultralight UAVs can perform many from photography photography the terrain terrain to carry carry outdifferent flight in intasks the identification state of plant in a wide range of tasks. Therefore we use state of plant in a wide range of tasks. Therefore we use from the to out flight the algorithm described in (Gawthrop, 1987) befrom photography the terrain to carry out flight in the mountains. The main condition of all this tasks is stable state of plant in a wide range of tasks. Therefore we identification algorithm described in (Gawthrop, 1987) befrom photography thecondition terrain toofcarry outtasks flightis in the identification algorithm described in (Gawthrop, 1987) use mountains. The main all this stable beidentification algorithm described in (Gawthrop, 1987) bemountains. The main condition of all this tasks is stable cause it is based on knowledge about input-output signals mountains. The main condition of all this tasks is stable and precise altitude holding. When aircraft is flying at an identification algorithm described in (Gawthrop, 1987) because it is based on knowledge about input-output signals mountains. The main condition of all this tasks is stable and precise altitude holding. When aircraft is flying at an cause it is based on knowledge about input-output signals cause it is based on knowledge about input-output signals and precise altitude holding. When aircraft is flying at an and about orders of transfer function only. Let us abbreand precise altitude holding. When aircraft is flying at an altitude of several thousand meters above terrain without cause it is based on knowledge about input-output signals about orders of transfer function only. Let us abbreand precise altitudethousand holding. meters When aircraft is flying at an and altitude of several above terrain without and about about orders of of transfer function only. only. Let us abbreabbreand orders transfer function us altitude of thousand meters terrain viate this method as CIA (Continuous-time Identification altitude of several several thousand meters above abovebecause terrain without without mountains there is no such requirement error of and about orders of transfer function only. Let Let us abbreviate this method as CIA (Continuous-time Identification altitude of several meters above terrain without mountains there is isthousand no such such requirement requirement because error of of Algorithm). viate this method as CIA (Continuous-time Identification viate this method as CIA (Continuous-time Identification mountains there no because error It is easy to see that it had a good approach mountains thereis is is no no such such requirement because error of Algorithm). several meters for this But when viate this method as CIA (Continuous-time Identification It is easy to see that it had aa good approach mountains there because error of several meters meters is not not critical criticalrequirement for this this altitude. altitude. But when Algorithm). It is easy to see that it had approach Algorithm). It is easy to see that it had aa good good approach several is for But when to identification for “LAAS Helicopter Benchmark” with several meters is not not critical critical for over this altitude. altitude. But when aircrafts flight trajectory lies low the ground error in Algorithm). It is easy to see that it had good approach to identification for “LAAS Helicopter Benchmark” with several meters is not critical for this altitude. But when aircrafts flight trajectory lies low over the ground error in to identification for “LAAS Helicopter Benchmark” with to identification for “LAAS Helicopter Benchmark” with aircrafts flight trajectory lies low over the ground error in a good quality, see, e.g. (Fradkov et al., 2010). This alaircrafts flight trajectory lies low over the ground error in several meters may cause a crash. to identification for “LAAS Helicopter Benchmark” with quality, see, e.g. (Fradkov et al., 2010). This alaircrafts flight trajectory lies low over the ground error in aa good several meters may cause a crash. good quality, see, e.g. (Fradkov et al., 2010). This ala good quality, see, e.g. (Fradkov et al., 2010). This alseveral meters may cause a crash. gorithm is described in Section 2. It would be helpful to several meters may cause a crash. a good quality, see, e.g. (Fradkov et al., 2010). This algorithm is described in Section 2. It would be helpful to several meters may cause a crash. gorithm is described in Section 2. It would be helpful to Identification problem can be solved by two ways: ongorithm is described in Section 2. It would be helpful to compare results with another way to identify plant model Identification problem can be solved by two ways: ongorithm is described in Section 2. It would be helpful to compare results with another way to identify plant model Identification problem can be solved by two ways: onIdentification problem can be solved by two ways: oncompare results with another way to identify plant model board in real time and after receiving all data results. compare results with another way to identify plant model Identification problem can be receiving solved byalltwo ways: on- parameters with using non-recursive least squares identiboard in real time and after data results. compare results with another way to identify plant model parameters with using non-recursive least squares identiboard in real time and after receiving all data results. board in real real by timefirst andway after receiving all much data results. results. parameters with using using non-recursive least squares identiidentiIdentification is assumed in papers parameters with non-recursive least squares board in time and after all data algorithm which is described in subsection 3.2 Identification by first way is receiving assumed in in much papers fication parameters with using non-recursive least squares identialgorithm which is described in subsection 3.2 Identification by way is papers Identification bytofirst first waydata is assumed assumed in much much papers fication fication algorithm which is described in subsection 3.2 but it is easier collect after flying and simulate fication algorithm which is described in subsection 3.2 Identification by first way is assumed in much papers but it is easier to collect data after flying and simulate fication algorithm which is described in subsection 3.2 but it is easier to collect data after flying and simulate Consider the following single-input–single-output (SISO) but it is easier to collect data after flying and simulate this identification process separately (how we make in this Consider the following single-input–single-output (SISO) but it is easier to collect data after flying and simulate this identification identification process process separately separately (how (how we we make make in in this this plant Consider the following single-input–single-output (SISO) Consider the following single-input–single-output (SISO) this model: this identification process separately (howStanculeanu we make make in in and this plant paper). Papers (Abas et al., 2011, 2012; Consider the following single-input–single-output (SISO) this identification process separately (how we this paper). Papers (Abas (Abas et al., al., 2011, 2012; 2012; Stanculeanu and plant model: model: plant model: paper). Papers et 2011, Stanculeanu and paper). Papers (Abas et al., 2011, 2012; Stanculeanu and Borangiu, 2011) are devoted to second way. plant model: paper). Papers et al., 2011, 2012;way. Stanculeanu and Borangiu, 2011)(Abas are devoted to second Borangiu, (n−1) Borangiu, 2011) 2011) are are devoted devoted to to second second way. way. yy (n) Borangiu, are devoted to second way. (n) (t) + a1 y (n−1) (t) + . . . + an y(t) = The work 2011) (n) (n−1) (t) + . . . + an y(t) = + a 1y (n) (t) was performed in IPME RAS, supported by RSF (grant y (t) + a y (t) + . .. .. + a 1 y (n−1) (n−1) yy (n)(m) (t) + a (t) + ann y(t) y(t) = = The work was performed in IPME RAS, supported by RSF (grant 1 (m−1) (t) + a y (t) (t) + ..+ . ..+ + = The work was performed in IPME RAS, supported by RSF (grant 1 n y(t) u (t) + b .. .. a+ bm u(t) (1) b 14-29-00142). The work was performed in IPME RAS, supported by RSF (grant (m) (m−1) The work was performed in IPME RAS, supported by RSF (grant 0 u(m) (t) + b1 u (m−1) u (t) + . + (1) b 14-29-00142). 0 u(m) (t) + b1 u(m−1) (t) + . . . + b m u(t) b u(t) (1) b 14-29-00142). 0 u(m) (t) + b1 1 u(m−1) (t) + . . . + bm m u(t) (1) b 14-29-00142). 0 (t) + . . . + b u(t) (1) b u (t) + b u 14-29-00142).
Konstantin Konstantin Konstantin Konstantin Konstantin
0
1
Copyright © 2015, IFAC IFAC 2015 (International Federation of Automatic Control) 237Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © Copyright © IFAC 2015 237 Copyright ©under IFAC responsibility 2015 237Control. Copyright 237 Peer review© of International Federation of Automatic Copyright © IFAC IFAC 2015 2015 237 10.1016/j.ifacol.2015.09.189
m
MICNON 2015 234 Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237 June 24-26, 2015. Saint Petersburg, Russia
where ai , i = 1, 2, . . . , n and bj , j = 1, 2, . . . , m are unknown parameters which are needed to be identified. y(t), u(t) is output and input respectively. Let us introduce two state filters ζ˙ = Ad ζ + bd y(t), σ˙ = Ad σ + bd u(t)
(2) where ζ(t), σ(t) ∈ R ; pair (Ad , bd ) has a regular canonical form and det(sIn −Dd ) = D(s) is an arbitrary Hurwitz polynomial D(s) = sn + d1 sn−1 + . . . + dn . N
Identification algorithm has the following form ˜ ˙ ρ(t)ζ(t) Ω(t) = −Γ(t)˜ ρ(t)˜ ρ(t)T Ω(t) + Γ(t)˜ ˙ Γ(t) = −Γ(t)˜ ρ(t)˜ ρ(t)T Γ(t) + αΓ(t)
(3) (4)
n+m+1
where Ω(t) ∈ R – vector of estimates for ai , i = 1, 2, . . . , n and bj , j = 1, 2, . . . , m; Γ(t) is gain matrix; ρ˜ = [ζn , . . . , ζ1 , σm+1 , . . . , σ1 ] is vector of filters (2) n outputs; ζ˜ = y(t) − dn−i+1 ζi and α > 0 is algorithm i=1
parameter. It is needed to take initial value Γ(0) = Γ(0)T and this parameter is responsible for identification speed. 3. QUADROTOR MODEL 3.1 Design modeling dynamics
For quadrocopters the problem of altitude hold is simple but requires precision and reliability. The quadrotor helicopter dynamics model may be derived using EulerLagrange formalism under the assumptions that the frame structure is rigid, the structure is symmetrical, the centre of gravity (CoG) and the body fixed frame origin are assumed to coincide and the propellers are rigid. There is significant amount of literature dedicated to quadrotors, where the model of platform dynamics is derived, see papers by Fantoni and Lozano (2002); Altug et al. (2002); Castillo et al. (2005); Das et al. (2008); Janusz et al. (2013) for mentioning a few. This model has the following form: Ix p˙ = (Iy − Iz )qr + Tϕ (t), Iy q˙ = (Iz − Ix )pr + Tθ (t), Iz r˙ = (Ix − Iy )pq + Tψ (t), (5) m¨ x = (sin ψ · sin ϕ + cos ψ · sin θ · cos ϕ)F (t), m¨ y = (sin ψ · sin θ cos ϕ − cos ψ · sin ϕ)F (t), m¨ z = cos θ · cos ϕ · F (t) − mg,
where ϕ, θ, ψ are the Euler angles (roll, pitch, yaw, respectively); p, q, r – angular velocities; (x, y, z) gives the quadrotor CoG position in the normal Earth reference frame, Ix , Iy , Iz stand for inertia moments; Tϕ (t), Tθ (t), Tψ (t) are external torques, acting to the quadrotor; m denotes the quadrotor mass; g is the gravity acceleration. F (t) denotes the following total thrust, produced by the rotors: F (t) = f1 (t) + f2 (t) + f3 (t) + f4 (t), (6) where fi (t), i = 1, . . . 4, are forces, applied from each one rotor. It is worth mentioning that fi (t) depend on the blades rotation speed, aerodynamic drag and external disturbances. 238
In the present paper we focus our attention to quadrotor movement in height, assuming that the quadrotor is horizontally positioned. Therefore angles ϕ, θ in (5) should be taken zero, yaw angle ψ may have an arbitrary value. This leads to the following equation for isolated vertical motion dynamics: m¨ z = F (t) − mg, (7) where z denotes the quadrotor altitude. Interaction between rotating blades and surrounding air has a complex nature, cf. (Leishman, 2006; Janusz et al., 2013). Based on (Janusz et al., 2013) we assume that this interaction may be modeled as a viscous friction of the vertical motion. Then the vertical motion dynamics may be described as follows: 4 ˙ = Ku Ui (t) − mg, (8) m¨ z (t) + Kz˙ z(t) i=1
where Ui denotes the rotation speed of ith propeller (i = 1, . . . , 4), Kz˙ is the viscous friction parameter, Ku may be called the propeller effectiveness gain.
In what follows, due to the assumption that the quadrotor is positioned horizontally, it is taken that U1 (t) = . . . U4 (t). Provisionally neglecting the brushless motors dynamics (which are fast in compare with the quadrotor body dynamics), we obtain the following model of the vertical motion dynamics z¨(t) + K z(t) ˙ = Ku Uz (t) − g, (9) 4 where Uz = Ui is the altitude control signal (the sum i=1
of PWM signals Ui , applied to the motors), K, Ku are normalized model parameters. K, Ku depend on various factors and they are difficult to calculate. Therefore we use an identification procedure to find them based on the experimental results
In this way altitude hold is needed in working of servos: their generated traction must compensate gravity. To make this it is available to use PID-controller as example. 3.2 Estimate by LSE Coordinates z, acceleration z¨ and control signal Uz were obtained according to the results of the flight from the autopilot dataflash. With the help of differentiation by Eulers method the velocities and the accelerations vector z¨ can be obtain on the basis of the vector z or z˙ can be obtained by integration z¨. Moreover, the production of the vector z˙ by differentiation practically coincides with the data obtained from the quadrotors accelerators. The equation describing the dynamics (9) can be rewritten in matrix form χ = Y, (10) where Y = row(¨ z (1), z¨(2), . . . , z¨(n)), χ has the form Uz (1) 1 z(1) ˙ ˙ Uz (2) 1 z(2) , (11) χ= .. .. ... . . Uz (n) 1 z(n) ˙ where n – number of measurements. – vector of unknown parameters. From the expression (10) vector of unknown parameters is found as = χ+ Y , where χ+ – pseudoinverse matrix to matrix χ.
MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237
235
Let us identify parameters of quadrotor with model (9). In this subsection quadrotor model was found only for test No. 1 to show a rough estimate of the identification parameters (parameters of quadrotor for this case can be seen in Tab. 1). Result of identification by this method for other experiments is shown in Tab. 2. Following parameters have been found: Ku = 0.0208, g = 9.5841, K = 0.4631, then the dynamic equation (9) takes the form (12) z¨ + 0.4631z˙ = 0.0208Uz − 9.5841. From (12) it can be seen that the obtained acceleration of gravity close to how it should be gr ≈ 9.81, so we can assume it g = 9.81. 4. THE PRACTICAL APPLICATION FOR QUADROCOPTER For experiments four configuration are used, all variants is shown in Table 1 Table 1. Configurations parameters of quadrocopter. Test No. 1 2 3 4
Battery, mAh 4000 2600 2600 4000
Prop, inch 10 × 4.5 10 × 4.5 8×6 8×6
Mass, g 1220 1120 1118 1218
Fig. 3. Velocity z˙ of quadrotor. where Kz Uz = F1 + F2 + F3 + F4 and Kz is conversion factor. This equation is given for moments in time, when z¨ = 0, z˙ = 0 and ϕa = ψa = 0 because we consider the movement in a vertical plane. The equation (13) allows us to estimate conversion factor Kz = mg Uz ≈ 0.018. This estimation coincides with that obtained by the least squares method in the model (12). In the light of the above, equation (9) takes the form (14) m¨ z + K z˙ = 0.018Uz − m9.81. Knowledge of the equation (14) simplify identification, because this is necessary to evaluate one parameter less. But it will be shown below that this parameter Kz can be estimate too, but with some problems. 5. RESULTS OF IDENTIFICATION
Altitude is taken as the output value. These measurements are red from GPS and barometer, then mean value is used. The input signal is sum of PWM signals which are sending on each of motors. This signal is nondimensional quantity: it takes values from 1118 to 1919 for each motor and converts in autopilot with the following law: at the minimum (1118) rotation speed is zero and at the maximum (1919) rotation speed is 920 rmp/V × 12V , linear dependence. Input signal for one case is presented on Figure 1, acceleration z¨ and velocity z˙ are presented on Fig. 2 and on Fig 3 respectively.
Using identification algorithm which is described in Section 2, for four variations of equipment gives us result which is shown in the Table 2. It should be noted too that gravitational constant is known and this knowledge can simplify the identification because we need to estimate only one parameter: coefficient K. Since we use the method which is described in Section 2 it is needed to choose parameters for identification algorithm. Let us use next settings: α = 0.001, Γ(0) = 1000, ρ˜(0) = col{0, 0, 0}. Since there was four different experiments with various configurations and it means that it is necessary to allocate equal time intervals and identify parameters for these intervals to compare setting process later. These time histories is presented on Figs. 4-6.
Fig. 1. Input signal Uz . Fig. 4. Time histories of identification parameter g. Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively.
Fig. 2. Acceleration z¨ of quadrotor. As we can see, when input signal Uz ≈ 5800, then z¨ = 0, z˙ = 0. It means that equation (9) can be represented as follows Ku Uz = mg (13) 239
It can be seen that the parameters identifications for all four cases are the same and almost equal. There are some problems with parameter K: identification has a small variation and does not converge as the two other parameters, but it can be estimated separately as it is stated above. Finally, it is needed to show how the identified model reproduces input signal for one of the cases (test No 2).
MICNON 2015 236 Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237 June 24-26, 2015. Saint Petersburg, Russia
by LSE. Quadrocopters model dynamics for motion in a vertical plane are obtained. The results obtained are compared with the theoretical statements and efficiency of the resulting model is demonstrated. The difficulties encountered in the identification are noted.
Fig. 5. Time histories of identification parameter K. Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively.
The future work intentions are decide difficulties with equipment and make possible to have more accurately onboard identification. Also it is planned to identify other parameters quadrocopter to receive full model. REFERENCES
Fig. 6. Time histories of identification parameter Ku . Solid, dashed, dotted and dash-dotted lines correspond tests No 1, 2, 3, 4 respectively. Let us use model (15) to show efficiency of the resulting model (15) z¨ + 0.804z˙ = 0.016Uz − 9.744. Result of representing can be seen on Fig. 7. There is a good match with modeling output signal with the resulting model and real output.
Fig. 7. Output signal z, estimate z end error – blue solid line, black solid line and dashed line respectively. In the Figure 7 solid blue line is quadrocopter real altitude which was taken on certain time interval for 40 seconds, solid blue line is estimation of the altitude, obtained using the identified model and dashed line is error e(t) = z(t) − z(t). There are some problems with identification: input signals were very noisy at some time intervals and sensors accumulated errors over time, making it difficult to identify UAVs dynamic at large intervals of time. This is the reason why problems with identification of parameter K were found. 6. CONCLUSION Parameter identification is studied experimentally by two methods: by continuous-time identification algorithm and Table 2. The estimates obtained by two methods Test No. 1 2 3 4
g, m/s2 CIA LSE 9.638 8.516 9.071 7.945 11.213 9.972 11.714 8.313
K, 1/s CIA LSE 0.424 0.758 0.301 0.704 0.406 0.456 0.863 0.848
Ku CIA LSE 0.015 0.014 0.012 0.013 0.016 0.014 0.017 0.013
240
Abas, N., Legowo, A., and Akmeliawati, R. (2011). Parameter identification of an autonomous quadrotor. In 2011 4th International Conference On Mechatronics (ICOM), 1–8. Abas, N., Legowo, A., Ibrahim, Z., Rahim, N., and Kassim, A.M. (2012). Modeling and system identification using extended kalman filter for a quadrotor system. In Proceedings of 2012 International Conference on Electrical Engineering and Applications (ICEEA 2012). Altug, E., Ostrowski, J., and Mahony, R. (2002). Control of a quadrotor helicopter using visual feedback. In Proc. IEEE Int. Conf. on Robotics and Automation, 72–76. IEEE Press, Washington, DC, USA. Amelin, K., Amelina, N., Granichin, O., Granichina, O., and Andrievsky, B. (2013). Randomized algorithm for uavs group flight optimization. In Proc. of 11th IFAC International Workshop on Adaptation and Learning in Control and Signal Processing, 205–208. Castillo, G.P., Lozano, R., and Dzul, A.E. (2005). Modelling and Control of Mini-Flying Machines. Series: Advances in Industrial Control. Springer-Verlag Ltd. Das, A., Subbarao, K., and Lewis, F. (2008). Dynamic inversion of quadrotor with zero-dynamics stabilization. In Proc. 17th IEEE Int. Conf. on Control Applications, Part of 2008 IEEE Multi-conference on Systems and Control (MSC 2008). IEEE Press, San Antonio, USA. Fantoni, I. and Lozano, R. (2002). Non-linear Control for Underactuated Mechanical Systems. Series: Communications and Control Engineering. Springer-Verlag Ltd., London. Fradkov, A.L., Andrievsky, B., and Peaucelle, D. (2010). Estimation and control under information constraints for laas helicopter benchmark. Control Systems Technology, IEEE Transactions on, 18(5), 1180–1187. Gawthrop, P. (1987). Continuous-Time Self-Tuning Control. Letchworth, U.K.: Research Studies Press. Janusz, W.C., Czyba, R., and Szafra´ nski, G. (2013). Model identification and data fusion for the purpose of the altitude control of the VTOL aerial robot. In Proc. 2nd IFAC Workshop on Research, Education and Development of Unmanned Aerial Systems (RED-UAS 2013), 263–269. IFAC, Compiegne, France. Kalouptsidis, N. and Theodoridis, S. (eds.) (1993). Adaptive system identification and signal processing algorithms. Prentice-Hall, Inc., Upper Saddle River, NJ, USA. Leishman, J. (2006). Principles of Helicopter Aerodynamics. 2nd ed. Cambridge University Press, New York. Ljung, L. (1999). System Identification: Theory for the User. Prentice-Hall information and system sciences series. PTR Prentice Hall.
MICNON 2015 June 24-26, 2015. Saint Petersburg, Russia Konstantin Amelin et al. / IFAC-PapersOnLine 48-11 (2015) 233–237
Stanculeanu, I. and Borangiu, T. (2011). Quadrotor blackbox system identification. In World Academy of Science, Engineering and Technology, volume 5, 276–279. Appendix A. TECHNICAL CHARACTERISTICS We use the ultralight UAV which is described in (Amelin et al., 2013). Autopilot Ardupilot Mega 2.6 with inertial system (gyroscope, accelerometer, magnetometer, barometer) and receiver signal from the satellite navigation system GPS. Autopilot Ardupilot Mega 2.6 is built on chips Atmels ATMEGA 2560, has 8 inputs and outputs for the PWM signal, three UART port, a built-in FTDI, 8 analog inputs, 4 Mb Flash memory for recording log data telemetry. It is equipped with Invensenses 6 DoF Accelerometer/Gyro MPU-6000, a barometer MS5611-01BA03 and 4 sonars MB1260 XL-MaxSonar-EZL0. Corps UAV – quadrocopter frame DJI 450. Driving body – radial-beam with the size of 450 mm in diameter. At the end of each of the four beams is a brushless motor with an external rotor and the turns controller (control mechanisms). The main UAV characteristics are as follows: • • • • •
The length of the beams, mm: 225; Flight weight, gram: varied in experiments; Flight speed, km/h: from 0 to 10; Engines, watts: 250 each of the four engines; Battery, mAh: 2600 or 4000 (depending on the experiment); • Autopilot and additional microcomputer: installed.
Fig. A.1. Photo of quadrocopter, which was used for the tests.
241
237