A New Position Control Strategy for VTOL UAVs using IMU and GPS measurements

# A New Position Control Strategy for VTOL UAVs using IMU and GPS measurements

[    [
###### Abstract

We propose a new position control strategy for VTOL-UAVs using IMU and GPS measurements. Since there is no sensor that measures the attitude, our approach does not rely on the knowledge (or reconstruction) of the system orientation as usually done in the existing literature. Instead, IMU and GPS measurements are directly incorporated in the control law. An important feature of the proposed strategy, is that the accelerometer is used to measure the apparent acceleration of the vehicle, as opposed to only measuring the gravity vector, which would otherwise lead to unexpected performance when the vehicle is accelerating (i.e. not in a hover configuration). Simulation results are provided to demonstrate the performance of the proposed position control strategy in the presence of noise and disturbances.

a]Andrew Roberts, a,b]Abdelhamid Tayebi

Department of Electrical and Computer Engineering, University of Western Ontario, London, Ontario, Canada, N6A 3K7.

Department of Electrical Engineering, Lakehead University, Thunder Bay, Ontario, Canada P7B 5E1

Key words:  VTOL, UAV, vector measurments, position control,

11footnotetext: This paper was not presented at any IFAC meeting. Corresponding author A. Tayebi. Tel. +1 (807) 343-8597. Fax +1 (807) 766-7243.

## 1 Introduction

The design of position controllers for Vertical take-off and landing (VTOL) unmanned airborne vehicles (UAVs) has been the focus of several research groups, which has resulted in significant breakthroughs in this field, for example see [Abdessameud and Tayebi, 2010], [Aguiar and Hespanha, 2007], [Frazzoli et al., 2000], [Hauser et al., 1992], [Hua et al., 2009], [Pflimlin et al., 2007] and [Roberts and Tayebi, 2011a]. Existing position controllers, usually require that the system states are accurately known or measured, namely the position, linear velocity, angular velocity and the orientation. For outdoor applications a global positioning system (GPS) mounted to the system can be used to provide the position and velocity measurements, while the angular velocity is obtained using a gyroscope which is included in the inertial measurement unit (IMU) in addition to an accelerometer and a magnetometer. However, there does not exist any sensor that provides directly the orientation of a rigid body. Motivated by this problem, the study of rigid-body attitude estimation has seen substantial breakthroughs due to the efforts of the research community (see, for instance, [Mahony et al., 2008] ). However, we are not aware of any work in the literature, providing a rigorous results for the combination of an attitude observer and a position controller for VTOL-UAVs.
To address this shortcoming, there has been some effort to design position control algorithms which do not directly require the measurement of the system attitude. For example, in [Roberts and Tayebi, 2011b] the authors propose a position control law which utilizes a number of vector measurements as a means to eliminate the requirement for the attitude measurement. By vector measurements we are referring to the body-referenced measurements of vectors whose coordinates are known in the inertial frame. Since the vector measurements contain information about the system orientation, it has been shown that they can be applied directly to the position controller thereby eliminating the need of the observer completely. Consequently, the resulting vector-measurement-based position control laws do not require the direct measurement of the system attitude, nor do they require an attitude observer which provides practitioners with a simpler, reduced order closed loop system, with accompanying proofs for stability.
Unfortunately, the vector-measurements based position control strategy can be susceptible to a problem associated with the lack of sensors which can provide suitable vector measurements. This shortcoming stems from the fact that the two sensors most commonly used to provide vector measurements are the magnetometer and accelerometer, used to provide body-referenced measurements of the Earths magnetic field and gravity vector, respectively. However, in order to satisfy the requirement that the accelerometer provides a measurement of the gravity vector only, one must assume that the body-fixed frame is non-accelerating. It is clear that this condition is not guaranteed to be satisfied in some applications involving VTOL UAVs.
Of course, this practical limitation is relevant to both vector-measurement-based position controllers and attitude observers which use accelerometers. Fortunately, this limitation has led to the development of a new class of attitude observers which uses the accelerometer (and magnetometer) to provide vector-measurements. This type of observer acknowledges the fact that the accelerometer measures a combination of the gravity vector and the acceleration of the rigid-body in the body-fixed frame. This combination of the gravity vector and linear acceleration in the inertial reference frame is commonly referred to as the apparent acceleration. This inertial vector violates the requirement of many of the vector-measurement-based attitude observers, since the system acceleration is not known in the inertial frame of reference. In order to deal with the fact that the inertial vector is unknown, this type of attitude observer uses the velocity of the rigid-body (assumed to be measurable using, for instance, a GPS) in addition to the signals obtained from an IMU. These attitude observers, which are often referred to as velocity-aided attitude observers, can be found in [Bonnabel et al., 2008], [Martin and Salaün, 2008] and [Martin and Salaün, 2010] with local stability proofs , and in [Hua, 2010] with almost semiglobal stability results.
In this paper we propose a new position control approach which obviates the requirement of the system attitude measurement by using the vector measurements directly in the control law. We specifically use a magnetometer and accelerometer to provide the two vector measurements. The accelerometer is used to measure the system apparent acceleration, rather than the gravity vector only. Using our proposed approach, we show that, upon a suitable choice of the control gains, all system states remain bounded, and the system position converges to a constant reference position. Our proposed control strategy 1) does not require direct measurement of the system attitude; 2) does not require the use of an attitude observer; 3) uses an accelerometer to provide a vector measurement without limiting the motion of the system to a near-hover state.

