Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements

ScienceDirect IFAC-PapersOnLine 48-9 (2015) 067–072 ScienceDirect Available online at www.sciencedirect.com 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Peer review under responsibility of International Federation of Automatic Control. 10.1016/j.ifacol.2015.08.061 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements  Ram´on A. Garc´ıa ∗Guilherme V. Raffo ∗∗ Manuel G. Ortega ∗ Francisco R. Rubio ∗ ∗Depto. Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Spain (e-mail: {ramongr, mortega, rubio}@us.es) ∗∗ Depto. Engenharia Eletrˆonica, Universidade Federal de Minas Gerais, Brazil (e-mail: [email protected]) Abstract: The application of two guaranteed estimation methods on the estimation of a quadrotor position is presented in this paper. A discrete time system model used for the position estimation is supposed, whereas GPS measurement is performed with greater sampling time. Firstly, guaranteed algorithms are applied in order to compute the feasible set of positions where the quadrotor is likely to be found. Subsequently, the estimation is corrected and improved by using the GPS sensor measurement. Keywords: Guaranteed estimation, intervalar arithmetic, zonotopes, quadrotor, GPS. 1. INTRODUCTION A quadrotor is an aerial vehicle which has four coplanar motors. The forces and torques that are applied over this kind of vehicles are shown in Fig. 1. Fig. 1. Sketch of main forces and torques applied on a quadrotor. It is a continuous system whose position measurement in outdoors environment is usually performed by a GPS in discrete mode. In case that the UAV position would like to The authors would like to thank the MCeI for funding this work through projects DPI2010-19154 and DPI2012-37580-C02-02, as well as FAPEMIG and Programa Institucional de Aux´ılio `a Pesquisa de Doutores Rec´emContratados of the PRPq/UFMG. be known at times lesser than the GPS sample period, a discrete observer could be used for its estimation. Thereby, two sampling times can be assumed, one due to discrete observer implementation, T, and the other due to GPS measurements, tsynchro, being T <t synchro. The proposed time frame is represented in Fig.2: at starting time t1, the estimation of the position through the equation model as well as the GPS signal measurement are available. Until tsynchro seconds later, the GPS measurement is not available again. Meanwhile, in the intermediate sample times, every Tseconds, the estimation based on the equation of system model is used. Fig. 2. Different sampling times for the position estimation. It is also interesting how to deal with the possible uncertainties, since the models have diverse uncertainties which may come from various sources, such as unmodelled system dynamics, inability to measure disturbances such as wind, ground effect, etc. Also it is likely to have uncertainProceedings of the 2015 IFAC Workshop on Advanced Control and Navigation for Autonomous Aerospace Vehicles June 10-12, 2015. Seville, Spain Copyright © 2015 IFAC 67 Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements  Ram´on A. Garc´ıa ∗Guilherme V. Raffo ∗∗ Manuel G. Ortega ∗ Francisco R. Rubio ∗ ∗Depto. Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Spain (e-mail: {ramongr, mortega, rubio}@us.es) ∗∗ Depto. Engenharia Eletrˆonica, Universidade Federal de Minas Gerais, Brazil (e-mail: r[email protected]) Abstract: The application of two guaranteed estimation methods on the estimation of a quadrotor position is presented in this paper. A discrete time system model used for the position estimation is supposed, whereas GPS measurement is performed with greater sampling time. Firstly, guaranteed algorithms are applied in order to compute the feasible set of positions where the quadrotor is likely to be found. Subsequently, the estimation is corrected and improved by using the GPS sensor measurement. Keywords: Guaranteed estimation, intervalar arithmetic, zonotopes, quadrotor, GPS. 1. INTRODUCTION A quadrotor is an aerial vehicle which has four coplanar motors. The forces and torques that are applied over this kind of vehicles are shown in Fig. 1. Fig. 1. Sketch of main forces and torques applied on a quadrotor. It is a continuous system whose position measurement in outdoors environment is usually performed by a GPS in discrete mode. In case that the UAV position would like to The authors would like to thank the MCeI for funding this work through projects DPI2010-19154 and DPI2012-37580-C02-02, as well as FAPEMIG and Programa Institucional de Aux´ılio `a Pesquisa de Doutores Rec´emContratados of the PRPq/UFMG. be known at times lesser than the GPS sample period, a discrete observer could be used for its estimation. Thereby, two sampling times can be assumed, one due to discrete observer implementation, T, and the other due to GPS measurements, tsynchro, being T <t synchro. The proposed time frame is represented in Fig.2: at starting time t1, the estimation of the position through the equation model as well as the GPS signal measurement are available. Until tsynchro seconds later, the GPS measurement is not available again. Meanwhile, in the intermediate sample times, every Tseconds, the estimation based on the equation of system model is used. Fig. 2. Different sampling times for the position estimation. It is also interesting how to deal with the possible uncertainties, since the models have diverse uncertainties which may come from various sources, such as unmodelled system dynamics, inability to measure disturbances such as wind, ground effect, etc. Also it is likely to have uncertainProceedings of the 2015 IFAC Workshop on Advanced Control and Navigation for Autonomous Aerospace Vehicles June 10-12, 2015. Seville, Spain Copyright © 2015 IFAC 67 Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements  Ram´on A. Garc´ıa ∗Guilherme V. Raffo ∗∗ Manuel G. Ortega ∗ Francisco R. Rubio ∗ ∗Depto. Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Spain (e-mail: {ramongr, mortega, rubio}@us.es) ∗∗ Depto. Engenharia Eletrˆonica, Universidade Federal de Minas Gerais, Brazil (e-mail: [email protected]) Abstract: The application of two guaranteed estimation methods on the estimation of a quadrotor position is presented in this paper. A discrete time system model used for the position estimation is supposed, whereas GPS measurement is performed with greater sampling time. Firstly, guaranteed algorithms are applied in order to compute the feasible set of positions where the quadrotor is likely to be found. Subsequently, the estimation is corrected and improved by using the GPS sensor measurement. Keywords: Guaranteed estimation, intervalar arithmetic, zonotopes, quadrotor, GPS. 1. INTRODUCTION A quadrotor is an aerial vehicle which has four coplanar motors. The forces and torques that are applied over this kind of vehicles are shown in Fig. 1. Fig. 1. Sketch of main forces and torques applied on a quadrotor. It is a continuous system whose position measurement in outdoors environment is usually performed by a GPS in discrete mode. In case that the UAV position would like to The authors would like to thank the MCeI for funding this work through projects DPI2010-19154 and DPI2012-37580-C02-02, as well as FAPEMIG and Programa Institucional de Aux´ılio `a Pesquisa de Doutores Rec´emContratados of the PRPq/UFMG. be known at times lesser than the GPS sample period, a discrete observer could be used for its estimation. Thereby, two sampling times can be assumed, one due to discrete observer implementation, T, and the other due to GPS measurements, tsynchro, being T <t synchro. The proposed time frame is represented in Fig.2: at starting time t1, the estimation of the position through the equation model as well as the GPS signal measurement are available. Until tsynchro seconds later, the GPS measurement is not available again. Meanwhile, in the intermediate sample times, every Tseconds, the estimation based on the equation of system model is used. Fig. 2. Different sampling times for the position estimation. It is also interesting how to deal with the possible uncertainties, since the models have diverse uncertainties which may come from various sources, such as unmodelled system dynamics, inability to measure disturbances such as wind, ground effect, etc. Also it is likely to have uncertainProceedings of the 2015 IFAC Workshop on Advanced Control and Navigation for Autonomous Aerospace Vehicles June 10-12, 2015. Seville, Spain Copyright © 2015 IFAC 67 Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements  Ram´on A. Garc´ıa ∗Guilherme V. Raffo ∗∗ Manuel G. Ortega ∗ Francisco R. Rubio ∗ ∗Depto. Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Spain (e-mail: {ramongr, mortega, rubio}@us.es) ∗∗ Depto. Engenharia Eletrˆonica, Universidade Federal de Minas Gerais, Brazil (e-mail: [email protected]) Abstract: The application of two guaranteed estimation methods on the estimation of a quadrotor position is presented in this paper. A discrete time system model used for the position estimation is supposed, whereas GPS measurement is performed with greater sampling time. Firstly, guaranteed algorithms are applied in order to compute the feasible set of positions where the quadrotor is likely to be found. Subsequently, the estimation is corrected and improved by using the GPS sensor measurement. Keywords: Guaranteed estimation, intervalar arithmetic, zonotopes, quadrotor, GPS. 1. INTRODUCTION A quadrotor is an aerial vehicle which has four coplanar motors. The forces and torques that are applied over this kind of vehicles are shown in Fig. 1. Fig. 1. Sketch of main forces and torques applied on a quadrotor. It is a continuous system whose position measurement in outdoors environment is usually performed by a GPS in discrete mode. In case that the UAV position would like to The authors would like to thank the MCeI for funding this work through projects DPI2010-19154 and DPI2012-37580-C02-02, as well as FAPEMIG and Programa Institucional de Aux´ılio `a Pesquisa de Doutores Rec´emContratados of the PRPq/UFMG. be known at times lesser than the GPS sample period, a discrete observer could be used for its estimation. Thereby, two sampling times can be assumed, one due to discrete observer implementation, T, and the other due to GPS measurements, tsynchro, being T <t synchro. The proposed time frame is represented in Fig.2: at starting time t1, the estimation of the position through the equation model as well as the GPS signal measurement are available. Until tsynchro seconds later, the GPS measurement is not available again. Meanwhile, in the intermediate sample times, every Tseconds, the estimation based on the equation of system model is used. Fig. 2. Different sampling times for the position estimation. It is also interesting how to deal with the possible uncertainties, since the models have diverse uncertainties which may come from various sources, such as unmodelled system dynamics, inability to measure disturbances such as wind, ground effect, etc. Also it is likely to have uncertainProceedings of the 2015 IFAC Workshop on Advanced Control and Navigation for Autonomous Aerospace Vehicles June 10-12, 2015. Seville, Spain Copyright © 2015 IFAC 67 68 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 ties due to sensor measurements, whose problem is shown in Fig.3. Fig.3a shows the case wherein a completely accurate measurement is considered, while the case in which an inaccurate sensor is assumed is depicted in Fig.3b. More in detail, the shadowed zones in these figures represent all the possible vehicle trajectories provided by the observer, while the points at the synchronization instants are the measure provided by the GPS sensor; finally, the dashed lines represent the actual vehicle trajectories. It can be observed in Fig.3a how the estimation is accurately updated with the refreshed GPS measure. However, Fig. 3b depicts how the measurement points do not coincide with the actual trajectory at those synchronization instants, due to sensor inaccuracy. (a) (b) Fig. 3. GPS measurement examples. A lot of possible methods for the position estimation of autonomous vehicles can be found in the literature (see, for example, Abbott and Powell (1999), Ra¨ısi et al. (2012), Farrelly and Wellstead (1996), Scholte and Campbell (2002) and Brunke and Campbell (2002), among others). Based in the classification carried out in Combastel (2003), three main approaches may be stood out: (1) The Luenberger observers approach, where uncertainties are not explicitly considered. (2) Estimators based on the statistical behavior, such as Kalman filter, that presupposes noise with some probability function. (3) Estimators based on guaranteed techniques, which do not assume any probability function. They are based on the assumption of bounded variables within a certain range, i.e., the variables behaves as compact sets (in Rnare bounded and closed). The use of some methodologies, such as interval arithmetics, can be considered as example of the latter type, being the result of state estimation a compact set in the state space representing an outer approximation of all the states that are consistent both with the uncertain model and the uncertain measurements. Thereby, a statistic study of the noise is not necessary, which is some times hard to compute. Returning to the general case, that kind of observers provide a compact set whose domain can be represented by many ways, such as ellipsoids, box (due to using interval arithmetic), parallelotopes or even limited complexity polytopes, i.e., with a limited number of vertexes and faces. According to the method used to represent the domain, the wrapping 1effect is appreciated in greater or lesser way, i.e., the method is more or less conservative in terms of how it delimits the feasible set. It is not the same to enclose it by a box with sides parallel to the axes, than by a parallelepiped or by a zonotope. In Fig.4, the set of possible states as well as two external bounded boxes can be observed: the black one, which consists in a box with sides parallel to the axes (assuming horizontal and vertical axes), and the green one, characterized by a parallelepiped whose sides are not parallel to the axes. It can be seen how the last one tries to adapt in a better way to the shape of the possible state set. The solution obtained by the first box is more conservative than the one obtained by the parallelepiped. Fig. 4. Example of fitting of possible state set. In this work, an observer that estimates the state of a quadrotor UAV (position and linear velocity components) for times between updated GPS measurements is presented. The remainder of the paper is organized as follows. Section 2 presents the model of the system. In Section 3, the algorithms used are introduced. Afterwards, in Section 4, simulation results are presented. And finally, the main conclusions are drawn in Section 5. 2. SYSTEM MODELING In general, let us consider a discrete time system with uncertainties both in the model and the measurement: xk+1 =f(xk,w k) yk=g(xk,v k)(1) where xk∈Rnis the state system vector, yk∈Rpis the output vector, wk∈Rnwrepresents modeling errors and disturbances and finally, vk∈Rpvrepresents the noise measurements vector. 1The wrapping effect is the fact to add more solutions than the real one. It is due to use a method to bound so conservative, and therefore a possible state set bigger is obtained than the real state set that the system could achieve. IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 68 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 69 ties due to sensor measurements, whose problem is shown in Fig.3. Fig.3a shows the case wherein a completely accurate measurement is considered, while the case in which an inaccurate sensor is assumed is depicted in Fig.3b. More in detail, the shadowed zones in these figures represent all the possible vehicle trajectories provided by the observer, while the points at the synchronization instants are the measure provided by the GPS sensor; finally, the dashed lines represent the actual vehicle trajectories. It can be observed in Fig.3a how the estimation is accurately updated with the refreshed GPS measure. However, Fig. 3b depicts how the measurement points do not coincide with the actual trajectory at those synchronization instants, due to sensor inaccuracy. (a) (b) Fig. 3. GPS measurement examples. A lot of possible methods for the position estimation of autonomous vehicles can be found in the literature (see, for example, Abbott and Powell (1999), Ra¨ısi et al. (2012), Farrelly and Wellstead (1996), Scholte and Campbell (2002) and Brunke and Campbell (2002), among others). Based in the classification carried out in Combastel (2003), three main approaches may be stood out: (1) The Luenberger observers approach, where uncertainties are not explicitly considered. (2) Estimators based on the statistical behavior, such as Kalman filter, that presupposes noise with some probability function. (3) Estimators based on guaranteed techniques, which do not assume any probability function. They are based on the assumption of bounded variables within a certain range, i.e., the variables behaves as compact sets (in Rnare bounded and closed). The use of some methodologies, such as interval arithmetics, can be considered as example of the latter type, being the result of state estimation a compact set in the state space representing an outer approximation of all the states that are consistent both with the uncertain model and the uncertain measurements. Thereby, a statistic study of the noise is not necessary, which is some times hard to compute. Returning to the general case, that kind of observers provide a compact set whose domain can be represented by many ways, such as ellipsoids, box (due to using interval arithmetic), parallelotopes or even limited complexity polytopes, i.e., with a limited number of vertexes and faces. According to the method used to represent the domain, the wrapping 1effect is appreciated in greater or lesser way, i.e., the method is more or less conservative in terms of how it delimits the feasible set. It is not the same to enclose it by a box with sides parallel to the axes, than by a parallelepiped or by a zonotope. In Fig.4, the set of possible states as well as two external bounded boxes can be observed: the black one, which consists in a box with sides parallel to the axes (assuming horizontal and vertical axes), and the green one, characterized by a parallelepiped whose sides are not parallel to the axes. It can be seen how the last one tries to adapt in a better way to the shape of the possible state set. The solution obtained by the first box is more conservative than the one obtained by the parallelepiped. Fig. 4. Example of fitting of possible state set. In this work, an observer that estimates the state of a quadrotor UAV (position and linear velocity components) for times between updated GPS measurements is presented. The remainder of the paper is organized as follows. Section 2 presents the model of the system. In Section 3, the algorithms used are introduced. Afterwards, in Section 4, simulation results are presented. And finally, the main conclusions are drawn in Section 5. 2. SYSTEM MODELING In general, let us consider a discrete time system with uncertainties both in the model and the measurement: xk+1 =f(xk,w k) yk=g(xk,v k)(1) where xk∈Rnis the state system vector, yk∈Rpis the output vector, wk∈Rnwrepresents modeling errors and disturbances and finally, vk∈Rpvrepresents the noise measurements vector. 1The wrapping effect is the fact to add more solutions than the real one. It is due to use a method to bound so conservative, and therefore a possible state set bigger is obtained than the real state set that the system could achieve. IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 68 For the methods to be employed in the resolution of the estimation, the uncertainties are assumed bounded between a certain values, regardless of the statistical distribution. Therefore, wk∈W,vk∈Vand x0∈X0 are assumed. By particularizing (1) in the case that a quadrotor spatial position is desired, starting from Newton’s second law, by considering all forces are grouped in one only component. Being analogous the expressions for the three axes, only the expression for the component xis considered, its expression is ¨x=Fx m. By moving to space state, two states can be defined, x1=xand x2=˙x A system described by the position and the velocity is given in R6:        ˙x1 ˙x2 ˙x3 ˙x4 ˙x5 ˙x6        =              ˙x ˙y ˙z Fx m Fy m Fz m              (2) where x, y, z are the position coordinates, ˙x, ˙y, ˙zare the linear velocities, mis the vehicle mass and Fx,F y,F zare the forces applied into each axis respectively. By expressing (2) in a discrete form:        x1k+1 x2k+1 x3k+1 x4k+1 x5k+1 x6k+1        =        x1k+Tx 4k x2k+Tx 5k x3k+Tx 6k x4k x5k x6k        +                    T2 2 Fx m T2 2 Fy m T2 2 Fz m TFx m TFy m TFz m                    (3) where Tis the sample time. The applied forces are shown in a schematic way, if a quadrotor rotational subsystem is included, the applied forces can be expressed as: Fx= (cos φcos ψsin θ+ sin φsin ψ)u+Ax Fy= (cos φsin ψsin θ−sin φcos ψ)u+Ay Fz= (cos φcos θ)u+Az (4) being φ,θ,ψ the Euler’s angles (roll, pitch and yaw respectively), uthe total thrust and Aithe disturbance in the correspondent axis i. However, in (3) a simplified notation is used where these angles or the thrust are not mentioned, because this task is the goal of the attitude controller that is not considered in this work. Therefore, the forces Fx,F y,F zneeded for the traslational movement are calculated and then, with theses forces the rotational controller calculates the angles and torques references. After introducing the system, the resolution of the problem is tackled. Given the system, its position trough a discrete time model and some measurements is estimated. The model has uncertainties and the measurements could have uncertainties or not, according on the case of study. 3. ALGORITHMS Two cases of study are shown below: interval arithmetic and zonotopes. The movement in the plane XY is considered, the state variables of (3) are realigned and the components corresponding to axis z are not used. 3.1 Intervalar Arithmetic The idea of using interval arithmetic as guaranteed estimation method arises from the knowledge of the existence of limits of the uncertainties. For instance, there are a maximum and minimum mass, but the value is not known exactly because it can change slightly, imagine that a radiocontol (rc) plane with a engine instead of a motor, the amount of gas will change during the flight. Once the application goal has been presented, the realigned system from (3) is given by:    x1k+1 x2k+1 x3k+1 x4k+1   =   x1k+Tx 3k x2k+Tx 4k x3k x4k   +            T2 2 Fx m T2 2 Fy m TFx m TFy m            (5) where x1is the position x,x2is the position y,x3is the velocity ˙xand x4is the velocity ˙y. The rest of the notation is the same as the one used in (3). To deal with the problem, an adaptation of the algorithm presented in Jaulin (2002) is used. This algorithm used interval arithmetic to estimate de state of an autonomous system (xk+1 =f(xk)). Being xa vector, [x] is a vector composed by intervals, it is called a box. f() is the inclusion function of f(). The algorithm can be sketched as follows: •Step 1: [ˆx](t2)=[x](t1)+(t2−t1)f([x](t1)) •Step 2: [v]=[x](t1)∪[ˆx](t2) •Step 3: [w]=inflated([v], αω[v]+β) •Step 4: [x](t2)=[x](t1)+(t2−t1)f([w]) where ω[v] is the length of its greater side and the inflated operation consists in given a box, [x], and a escalar, : inflated([x],)=[x1−, x1+]×... ×[xn−, xn+] As a negative result of using this kind of so conservative method is the wrapping effect. A black box bounding the function fis observed in Fig.4, this box clearly has areas which do not belong to the image of function f. IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 69 70 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 The source of model uncertainty is assumed as the fact that the mass can change, m. Due to this fact, an interval is obtained, in which the different states can be found. The inclusion function 2of (5) is computed as follows (for simplicity the subscripts kare omitted ): f([x]) =            [x1, x1]+T[x3, x3]+T2 2[1 m,1 m]Fx [x2, x2]+T[x4, x4]+T2 2[1 m,1 m]Fy [x3, x3]+T[1 m,1 m]Fx [x4, x4]+T[1 m,1 m]Fy            (6) being Tthe sample time, α=0.1 and β=0.0001. In the calculation of ω[v], states are separated according to coordinate represented, i.e., x1and x3on one side and x2 and x4on the other one. Then it is applied to each pair separately, by using the larger of each pair. 3.2 Zonotopes The idea is to use another procedure that let a guaranteed estimation whose result can be less conservative enclosure. Once again the knowledge of a limit in the uncertainty is the initial assumption. However, this time a rectangular shape or box is not assumed. In this case a more complex shape is allowed, specifically zonotopes. In this approach, the following system obtained from (3) is used:    x1k+1 x2k+1 x3k+1 x4k+1   =   x1k+Tx 3k x2k+Tx 4k x3k x4k   +            T2 2 Fx m T2 2 Fy m TFx m TFy m            (7) where x1is the position x,x2is the position y,x3is the velocity ˙x, and x4is the velocity ˙y. The rest of the notation is the same as the one used in (3). According to the algorithm shown in Bravo (2004), the following steps are considered: •Step 1: Use an inclusion function in order to narrow the trajectory of the uncertainty nonlinear system. •Step 2: Calculate a bound of the consistent state set. •Step 3: Calculate a tight bound of the intersection set. The K¨uhn method K¨uhn (1998) is used for step 1. This method provides a zonotope that encloses the possible system state. A strip for variable set of feasible states is estimated from the measurements obtained in the second step. And finally, in the step 3, intersection between 2According to other authors can be represented as [f] zonotope and strip is calculated. Thereby, how to calculate a strip for a possible states is shown below. The way to obtain a state set through a measurement is more detailed in Bravo (2004). As many strips as output component are necessary, output is defined as yk=g(xk,v k), where xk∈Rnis the system state and vk∈Rpvis the measurement noise. A measurement yk∈Rpis assumed, a consistent states set with the measurement is given and denoted as Xyk: Xyk={x∈Rn:yk∈g(x, V )} Xyk(i) is defined as consistent states set with the i-th component of yk: Xyk(i)={x∈Rn:yk(i)∈gi(x, V )} where gi(x, V ) denote the i-th component of g(x, V ). Taking into account the preceding definitions the following hold: Xyk⊆ p  i=1 Xyk(i) The reason for using a strip to outside bound of an interval is demonstrated in Bravo (2004). This paper only focuses in the result and the procedure to obtain a strip from a measurement. A zonotope ˆ Xkand a measurement ykare assumed. The middle of an interval or its center is defined by mid(). By using interval arithmetic a vector ci∈Rnand the scalars si,σ i∈Rmust be obtained such that: •ci=mid ∇xgi(ˆ Xk,V) •cT iˆ Xk−gi(ˆ Xk,V)⊆[si−σi,s i+σi] then, if Xe yk(i)=x:cT ix−yk(i)−si≤σi, the intersection between the zonotope ˆ Xkand the i-th consistent states set belongs to the intersection between the zonotope ˆ Xkand the i-th strip Xe yk(i): ˆ XkXyk(i)⊆ˆ XkXe yk(i) 4. SIMULATIONS RESULTS Once the estimate algorithms have been introduced, in this section some simulation results are presented. 4.1 Simulation results with Interval Arithmetic One advantage of the use of intervals is that the intersection between two intervals is still an interval. However, this fact does not happen with zonotopes. For the first case of study, the applied force consists in a rising step, followed by a falling step of the same magnitude in order to get the initial control signal. Fig.5 shows the simulation results when non uncertainties measurement is considered. The simulation has two different IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 70 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 71 The source of model uncertainty is assumed as the fact that the mass can change, m. Due to this fact, an interval is obtained, in which the different states can be found. The inclusion function 2of (5) is computed as follows (for simplicity the subscripts kare omitted ): f([x]) =            [x1, x1]+T[x3, x3]+T2 2[1 m,1 m]Fx [x2, x2]+T[x4, x4]+T2 2[1 m,1 m]Fy [x3, x3]+T[1 m,1 m]Fx [x4, x4]+T[1 m,1 m]Fy            (6) being Tthe sample time, α=0.1 and β=0.0001. In the calculation of ω[v], states are separated according to coordinate represented, i.e., x1and x3on one side and x2 and x4on the other one. Then it is applied to each pair separately, by using the larger of each pair. 3.2 Zonotopes The idea is to use another procedure that let a guaranteed estimation whose result can be less conservative enclosure. Once again the knowledge of a limit in the uncertainty is the initial assumption. However, this time a rectangular shape or box is not assumed. In this case a more complex shape is allowed, specifically zonotopes. In this approach, the following system obtained from (3) is used:    x1k+1 x2k+1 x3k+1 x4k+1   =   x1k+Tx 3k x2k+Tx 4k x3k x4k   +            T2 2 Fx m T2 2 Fy m TFx m TFy m            (7) where x1is the position x,x2is the position y,x3is the velocity ˙x, and x4is the velocity ˙y. The rest of the notation is the same as the one used in (3). According to the algorithm shown in Bravo (2004), the following steps are considered: •Step 1: Use an inclusion function in order to narrow the trajectory of the uncertainty nonlinear system. •Step 2: Calculate a bound of the consistent state set. •Step 3: Calculate a tight bound of the intersection set. The K¨uhn method K¨uhn (1998) is used for step 1. This method provides a zonotope that encloses the possible system state. A strip for variable set of feasible states is estimated from the measurements obtained in the second step. And finally, in the step 3, intersection between 2According to other authors can be represented as [f] zonotope and strip is calculated. Thereby, how to calculate a strip for a possible states is shown below. The way to obtain a state set through a measurement is more detailed in Bravo (2004). As many strips as output component are necessary, output is defined as yk=g(xk,v k), where xk∈Rnis the system state and vk∈Rpvis the measurement noise. A measurement yk∈Rpis assumed, a consistent states set with the measurement is given and denoted as Xyk: Xyk={x∈Rn:yk∈g(x, V )} Xyk(i) is defined as consistent states set with the i-th component of yk: Xyk(i)={x∈Rn:yk(i)∈gi(x, V )} where gi(x, V ) denote the i-th component of g(x, V ). Taking into account the preceding definitions the following hold: Xyk⊆ p  i=1 Xyk(i) The reason for using a strip to outside bound of an interval is demonstrated in Bravo (2004). This paper only focuses in the result and the procedure to obtain a strip from a measurement. A zonotope ˆ Xkand a measurement ykare assumed. The middle of an interval or its center is defined by mid(). By using interval arithmetic a vector ci∈Rnand the scalars si,σ i∈Rmust be obtained such that: •ci=mid ∇xgi(ˆ Xk,V) •cT iˆ Xk−gi(ˆ Xk,V)⊆[si−σi,s i+σi] then, if Xe yk(i)=x:cT ix−yk(i)−si≤σi, the intersection between the zonotope ˆ Xkand the i-th consistent states set belongs to the intersection between the zonotope ˆ Xkand the i-th strip Xe yk(i): ˆ XkXyk(i)⊆ˆ XkXe yk(i) 4. SIMULATIONS RESULTS Once the estimate algorithms have been introduced, in this section some simulation results are presented. 4.1 Simulation results with Interval Arithmetic One advantage of the use of intervals is that the intersection between two intervals is still an interval. However, this fact does not happen with zonotopes. For the first case of study, the applied force consists in a rising step, followed by a falling step of the same magnitude in order to get the initial control signal. Fig.5 shows the simulation results when non uncertainties measurement is considered. The simulation has two different IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 70 sampling times, a synchronization time of 1 second (the GPS refreshing measurement time) and a sampling time of 20 milliseconds, which corresponds to the observer period. Fig. 5. Coordinate X without uncertainties in GPS. For the second case of study, GPS measurement is considered unreliable. In this case, the use of a differential GPS with a accuracy of ±0.2mis assumed. Simulation results can be observed in Fig.6. As well as in the first case, the synchronization time is also 1 second and the observer sampling time is 20 milliseconds. Fig. 6. Coordinate X with uncertainties in GPS. 4.2 Simulation results with Zonotopes As in the preceding subsection, only results on axis xare displayed. The applied force also consists in a rising step, followed by a falling step of the same magnitude in order to get the initial control signal, reaching then a steady state. Simulations have been carried out considering two different conditions too: an accurate GPS signal and a GPS accuracy of 0.2 meters. The results are shown in Fig.7 and 8 respectively. The behavior of the zonotope could be appreciated in Fig.7. At the beginning the region defined by the zonotope grows until the reliable measurements are received. In those moments the zonotope size is reduced. In the final case of study, an inaccurate GPS measurement is considered. The results are shown in Fig. 8, where a non drastic reduction is observed. 5. CONCLUSIONS In this paper an application of two guaranteed estimators, one based on interval arithmetics and other based on zonotopes, on a quadrotor has been presented. The application of guaranteed algorithms is due to the fact that a continuous system is considered together with a discrete position measurement, whose sampling times are greater than the sampling time used in the system model. Thus, if only the sensor information is used, there will be time periods where the system will be in openloop, which may be dangerous taking into account that the system is unstable. In order to have a position estimation during these periods of time, an uncertainty model has been considered and the feasible states that system can reach have been contemplated. In this case the estimation with interval arithmetic looks a better performance than the one obtained with the zonotope approach. It could be due to the model used, the uncertainties considered or even how the algorithms deal with those uncertainties. As future work will aim to improve the algorithm with zonotopes in order to get a better performance than the one presented in this work. REFERENCES Abbott, E. and Powell, D. (1999). Land-vehicle navigation using gps. Proceedings of the IEEE, 87(1). Bravo, J. (2004). Control predictivo no lineal robusto basado en t´ecnicas intervalares. Ph.D. thesis. Brunke, S. and Campbell, M. (2002). Estimation architecture for future autonomous vehicles. Proceedings of the American Control Conference. Combastel, C. (2003). A state bounding observer based on zonotopes. Proceedings of European Control Conference. Farrelly, J. and Wellstead, P. (1996). Estimation of vehicle lateral velocity. Proceedings of the 1996 IEEE International Conference on Applications. Jaulin, L. (2002). Nonlinear bounded-error state estimation of continuous-time systems. Automatica, 38, 1079– 1082. K¨uhn, W. (1998). Rigorously computed orbits of dynamical systems without the wrapping effect. J. Brit. Med. Assoc. Ra¨ısi, T., Efimov, D., and Zolghadri, A. (2012). Interval state estimation for a class of nonlinear systems. IEEE Transactions on Automatic Control, 57(1). Scholte, E. and Campbell, M. (2002). On-line nonlinear guaranteed estimation with application to a high performance aircraft. Proceedings of the American Control Conference. IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 71 72 Ramón A. García et al. / IFAC-PapersOnLine 48-9 (2015) 067–072 Fig. 7. Zonotope X ˙ Xwithout uncertainties in GPS. Fig. 8. Zonotope X ˙ Xwith uncertainties in GPS IFAC ACNAAV 2015 June 10-12, 2015. Seville, Spain 72