## 2 Background

In this section we present some of the necessary mathematical details we use throughout the paper. In section 2.1 we describe two commonly used attitude representations (rotation matrices and unit-quaternion). In section 2.2 we define functions which are necessary in developing the proposed control laws.

### 2.1 Attitude Representation

To represent the orientation of the aircraft (rigid-body), we define two reference frames: An inertial frame , which is rigidly attached the Earth (assumed flat), and a body frame which is rigidly attached to the aircraft center of gravity (COG). The orthonormal basis of is taken such that the axis is directed towards the front of the aircraft (or rigid body), the axis is taken towards the starboard (right) side, and the axis is directed downwards (opposite the direction of the system thrust).
Throughout the paper we often refer to the orientation of the rigid-body, by which we mean the relative angular position of with respect to . The goal of the attitude representation is to mathematically describe the orientation of the rigid-body. The unit-quaternion, which is a unit vector on , is given by , where is the quaternion-scalar and is the quaternion-vector, and is the set of unit-quaternion defined by

 Q≡{Q∈R×R3, ∥Q∥=1}. (1)

Let , denote two unit-quaternion; then the quaternion product of and , denoted by is defined by the following operation

 Q3=Q1⊙Q2=(η1η2−qT1q2,η1q2+η2q1+S(q1)q2). (2)

The set of unit-quaternion forms a group with the quaternion multiplication operation , with the quaternion inverse , and identity element .  The unit-quaternion is an over-parameterization of the the special group of orthogonal matrices of dimension three , defined as

 (3)

that is, the transformation from the quaternion space to , given by the following Rodrigues formula:

 R(Q)=I+2S(q)2−2ηS(q). (4)

where is a skew-symmetric matrix given in the next section, is a two-to-one map, i.e., .

### 2.2 Skew symmetric matrices and bounded functions

Let . We define the skew-symmetric matrix such that , where denotes the vector cross product. Several useful properties of this skew-symmetric matrix are given below:

 S(x)2 =xxT−xTxI3×3, (5) S(Rx) =RS(x)RT,R∈SO(3), (6) S(x)y =−S(y)x = x×y, (7) λ(S(x)2) =[0,−∥x∥2,−∥x∥2], (8)

where denotes the eigenvalues of the matrix .
Consider the bounded, differentiable function, denoted as , which satisfies the following properties:

 uTh(u)>0∀u∈R3,∥u∥∈(0,∞),0≤∥h(u)∥<10<∥ϕ(u)∥≤1}∀u∈R3,∥u∥∈[0,∞), (9)

where . Throughout the paper we make use of one particular example given by . Using this definition for one can derive the expression .

## 3 Position Control Using GPS and IMU Measurements

Using the above mathematical background, we will now proceed to formulate the problem and define the position control laws. A number of steps are taken which are grouped into various sections. In section 3.1 we define the system model. In section 3.2 we formulate the problem and state some necessary assumptions. Section 3.3 provides an attitude extraction algorithm which allows us to specify a desired system attitude based upon the values of the position and velocity error, and section 3.4 defines attitude error functions. Finally, in section 3.5 we describe the position control laws.

### 3.1 Equations of Motion

To model the system translational dynamics, we let denote the position and velocity, respectively, of the vehicle COG expressed in the inertial frame . For this problem we assume that the body-referenced angular velocity vector is available as a control input. We consider the following VTOL UAV model:

 ˙p =v, (10) ˙v =μ+δ,μ=ge3−utRTe3, (11) ˙Q =12[−qTηI3×3+S(q)]ω, (12)

where , is the system thrust, is the system mass, , is the gravitational acceleration, and is a disturbance which is dependent on aerodynamic drag forces. The control input of the system is defined as . The system output is defined as where is the signal obtained using an accelerometer, is a signal obtained using a magnetometer, and is the magnetic field of the surrounding environment (assumed constant). Note that the system attitude (or ) is not assumed to be a known output of the system.

We consider a well-known model for the accelerometer model (which includes forces due to linear acceleration ) which is given by

 b2=R(˙v−ge3)=−ute3+Rδ=Rr2, (13)

where is the inertial referenced system apparent acceleration, which satisfies

 ˙v = ge3+r2,r2 = −utRTe3+δ. (14)

In the development of attitude observers, it is often assumed that the system is near hover (or ) in order to assume that the accelerometer measures the direction of the gravity vector. Also, in most situations the aerodynamic disturbance vector is not included in the model. However, for the VTOL UAV model, one can easily see that if the aerodynamic disturbance is neglected, or we assume that , then the accelerometer signal provides the measurement , which is the constant vector multiplied by the system thrust. In this case the use of the accelerometer seems trivial since its measurement is known a priori and does not contain any information about the system attitude. Therefore, we see that for the VTOL UAV system, the assumption that the accelerometer measures only the gravity vector may be a dangerous assumption which may lead to unexpected performance, even in the case where . In fact, it seems that the utility of the accelerometer measurements is related to the measurement of the vector since the accelerometer measures . For this reason we believe that it is important to include a model of the aerodynamic disturbances.

### 3.2 Problem Formulation

Let denote a desired reference position, which is assumed to be constant (or slowly-varying), and let . Our main objective is to develop a control law for the system inputs and , using the available system outputs , such that the system states and are bounded and . For the position control design we first require the following assumptions are satisfied.

###### Assumption 1

There exist positive constants and such that and

###### Assumption 2

Given two positive constants, and , there exists a positive constant such that where

The second assumption is satisfied if is non-vanishing and is not collinear to the magnetic field vector . In the case where , the system velocity dynamics become (which corresponds to the rigid body being in a free-fall state) which is not likely in normal circumstances. When this assumption is satisfied, it follows that is positive definite. Furthermore, if this assumption is satisfied, the value of can be arbitrarily increased by increasing the values of and .
In addition to this assumption, we also require some conditions on the aerodynamic force vector .

###### Assumption 3 (Aerodynamic Forces)

In light of the fact that the disturbance force is due to aerodynamic forces exerted on the vehicle we make the following simplifying assumptions:

1. The aerodynamic disturbance is dissipative with respect to the system translational kinetic energy and satisfies .

2. The aerodynamic disturbance force is only dependent on the system translational velocity, and there exist a positive constant such that

3. There exists positive constants and such that .

Assumption 3(3) and 3(3) can be realized when the system is operating in an environment where the exogenous airflow is negligible (no wind). Assumption 3(3) can be satisfied when the system geometry is sufficiently symmetrical such that the system aerodynamic forces do not significantly depend on the system orientation. Although this assumption may be reasonable for certain VTOL type aircraft, for example the ducted-fan, this assumption may not be the case with certain systems, for example fixed wing aircraft, where the system aerodynamics depend largely on the orientation of the vehicle. Now that we have established the required assumptions, let us consider the model for the system acceleration from (3.1). Due to the underactuated nature of this system, the translational acceleration is driven by the system thrust and orientation . That is, if was a control input, setting would satisfy the objectives (since ). However, since is a function of the system state, we define as the desired acceleration, and introduce the new error signal

 ~μ=μ−μd. (15)

Subsequently, a new objective is to force in order to obtain the desired translational dynamics. Since the signal is dependent on the system thrust and attitude, based upon the value of the desired acceleration we wish to obtain a suitable desired attitude, denoted as , and system thrust , such that the following equation is satisfied

 μd=ge3−utRTde3, (16)

where is the rotation matrix corresponding to the unit-quaternion , as defined by (2.1). An extraction method satisfying these requirements, which has been previously given in [Roberts and Tayebi, 2011a], is described in the following section.

### 3.3 Desired Attitude and Thrust Extraction

In this section, given a value of the desired acceleration , we seek to obtain the value of the desired orientation (or equivalently in terms of the unit-quaternion ) such that equation (3.2) is satisfied. To solve this problem we use an attitude and thrust extraction algorithm which has been previously proposed by [Roberts and Tayebi, 2011a]: Given where ,

 L={μd∈R3;μd=col[0,0,μd3];μd3∈[g,∞)}, (17)

then, one solution for the thrust and attitude where , which satisfies (3.2) is given by

 ut =∥μd−ge3∥, (18) ηd =(12(1+g−eT3μd∥μd−ge3∥))1/2, (19) qd =12∥μd−ge3∥ηdS(μd)e3. (20)

The extracted attitude has the time-derivative

 ˙Qd=12[−qTdηdI3×3+S(qd)]ωd, (21)

where the desired angular velocity is given by

 ωd = M(μd)˙μd, (22) M(μd) = 14η2du4t(−4S(μd)e3eT3+4η2dutS(e3)+2S(μd) (23) −2eT3μdS(e3))S(μd−ge3)2.

### 3.4 Attitude Error

To represent the relative orientation of the desired attitude with respect to the actual attitude , we let and denote the unknown attitude error which is defined by

 ~Q=Q⊙Q−1d,~R=R(~Q)=RTdR, (24)

where is the unit quaternion obtained using (3.3) and (3.3). In light of and , as defined by (3.1) and (3.3), respectively, the time derivative of the attitude error is found to be

 ˙~Q=12[−~qT~ηI+S(~q)]~ω,˙~R=−S(~ω)~R, (25)
 ~ω=RTd(ω−ωd), (26)

where is the desired angular velocity as defined by (3.3). One of the objectives of the control design is to force the system orientation to the desired attitude, or in terms of the rotation matrices, to force (and therefore ), in order to obtain the desired translational dynamics. As mentioned in section 2.1, this corresponds to two possible solutions for the unit-quaternion which are given by . The multiplicity of equilibrium solutions is manageable since our objectives are satisfied for both values of the unit-quaternion.

### 3.5 Position Controller

The position controller design is based upon a value of the desired system translational acceleration, which is specified by the virtual control law . Using the calculated values for the position error , and the system velocity , the value of the desired acceleration is obtained. This desired acceleration is directly related to a corresponding desired rigid-body orientation and thrust, denoted by and , respectively, which is obtained using the attitude and thrust extraction method described in section 3.3. The desired attitude given in the parametrization, denoted as , is subsequently obtained using with (2.1). Since the system attitude is not known, we incorporate the use of a special filter which is driven by the value of the linear velocity . We let denote the filter state variable which corresponds to the system velocity , and define the error function .
Although the system linear velocity is known, the use of the signal through the error function , for an appropriate choice of the estimation law can can be viewed as a function of the system acceleration in terms of the unknown signal . Since this vector is known in the body fixed frame (measured using an accelerometer, ), the filter variable through the error function can be used with the accelerometer to provide information related to the system attitude. After these steps, the remaining control design is focused on forcing the actual system attitude to the desired attitude using the control input . The proposed control law is given as follows:

 ω =M(μd)(fμd−kvϕ(v)RTd(b2+ute3))+ψ, (27) fμd =−kpϕ(ep)v+kvϕ(v)(kph(ep)+kvh(v)), (28) ψ =γ1S(Rdr1)b1+γ2k1S(Rd(v−^v))b2, (29) ˙^v =ge3+RTdb2+k1(v−^v)+1k1RTdS(b2)ψ, (30) μd =−kph(ep)−kvh(v), (31)

where , is the function defined by (23), is the bounded function defined in section 2.2, , and is obtained from the value of using the attitude extraction algorithm defined in section 3.3.

###### Theorem 1

Consider the system given by (3.1)-(3.1), where we apply the control laws and as defined by (27), where and are chosen such that . Let assumptions 2 and 3 be satisfied. Then the system thrust is bounded and non-vanishing such that

 0

and for all initial conditions (or equivalently ), there exists positive constants such that for , , , the system states and are bounded and .

Proof: We begin by first proving the upper and lower bounds on the thrust control input . Since the function is bounded by unity, the norm of the virtual control law is bounded by . Since the thrust control input is given by , and and are chosen such that , one easily arrives at the lower and upper bounds for described in the theorem. A nice consequence of the boundedness of , is that the function defined by (23), which is used in the expression for the desired angular velocity , is also bounded. In fact, in [Roberts and Tayebi, 2011a] the authors show that the norm of this matrix satisfies

 ∥M(μd)∥ ≤ √2/c–t. (33)

We now focus our attention on the dynamics of the position error and the system velocity . Let , where is the function defined by (3.1). In light of the choice for , the derivatives of the position error and velocity can be written as

 ˙ep = v,˙v = −kph(ep)−kvh(v)+~μ+δ. (34)

As previously mentioned, the velocity observer error is considered as a function of the apparent acceleration vector . In fact, we define an error function associated to the apparent acceleration which is given by

 ~r2=k1~v−(I−~R)r2. (35)

Another important error function which we will focus on is the attitude error function , or equivalently , which defines the relative orientation between the actual system attitude and the desired attitude. To prove the theorem, we will construct a Lyapunov function in terms of the error functions , , and , in order to show that all of these states tend to zero. Since the dynamics of (or equivalently ), and are somewhat complicated, we will begin by first simplifying the expressions for their derivatives. In order to analyze the dynamics of the attitude error, it is sufficient to study the derivative of the quaternion-scalar . This is also desired since the derivative of the quaternion scalar can be less complicated than the derivative of the quaternion vector. As a starting point, the derivative of can be found from (25) to be where and . To find a result for the desired angular velocity we first use the results (34), in addition to the derivative of the bounded function , denoted as as defined in section (2.2), to differentiate the virtual control law to obtain Simplifying this result, we obtain the following expression for the desired angular velocity

 ωd=M(μd)(fμd−kvϕ(v)δ−kvϕ(v)~μ). (36)

Recall the control input uses the function , given by (29). Using (35), the property (2.2) and the fact , can be rewritten as

 ψ=Rd(γ1S(r1)~Rr1+γ2S(r2)~Rr2+γ2S(~r2)~Rr2). (37)

Finally, using the expression for the control input , the error function , in addition to (36), (37) and the fact , we find the derivative To further simplify this result, we first recognize that in light of the definition of the rotation matrix from (2.1) and the property , one can find . Therefore, using the expression for the matrix defined by assumption 2, we obtain

 ˙~η =~η~qTW~q−γ22~qTS(~r2)~Rr2 −kv2~qTRTdM(μd)ϕ(v)((I−~R)δ+~μ). (38)

Note that due to assumption 2, the matrix is positive-definite. We now shift our focus to study the dynamics of the error function . In light of the expression for from (3.1), the expression for from (30), the attitude error dynamics from (25)-(26), the expressions (27), (36), and using the fact that , we obtain

 ˙~r2 =−k1~r2−(I−~R)˙r2 +kvRTdS(b2)M(μd)ϕ(v)((I−~R)δ+~μ). (39)

A commonality between the dynamic equations for and , is that they both depend on the error functions and . These two error functions can both be expressed in terms of the attitude error using the quaternion vector part , which will be a useful characteristic later in the Lyapunov analysis. To describe this relationship we define two functions, such that

 ~μ = f1(ut,~η,~q)~q,(I−~R)x = f2(x,~η,~q)~q, (40)

where . Using the definition of , in addition to the expressions for and from (3.1) and (3.2), respectively, one can find and . Based upon these definitions and the fact that , we find the following upper bounds for these two functions

 ∥f1(ut,~η,~q)∥≤2¯ct, ∥f2(x,~η,~q)∥≤2∥x∥, (41)

We now propose the following Lyapunov function candidate:

 V=γkp(√1+eTpep−1)+γ2vTv+γkr2~rT2~r2+γq(1−~η2), (42)

where and are positive constants. In light of (34), (38), (39), we have

 ˙V =−γkvvTh(v)+γvTδ−γkrk1~rT2~r2−2γq~η2~qTW~q +γkvkr~rT2RTdS(b2)M(μd)ϕ(v)(f1(ut,~η,~q) +f2(δ,~η,~q))~q−γkr~rT2f2(˙r2,~η,~q)~q+γvTf1(ut,~η,~q)~q +γqkv~η~qTRTdM(μd)ϕ(v)(f1(ut,~η,~q)+f2(δ,~η,~q))~q +γ2γq~η~qTS(~r2)~Rr2. (43)

Now, we wish to show that for an appropriate choice of the control gains, is guaranteed to be non-positive. However, this objective is a bit involved, and therefore requires we study the bound of several functions used in the expression of . We begin this analysis by defining the function scalar function .
Based upon the definition of from (42), the states and are bounded by as follows . Therefore, in light of assumption 3(3), one can conclude that

 ∥δ(v)∥≤c1σ(t)2/γ. (44)

Due to the bounds of the functions and from (41), and the definition of from (3.1) we also find

 ∥f1(ut,~η,~q)+f2(δ,~η,~q)∥ ≤2(γ¯ct+c1σ(t)2)/γ (45) ∥b2∥ ≤(γ¯ct+c1σ(t)2)/γ. (46)

Given these bounds, we now apply Young’s inequality to a number of the undesired terms in the expression for :

 (47)
 (48)

where the norm of is given by (33). To determine the bound of the term involving the time-derivative of , we first derive the expression for to be

 ˙r2=−˙utRTe3+utRTS(e3)ω+˙δ=−1ut(μd−ge3)T(−kpϕ(ev)v−kvϕ(v)f1(ut,~η,~q)~q+kvϕ(v)(kph(ep)+kvh(v))−kvϕ(v)δ)RTe3+utRTS(e3)(M(μd)(−kpϕ(ev)v−kvϕ(v)~Rδ+kvϕ(v)(kph(ep)+kvh(v)))+γ1S(Rdr1)b1+γ2RdS(~r2)~Rr2+γ2RdS(r2)~Rr2)+˙δ. (49)

Due to the bounds of the functions , , the (upper and lower) bounds of the thrust control input , the bound of from assumption 3(3), the bound of from (46) (same as the bound of ), and the bound of from (44), we find that there exists five positive constants , such that the norm of is bounded by . However, for the sake of simplicity, from this result we further conclude that there exists positive constants and such that As a result of this analysis, we again use Young’s inequality to establish the following bounds:

 γkr~rT2f2(˙r2,~η,~q)~q≤γkrϵ32~rT2~r2+2γkrϵ3(c3+c4σ(t)4)2~qT~q, (50)
 γ2γq~η~qTS(~r2)~Rr2≤γ22γqϵ42~rT2~r2+γq2γ2ϵ4(¯ctγ+c1σ(t)2)2~η2~qT~q, (51)
 γqkv~η~qTRTdM(μd)ϕ(v)(f1(ut,~η,~q)+f2(δ,~η,~q))~q≤γqkv(√2c–t)(2(γ¯ct+c1σ(t)2)γ)|~η|~qT~q≤2√2γqkv(γ¯ct+c1σ(t)2)γc–t|~η|~qT~q, (52)

Recall from assumption 2 that the norm of the matrix has a lower bound which is denoted as . Therefore, in light of the lower bounds defined above, and assumption 3(3) we find the expression is bounded by

 ˙V(t)≤−γvTh(v)(kv−ϵ1/2)−γkr~rT2~r2(k1−ϵ2kv+ϵ32−ϵ4γ22γq2γkr)−γq~η2~qT~q(2cw−1~η2(α1(t)ϵ1+α2(t)ϵ2+α3(t)ϵ3)−α4(t)ϵ4−2√γ¯c2t(γ+σ(t)2)1/2~η2−2√2kv(γ¯ct+c1σ(t)2)γc–t|~η|), (53)
 α1(t) =2√γ¯c2t√γ+σ(t)2/γ