arXiv is now an independent nonprofit! Learn more
License: arXiv.org perpetual non-exclusive license
arXiv:1809.00037v1 [cs.RO] 31 Aug 2018

Estimation for Quadrotors

Stefanie Tellex Brown University Andy Brown Udacity, Inc. Sergei Lupashin Fotokite

This document describes standard approaches for filtering and estimation for quadrotors, created for the Udacity Flying Cars course. We assume previous knowledge of probability and some knowledge of linear algebra. We do not assume previous knowledge of Kalman filters or Bayes filters. This document derives an EKF for various models of drones in 1D, 2D, and 3D. We use the EKF and notation as defined in Thrun et al. 2005. We also give pseudocode for the Bayes filter, the EKF, and the Unscented Kalman filter [14]. The motivation behind this document is the lack of a step-by-step EKF tutorial that provides the derivations for a quadrotor helicopter. The goal of estimation is to infer the drone’s state (pose, velocity, accelleration, and biases) from its sensor values and control inputs. This problem is challenging because sensors are noisy. Additionally, because of weight and cost issues, many drones have limited on-board computation so we want to estimate these values as quickly as possible. The standard method for performing this method is the Extended Kalman filter, a nonlinear extension of the Kalman filter which linearizes a nonlinear transition and measurement model around the current state. However the Unscented Kalman filter is better in almost every respect: simpler to implement, more accurate to estimate, and comparable runtimes.

1 Averaging

When performing estimation, the first thing one might think of is averaging the sensor measurements. We consider averaging for the 1D case, where the range sensor value, ztz_{t} is directly observing the state (the height of the drone), xtx_{t}, both of which are scalars. The notation xt^\hat{x_{t}} means our estimate of the true state xtx_{t}. Then the average is:

xt^\displaystyle\hat{x_{t}} =1tΣ0tzt\displaystyle=\frac{1}{t}\Sigma_{0}^{t}z_{t} (1)

However this form requires saving all sensor measurements from the beginning of time, requiring memory and computation that grows linearly in the number of observations. Instead we rewrite the update recursively in terms of our previous estimate, x^t1\hat{x}_{t-1}, allowing us to transform x^t1\hat{x}_{t-1} and ztz_{t} to x^\hat{x}:

xt^\displaystyle\hat{x_{t}} =1t[x^t1×(t1)+zt].\displaystyle=\frac{1}{t}\left[\hat{x}_{t-1}\times(t-1)+z_{t}\right]. (2)

The recursive form, where the next estimate is written in terms of the previous estimate, allows us to perform a constant time update, and requires us to store a constant amount of information about past sensor values (only the previous estimate, x^t1\hat{x}_{t-1} - in other words, our estimate is Markov). However this form updates slowly when the drone’s state changes. Imagine the drone is stationary for one minute, and then starts moving. Then 1t\frac{1}{t} will be very small, and the estimate will take a long time to move to the new value (longer and longer, the longer the drone is running). Instead, we want an average that moves more quickly. One way to achieve a faster updating average is to simply average the old estimate with the new sensor value:

xt^=x^t1+zt2\displaystyle\hat{x_{t}}=\frac{\hat{x}_{t-1}+z_{t}}{2} (3)

A more general formulation is a weighted average of the old value and the new observation:

xt^=(α)x^t1+(1α)zt\displaystyle\hat{x_{t}}=(\alpha)\hat{x}_{t-1}+(1-\alpha)z_{t} (4)

This formulation allows the designer to tune the parameter α\alpha, which weights how much to weight the old estimate compared to the new sensor value.

2 Bayes Filter

None of the above formulations update the estimate based on the drone’s movements. Intuitively, if we have told the drone to go up, for instance, then our belief about the drone’s position should also go up. The Bayes Filter gives us a way to incorporate motion prediction into our state estimate. First we predict the next state, given a control input, the current state, and a model of how the system evolves over time. We do not maintain a point estimate but rather a belief or distribution of the estimate. Then, we revise our prediction with an update from an observation. The update method takes the previous estimate from prediction, and an observed sensor value. It returns a new distribution that takes into account the sensor value, using a model of how the sensor works. Below you will see pseudocode for the prediction and update steps for the filter, following Thrun et al. 2005. In typical use one would call predict after determining the control input, and update after reading a sensor value. The BayesFilter function is illustrative only; in real life, one might call predict in the sensor callback for example.

Algorithm 1 General Bayes Filter algorithm. For specific filters such as the Kalman Filter or the particle filter, the representation for belbel and bel¯\bar{bel} changes, and the corresponding mathematical updates take specific computational forms.
1: function Predict(bel(xt1),ut,Δtbel(x_{t-1}),u_{t},\Delta t)
2:   bel¯(xt)=p(xt|ut,xt1)bel(xt1)dxt1\bar{bel}(x_{t})=\int p(x_{t}|u_{t},x_{t-1})bel(x_{t-1})dx_{t-1}
3:   return bel¯(xt)\bar{bel}(x_{t})
4: function Update(bel¯(xt),zt\bar{bel}(x_{t}),z_{t})
5:   bel(xt)=ηp(zt|xt)bel(xt)bel(x_{t})=\eta p(z_{t}|x_{t})bel(x_{t})
6:   return bel(xt)bel(x_{t})
7: function BayesFilter
8:   ut=ComputeControl(bel(xt1))u_{t}=\textsc{ComputeControl}(bel(x_{t-1}))
9:   bel¯(xt)=Predict(bel(xt1),ut,Δt)\bar{bel}(x_{t})=\textsc{Predict}(bel(x_{t-1}),u_{t},\Delta t)
10:   zt=ReadSensor()z_{t}=\textsc{ReadSensor}()
11:   bel(xt)=Update(bel¯(xt),zt)bel(x_{t})=\textsc{Update}(\bar{bel}(x_{t}),z_{t})

3 E(KF)

The Kalman Filter and Extended Kalman Filter make the assumption that the distributions over belief state are Gaussian, represented as a mean and covariance matrix. Compare to the Bayes’ filter, the distributions belbel and bel¯\bar{bel} are represented as a mean and covariance matrix. The KF assumes that the transition and observation models are linear, and can be defined by a matrix. The EKF is the extension to the nonlinear case, where we take use the Jacobian matrix of the transition and observation functions to compute a point-wise linear estimate, and then do the same updates as the Kalman Filter. We define the Extended Kalman Filter (EKF) algorithm following Thrun et al. 2005. We refactor it to include separate Predict and Update methods and to use our notation. We also unify the KF and EKF algorithm pseudocode. The transition or prediction covariance is QtQ_{t}; the measurement covariance is RtR_{t}. These matrices are often taken to be constant, but also sometimes people change them over time depending on the sensor model.

The EKF and KF are closely related. For the KF, the matrix GtG_{t} is constant every iteration of the function and does not need to be recomputed each time (except for Δt\Delta t). Another way to say it is the implementation of gg^{\prime} ignores its state input. Similarly, for the KF, the matrix HtH_{t} is constant every iteration of the function and does not need to be recomputed. Another way to say it is the function hh^{\prime} ignores its state input. For the EKF, these matrices change each iteration, because it linearizes around the current state.

Algorithm 2 E(KF) algorithm.
1: function Predict(μt1,Σt1,ut,Δt\mu_{t-1},\Sigma_{t-1},u_{t},\Delta t)
2:   μ¯t=g(ut,μt1)\bar{\mu}_{t}=g(u_{t},\mu_{t-1})
3:   Gt=g(ut,xt,Δt)G_{t}=g^{\prime}(u_{t},x_{t},\Delta t)
4:   Σ¯t=GtΣt1GtT+Qt\bar{\Sigma}_{t}=G_{t}\Sigma_{t-1}G_{t}^{T}+Q_{t}
5:   return μ¯t,Σ¯t\bar{\mu}_{t},\bar{\Sigma}_{t}
6: function Update(μ¯t,Σ¯t,zt\bar{\mu}_{t},\bar{\Sigma}_{t},z_{t})
7:   Ht=h(μ¯t)H_{t}=h^{\prime}(\bar{\mu}_{t})
8:   Kt=Σ¯tHtT(HtΣ¯tHtT+Rt)1K_{t}=\bar{\Sigma}_{t}H_{t}^{T}(H_{t}\bar{\Sigma}_{t}H_{t}^{T}+R_{t})^{-1}
9:   μt=μ¯t+Kt(zth(μ¯t))\mu_{t}=\bar{\mu}_{t}+K_{t}(z_{t}-h(\bar{\mu}_{t}))
10:   Σt=(IKtHt)Σ¯t\Sigma_{t}=(I-K_{t}H_{t})\bar{\Sigma}_{t}
11:   return μt,Σt\mu_{t},\Sigma_{t}
12: function ExtendedKalmanFilter
13:   ut=ComputeControl(μt1,Σt1)u_{t}=\textsc{ComputeControl}(\mu_{t-1},\Sigma_{t-1})
14:   μ¯t,Σ¯t=Predict(μt1,Σt1,ut,Δt)\bar{\mu}_{t},\bar{\Sigma}_{t}=\textsc{Predict}(\mu_{t-1},\Sigma_{t-1},u_{t},\Delta t)
15:   zt=ReadSensor()z_{t}=\textsc{ReadSensor}()
16:   μt,Σt=Update(μ¯t,Σ¯t,zt)\mu_{t},\Sigma_{t}=\textsc{Update}(\bar{\mu}_{t},\bar{\Sigma}_{t},z_{t})

4 Unscented Kalman Filter

The Unscented Kalman Filter [14] is similar to the EKF in that it handles nonlinear transition models gg and measurement models hh. However there are no Jacobians! Instead of linearizing around the current estimate, the Unscented Kalman Filter picks magic “sigma points” which are sample points chosen according to the current state estimate and covariance. Then these sigma points are passed through the nonlinear transition or observation function, and the sample mean and sample covariance of the sigma points is used to construct a new Gaussian μ\mu and σ\sigma. The pseudocode here follows Kandepu et al. 2008 but removes the augmentation for tracking moving prediction and measurement covariance, and uses our notation.

Sigma points are computed using the matrix SS which is defined from the covariance matrix, Σt\Sigma_{t}. SiS_{i} denotes the ithith colum of the matrix SS.

S=Σt\displaystyle S=\sqrt{\Sigma_{t}} (5)

We define the sigma points Xi,tXtX_{i,t}\in X_{t} as follows:

Xi,t={=μt,i=0=μt+γSi,i=1,,N=μtγSiN,i=N+1,,2N\displaystyle X_{i,t}=\left\{\begin{array}[]{lll}&=\mu_{t},&i=0\\ &=\mu_{t}+\gamma S_{i},&i=1,\dots,N\\ &=\mu_{t}-\gamma S_{i-N},&i=N+1,\dots,2N\end{array}\right.

Note that they are defined using the mean and covariance matrix of the distributions; by picking several representative points SiS_{i} away from the mean, we can use a relatively small set of points to represent the entire distribution. The weights when computing the sample mean, wimw_{i}^{m} are:

wim={=λN+λi=0=12(N+λ),i=1,,2N\displaystyle w_{i}^{m}=\left\{\begin{array}[]{lll}&=\frac{\lambda}{N+\lambda}&i=0\\ &=\frac{1}{2(N+\lambda)},&i=1,\dots,2N\\ \end{array}\right.

The weights when computing the sample covariance, wicw_{i}^{c} are:

wic={=λN+λ+(1α2+β)i=0=12(N+λ),i=1,,2N\displaystyle w_{i}^{c}=\left\{\begin{array}[]{lll}&=\frac{\lambda}{N+\lambda}+(1-\alpha^{2}+\beta)&i=0\\ &=\frac{1}{2(N+\lambda)},&i=1,\dots,2N\\ \end{array}\right.

The parameters are defined as:

γ=N+λ\displaystyle\gamma=\sqrt{N+\lambda} (13)
λ=α2(N+κ)N\displaystyle\lambda=\alpha^{2}(N+\kappa)-N (14)

See Kandepu et al. 2008 for tuning suggestions when implementing the filter.

The pseudocode for the Unscented Kalman Filter is given in Algorithm 3. Compared to the Bayes’ filter, we use the computed sigma points to represent the distribution bel¯\bar{bel}. However we use a mean and covariance to represent the distribution belbel at the beginning and end of each iteration of the filter.

Algorithm 3 Unscented Kalman Filter.
1: function ComputeSigmas(μt,Σt\mu_{t},\Sigma_{t})
2:   return X0,t,,X2N,tX_{0,t},\dots,X_{2N,t} following Equation 4.
3: function Predict(μt1,Σt1,ut,Δt\mu_{t-1},\Sigma_{t-1},u_{t},\Delta t)
4:   Xt1=ComputeSigmas(μt1,Σt1)X_{t-1}=\textsc{ComputeSigmas}(\mu_{t-1},\Sigma_{t-1})
5:   i=02NX¯i,t=g(Xi,t1,ut,Δt)\forall_{i=0}^{2N}\bar{X}_{i,t}=g(X_{i,t-1},u_{t},\Delta t)
6:   return X¯t\bar{X}_{t}
7: function Update(X¯t,zt\bar{X}_{t},z_{t})
8:   μ¯t=i=02N(wimX¯i,t)\bar{\mu}_{t}=\sum_{i=0}^{2N}(w_{i}^{m}\bar{X}_{i,t})
9:   Σ¯t=i=02Nwic(Xi,tμ¯t)(Xi,tμ¯t)T+Qt\bar{\Sigma}_{t}=\sum_{i=0}^{2N}w_{i}^{c}(X_{i,t}-\bar{\mu}_{t})(X_{i,t}-\bar{\mu}_{t})^{T}+Q_{t}
10:   i=12NZi,t=h(X¯i,t)\forall_{i=1}^{2N}Z_{i,t}=h(\bar{X}_{i,t})
11:   μz=Σi=02NwimZi,t\mu^{z}=\Sigma_{i=0}^{2N}w_{i}^{m}Z_{i,t}
12:   Σtz=Σi=02Nwic(Zi,tμz)(Zi,tμz)T+Rt\Sigma_{t}^{z}=\Sigma_{i=0}^{2N}w_{i}^{c}(Z_{i,t}-\mu^{z})(Z_{i,t}-\mu^{z})^{T}+R_{t}
13:   Σtxz=Σ02Nwic(X¯i,tμ¯t)(Zi,tμz)T\Sigma_{t}^{xz}=\Sigma_{0}^{2N}w_{i}^{c}(\bar{X}_{i,t}-\bar{\mu}_{t})(Z_{i,t}-\mu^{z})^{T}
14:   Kt=Σtxz(Σtz)1K_{t}=\Sigma_{t}^{xz}(\Sigma_{t}^{z})^{-1}
15:   μt=μt¯+Kt(ztμz)\mu_{t}=\bar{\mu_{t}}+K_{t}(z_{t}-\mu^{z})
16:   Σt=Σ¯tKtΣtzKtT\Sigma_{t}=\bar{\Sigma}_{t}-K_{t}\Sigma_{t}^{z}K_{t}^{T}
17:   return μt,Σt\mu_{t},\Sigma_{t}
18: function UnscentedKalmanFilter
19:   ut=ComputeControl(μt1,Σt1)u_{t}=\textsc{ComputeControl}(\mu_{t-1},\Sigma_{t-1})
20:   X¯t=Predict(μt1,Σt1,ut,Δt)\bar{X}_{t}=\textsc{Predict}(\mu_{t-1},\Sigma_{t-1},u_{t},\Delta t)
21:   zt=ReadSensor()z_{t}=\textsc{ReadSensor}()
22:   μt,Σt=Update(X¯t,zt)\mu_{t},\Sigma_{t}=\textsc{Update}(\bar{X}_{t},z_{t})

5 One Dimensional Quad

To implement the filters on specific vehicles, we need to define the state transition function, gg and the measurement model, hh. We will do this three times for increasingly more realistic models of a quadrotor. First, we define a 1D quad model, where the quadroter is moving in zz, but not in xx or yy. It has one control input, the downward pointing thrust, and a noisy range sensor. The intention is that this is identical to the 1D quad used in the controls lesson. The state is then the position, zz and velocity, z˙\dot{z}:

xt=[z˙z]\displaystyle x_{t}=\left[\begin{array}[]{c}\dot{z}\\ z\end{array}\right]

We define the control input as directly setting the accelleration, z¨\ddot{z}:

ut=[z¨]\displaystyle u_{t}=\left[\begin{array}[]{c}\ddot{z}\end{array}\right]

5.1 Transition Model

Then we define a transition function g(xt,ut,Δt)g(x_{t},u_{t},\Delta t) which returns a new xt+1x_{t+1}:

g(xt,ut,Δt)\displaystyle g(x_{t},u_{t},\Delta t) =[xt,z˙+ut,z¨×Δtxt,z+xt,z˙×Δt]\displaystyle=\left[\begin{array}[]{c}x_{t,\dot{z}}+u_{t,\ddot{z}}\times\Delta t\\ x_{t,z}+x_{t,\dot{z}}\times\Delta t\end{array}\right]
=[10Δt1][z˙z]+[Δt0][z¨]\displaystyle=\left[\begin{array}[]{cc}1&0\\ \Delta t&1\end{array}\right]\left[\begin{array}[]{c}\dot{z}\\ z\end{array}\right]+\left[\begin{array}[]{cc}\Delta t\\ 0\end{array}\right]\left[\begin{array}[]{c}\ddot{z}\end{array}\right]
We can rewrite it in terms of the AtA_{t} and BtB_{t} matrix, as in a conventional Kalman filter. In this form we see that the control update is linear because it can be written in this form.
=Atxt+Btut\displaystyle=A_{t}x_{t}+B_{t}u_{t} (27)

Then xt2x_{t}\in\mathbb{R}^{2} and g(xt,ut,Δt)2g(x_{t},u_{t},\Delta t)\in\mathbb{R}^{2}. So g(xt,ut)g^{\prime}(x_{t},u_{t}) is a 2×22\times 2 matrix, defined as the partial derivative of gg with respect to xtx_{t} for each component in xtx_{t}.

g(xt,ut,Δt)\displaystyle g^{\prime}(x_{t},u_{t},\Delta t) =[xt,z˙gz˙(xt,ut,Δt)xt,zgz˙(xt,ut,Δt)xt,z˙gz(xt,ut,Δt)xt,zgz(xt,ut,Δt)]\displaystyle=\left[\begin{array}[]{cccc}\frac{\partial}{\partial x_{t,\dot{z}}}g_{\dot{z}}(x_{t},u_{t},\Delta t)&\frac{\partial}{\partial x_{t,z}}g_{\dot{z}}(x_{t},u_{t},\Delta t)\\ \frac{\partial}{\partial x_{t,\dot{z}}}g_{z}(x_{t},u_{t},\Delta t)&\frac{\partial}{\partial x_{t,z}}g_{z}(x_{t},u_{t},\Delta t)\end{array}\right]
=[10Δt1]\displaystyle=\left[\begin{array}[]{cccc}1&0\\ \Delta t&1\end{array}\right]

Since this function is linear, the Jacobian is a constant matrix except for Δt\Delta t, and just reduces to the AtA_{t} matrix.

5.2 Measurement Model

Next we assume the drone has a range sensor pointed downwards at the ground, at z=0z=0. Then ztz_{t} is the range value, rr:

zt=[r]\displaystyle z_{t}=\left[\begin{array}[]{c}r\end{array}\right]

Then we define a measurement function h(xt)h(x_{t}) which returns a new ztz_{t}:

h(xt)\displaystyle h(x_{t}) =[xt,z]\displaystyle=\left[\begin{array}[]{c}x_{t,z}\end{array}\right]
=[01][z˙z]\displaystyle=\left[\begin{array}[]{cc}0&1\end{array}\right]\left[\begin{array}[]{c}\dot{z}\\ z\end{array}\right]
=Ct[z˙z]\displaystyle=C_{t}\left[\begin{array}[]{c}\dot{z}\\ z\end{array}\right]
=Ctxt\displaystyle=C_{t}x_{t} (40)

Finally we define h(xt)h^{\prime}(x_{t}), the Jacobian of hh with respect to xtx_{t}. The Jacobian is a 1×21\times 2 matrix.

h(xt)\displaystyle h^{\prime}(x_{t}) =[xt,z˙hr(xt)xt,zhr(xt)]\displaystyle=\left[\begin{array}[]{cccc}\frac{\partial}{\partial x_{t,\dot{z}}}h_{r}(x_{t})&\frac{\partial}{\partial x_{t,z}}h_{r}(x_{t})\end{array}\right]
=[01]\displaystyle=\left[\begin{array}[]{cccc}0&1\end{array}\right]

Note that in this case the Jacobian does not depend at all on the input xtx_{t}. This is because this system is linear, so the EKF linearization will boil back down to a regular Kalman filter.

6 Two Dimensional Quad

We define a 2D quad model where the quadrotor is operating at a fixed height zz. It has a range sensor pointing sideways, and it can move by rotating about its center, the angle ϕ\phi. Additionally it can move right and left in yy. Then we define the state xtx_{t} as the state transition function as follows:

xt=[ϕy˙y]\displaystyle x_{t}=\left[\begin{array}[]{c}\phi\\ \dot{y}\\ y\end{array}\right]

We define the control input as directly setting the angle, ϕ\phi.

ut=[ϕ]\displaystyle u_{t}=\left[\begin{array}[]{c}\phi\end{array}\right]

6.1 Transition Model

Then we define g(xt,ut,Δt)g(x_{t},u_{t},\Delta t) and returns a new xt+1x_{t+1}:

g(xt,ut,Δt)=[ut,ϕxt,y˙sin(xt,ϕ)×Δtxt,y+xt,y˙×Δt]\displaystyle g(x_{t},u_{t},\Delta t)=\left[\begin{array}[]{c}u_{t,\phi}\\ x_{t,\dot{y}}-\sin(x_{t,\phi})\times\Delta t\\ x_{t,y}+x_{t,\dot{y}}\times\Delta t\end{array}\right]

Then xt3x_{t}\in\mathbb{R}^{3} and g(xt,ut,Δt)3g(x_{t},u_{t},\Delta t)\in\mathbb{R}^{3}. So g(xt,ut)g^{\prime}(x_{t},u_{t}) is a 3×33\times 3 matrix.

Finally, we take the partial derivative of gg with respect to xtx_{t} for each component in xtx_{t}.

g(xt,ut,Δt)\displaystyle g^{\prime}(x_{t},u_{t},\Delta t) =[xt,ϕgϕxt,y˙gϕxt,ygϕxt,ϕgy˙xt,y˙gy˙xt,ygy˙xt,ϕgyxt,y˙gyxt,ygy]\displaystyle=\left[\begin{array}[]{cccc}\frac{\partial}{\partial x_{t,\phi}}g_{\phi}&\frac{\partial}{\partial x_{t,\dot{y}}}g_{\phi}&\frac{\partial}{\partial x_{t,y}}g_{\phi}\\ \frac{\partial}{\partial x_{t,\phi}}g_{\dot{y}}&\frac{\partial}{\partial x_{t,\dot{y}}}g_{\dot{y}}&\frac{\partial}{\partial x_{t,y}}g_{\dot{y}}\\ \frac{\partial}{\partial x_{t,\phi}}g_{y}&\frac{\partial}{\partial x_{t,\dot{y}}}g_{y}&\frac{\partial}{\partial x_{t,y}}g_{y}\end{array}\right]
=[xt,ϕut,ϕxt,y˙ut,ϕxt,yut,ϕxt,ϕxt,y˙sin(xt,ϕ)×Δtxt,y˙xt,y˙sin(xt,ϕ)×Δtxt,yxt,y˙sin(xt,ϕ)×Δtxt,ϕxt,y+xt,y˙×Δtxt,y˙xt,y+xt,y˙×Δtxt,yxt,y+xt,y˙×Δt]\displaystyle=\left[\begin{array}[]{cccc}\frac{\partial}{\partial x_{t,\phi}}u_{t,\phi}&\frac{\partial}{\partial x_{t,\dot{y}}}u_{t,\phi}&\frac{\partial}{\partial x_{t,y}}u_{t,\phi}\\ \frac{\partial}{\partial x_{t,\phi}}x_{t,\dot{y}}-\sin(x_{t,\phi})\times\Delta t&\frac{\partial}{\partial x_{t,\dot{y}}}x_{t,\dot{y}}-\sin(x_{t,\phi})\times\Delta t&\frac{\partial}{\partial x_{t,y}}x_{t,\dot{y}}-\sin(x_{t,\phi})\times\Delta t\\ \frac{\partial}{\partial x_{t,\phi}}x_{t,y}+x_{t,\dot{y}}\times\Delta t&\frac{\partial}{\partial x_{t,\dot{y}}}x_{t,y}+x_{t,\dot{y}}\times\Delta t&\frac{\partial}{\partial x_{t,y}}x_{t,y}+x_{t,\dot{y}}\times\Delta t\end{array}\right]
=[000cos(xt,ϕ)Δt100Δt1]\displaystyle=\left[\begin{array}[]{cccc}0&0&0\\ -\cos(x_{t,\phi})\Delta t&1&0\\ 0&\Delta t&1\end{array}\right]

6.2 Measurement Model

Next we assume the drone has a range sensor pointed sideways at a wally=5\text{wall}_{y}=5. Then ztz_{t} is the range value, rr:

zt=[r]\displaystyle z_{t}=\left[\begin{array}[]{c}r\end{array}\right]

Then we define h(xt)h(x_{t}) and returns a new ztz_{t}:

h(xt)=[wallyxt,ycos(xt,ϕ)]\displaystyle h(x_{t})=\left[\begin{array}[]{c}\frac{\text{wall}_{y}-x_{t,y}}{\cos(x_{t,\phi})}\end{array}\right]

Finally we define h(xt)h^{\prime}(x_{t}), the Jacobian of hh with respect to xtx_{t}.

h(xt)\displaystyle h^{\prime}(x_{t}) =[xt,ϕhr(xt)xt,y˙hr(xt)xt,yhr(xt)]\displaystyle=\left[\begin{array}[]{cccc}\frac{\partial}{\partial x_{t,\phi}}h_{r}(x_{t})&\frac{\partial}{\partial x_{t,\dot{y}}}h_{r}(x_{t})&\frac{\partial}{\partial x_{t,y}}h_{r}(x_{t})\end{array}\right]
=[[wallyxt,y]sin(xt,ϕ)cos2(xt,ϕ)01cos(xt,ϕ)]\displaystyle=\left[\begin{array}[]{cccc}\frac{\left[\text{wall}_{y}-x_{t,y}\right]\sin(x_{t,\phi})}{\cos^{2}(x_{t,\phi})}&0&\frac{-1}{\cos(x_{t,\phi})}\end{array}\right]

Using the above math, we have implemented an EKF in Python for this model, and showed it is able to estimate position and velocity using the drone’s simulated noisy range sensor.

7 Three Dimensional Quad

Variable Description
utu_{t} The control input at time tt.
xtx_{t} The state at time tt.
ztz_{t} The observation at time tt.
xt,variablex_{t,\mbox{variable}} The scalar value of the state vector at the index corresponding to variable. Similar notation for utu_{t} and ztz_{t}.
Δt\Delta t The elapsed time between updates in seconds.
QtQ_{t} Transition model covariance
RtR_{t} Measurement model covariance
ϕ\phi Roll
θ\theta Pitch
ψ\psi Yaw
RbgR_{bg} Rotation matrix from body to global
RgyR_{gy} Rotation matrix from global to yaw frame
xyx^{y} the xx coordinate in the yaw frame.
xx The xx coordinate in the global frame.
xbx^{b} The xx coordinate in the body frame.
g(xt,ut,Δt)g(x_{t},u_{t},\Delta t) The transition function.
h(xt)h(x_{t}) The observation function.
Table 1: Table of variables.

The state will be position and velocity which we will estimate with the GPS. We will use the magnetometer to estimate yaw. We will represent position and velocity in the global frame, and yaw that goes from body to local in the state. We would also then use a multirate Kalman Filter [4, 12] to perform updating, with separate updates for the GPS, the magnetometer, and the IMU (accelerometer plus rate gyro).

We track yaw as the heading from magnetic north. So it is the reading one would get if one reads from a compass.

xt=[xyzx˙y˙z˙ψ]\displaystyle x_{t}=\left[\begin{array}[]{c}x\\ y\\ z\\ \dot{x}\\ \dot{y}\\ \dot{z}\\ \psi\end{array}\right]

Then utu_{t} is the acceleration in the body frame, where ψ˙\dot{\psi} is global frame yaw.

ut\displaystyle u_{t} =[x¨by¨bz¨bψ˙]\displaystyle=\left[\begin{array}[]{c}\ddot{x}^{b}\\ \ddot{y}^{b}\\ \ddot{z}^{b}\\ \dot{\psi}\end{array}\right]

7.1 Attitude Filter

Markley 2003 gives an EKF in quaternions for attitude estimation. Higgins 1975 compares complementary filters to Kalman filters, in a way not specific to quads or attitude estimation. Quan 2017 says to use either a complementary filter or Kalman filter for attitude estimation and gives very terse, hard to understand math. Johansen and Kristiansen 2017 describes the MEKF model used in a quadrotor with adaptive fading. Crassidis et al. 2007 gives a survey of nonlinear attitude estimation including the MEKF. Nowicki et al. 2015 compares the complementary filter with an EKF for attitude estimation on mobile phones. They find that the complementary filter is simpler to implement, but the EKF is able to achieve in most cases better accuracy. Both have comparable processor loads.

We assume the state we are tracking is the vehicle’s attitude, that is roll ϕ\phi and pitch, θ\theta. Then the observation, ztz_{t} consists of the gyro angular velocity and pitch and roll angles as estimated from the accelerometer.

zt=[θϕpq]\displaystyle z_{t}=\left[\begin{array}[]{c}\theta\\ \phi\\ p\\ q\end{array}\right]

The accelerometer estimates for θ\theta and ϕ\phi are in the global frame, but the velocities from the gyro are in the body frame. Our approach follows the Linear Complementary Filter from Quan 2017. We assume that θ\theta and ϕ\phi are small, so that the turn rates measured by the gyro in the body frame approximate the global turn rates. In other words,

[ϕ˙θ˙ψ˙][pqr]\displaystyle\left[\begin{array}[]{c}\dot{\phi}\\ \dot{\theta}\\ \dot{\psi}\end{array}\right]\approx\left[\begin{array}[]{c}p\\ q\\ r\end{array}\right]

The state is the roll angle and pitch angle:

xt=[θϕ]\displaystyle x_{t}=\left[\begin{array}[]{c}\theta\\ \phi\end{array}\right]

7.1.1 Linear Complementary Filter

We define a linear complementary filter following Quan 2017. Here τ\tau is a time constant and TsT_{s} is the filter sampling period:

θ^t=ττ+Ts(θ^t1+Tszt,θ˙)+Tsτ+Tszt,θ\displaystyle\hat{\theta}_{t}=\frac{\tau}{\tau+T_{s}}\left(\hat{\theta}_{t-1}+T_{s}z_{t,\dot{\theta}}\right)+\frac{T_{s}}{\tau+T_{s}}z_{t,\theta} (86)

Similarly for roll:

ϕ^t=ττ+Ts(ϕ^t1+Tszt,ϕ˙)+Tsτ+Tszt,ϕ\displaystyle\hat{\phi}_{t}=\frac{\tau}{\tau+T_{s}}\left(\hat{\phi}_{t-1}+T_{s}z_{t,\dot{\phi}}\right)+\frac{T_{s}}{\tau+T_{s}}z_{t,\phi} (87)

We do not estimate the yaw with a complementary filter because we will use the magnetometer and do it in the GPS.

The above math assumes that the angular velocity (which is in body frame) can be used directly as angular rotation.

7.1.2 Nonlinear Complementary Filter

For the nonlinear complementary filter, following Quan 2017 Section 9.1.3, we use the state to define a quaternion, qtq_{t}, for the euler angles for ϕ\phi, θ\theta and ψ\psi. Then we can define dqdq to be the quaternion that consists of the measurement of the angular rates from the IMU in the body frame, following Equation 84 in Diebel 2006. Using these two, we can define a predicted quaternion, q¯t\bar{q}_{t} as follows:

q¯t=dqqt\displaystyle\bar{q}_{t}=dq*q_{t} (88)

Finally we can define θ¯t\bar{\theta}_{t} and ϕ¯t\bar{\phi}_{t} as follows:

θ¯t=Pitch(q¯t)\displaystyle\bar{\theta}_{t}=Pitch(\bar{q}_{t}) (89)
ϕ¯t=Roll(q¯t)\displaystyle\bar{\phi}_{t}=Roll(\bar{q}_{t}) (90)

Using these predicated estimates, we can compute the non-linear complementary filter as above.

θ^t=ττ+Ts(θ¯t1+Tszt,θ˙)+Tsτ+Tszt,θ\displaystyle\hat{\theta}_{t}=\frac{\tau}{\tau+T_{s}}\left(\bar{\theta}_{t-1}+T_{s}z_{t,\dot{\theta}}\right)+\frac{T_{s}}{\tau+T_{s}}z_{t,\theta} (91)

Similarly for roll:

ϕ^t=ττ+Ts(ϕ¯t1+Tszt,ϕ˙)+Tsτ+Tszt,ϕ\displaystyle\hat{\phi}_{t}=\frac{\tau}{\tau+T_{s}}\left(\bar{\phi}_{t-1}+T_{s}z_{t,\dot{\phi}}\right)+\frac{T_{s}}{\tau+T_{s}}z_{t,\phi} (92)

7.2 Transition Model

We define the transition function in terms of the rotation matrix RbgR_{bg} which rotates from the body frame to the global frame. As described in Diebel 2006, there are 12 different orders one could perform the rotation; we follow the convention from aerospace of using the 1,2,31,2,3 order for roll, pitch, and yaw.

This matrix is defined as follows, taken from the transpose (or inverse) of Diebel 2006, equation 67.

Rbg=[cosθcosψsinϕsinθcosψcosϕsinψcosϕsinθcosψ+sinϕsinψcosθsinψsinϕsinθsinψ+cosϕcosψcosϕsinθsinψsinϕcosψsinθcosθsinϕcosθcosϕ]\displaystyle R_{bg}=\left[\begin{array}[]{ccc}\cos\theta\cos\psi&\sin\phi\sin\theta\cos\psi-\cos\phi\sin\psi&\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi\\ \cos\theta\sin\psi&\sin\phi\sin\theta\sin\psi+\cos\phi\cos\psi&\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi\\ -\sin\theta&\cos\theta\sin\phi&\cos\theta\cos\phi\end{array}\right]

Then the transition function is:

g(xt,ut,Δt)\displaystyle g(x_{t},u_{t},\Delta t) =[xt,x+xt,x˙Δtxt,y+xt,y˙Δtxt,z+xt,z˙Δtxt,x˙xt,y˙xt,z˙gΔtxt,ψ]+[000000000000Rbg[0:]0Rbg[1:]0Rbg[2:]00001]utΔt\displaystyle=\left[\begin{array}[]{c}x_{t,x}+x_{t,\dot{x}}\Delta t\\ x_{t,y}+x_{t,\dot{y}}\Delta t\\ x_{t,z}+x_{t,\dot{z}}\Delta t\\ x_{t,\dot{x}}\\ x_{t,\dot{y}}\\ x_{t,\dot{z}}-g\Delta t\\ x_{t,\psi}\\ \end{array}\right]+\left[\begin{array}[]{cccc}0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ R_{bg}[0:]&&&0\\ R_{bg}[1:]&&&0\\ R_{bg}[2:]&&&0\\ 0&0&0&1\end{array}\right]u_{t}\Delta t

Then we take the Jacobian:

g(xt,ut,Δt)\displaystyle g^{\prime}(x_{t},u_{t},\Delta t) =[100Δt0000100Δt0000100Δt0000100xt,ψ(xt,x˙+Rbg[0:]ut[0:3]Δt)000010xt,ψ(xt,y˙+Rbg[1:]ut[0:3]Δt)000001xt,ψ(xt,z˙+Rbg[2:]ut[0:3]Δt)0000001]\displaystyle=\left[\begin{array}[]{ccccccc}1&0&0&\Delta t&0&0&0\\ 0&1&0&0&\Delta t&0&0\\ 0&0&1&0&0&\Delta t&0\\ 0&0&0&1&0&0&\frac{\partial}{\partial x_{t,\psi}}\left(x_{t,\dot{x}}+R_{bg}[0:]u_{t}[0:3]\Delta t\right)\\ 0&0&0&0&1&0&\frac{\partial}{\partial x_{t,\psi}}\left(x_{t,\dot{y}}+R_{bg}[1:]u_{t}[0:3]\Delta t\right)\\ 0&0&0&0&0&1&\frac{\partial}{\partial x_{t,\psi}}\left(x_{t,\dot{z}}+R_{bg}[2:]u_{t}[0:3]\Delta t\right)\\ 0&0&0&0&0&0&1\\ \end{array}\right]
=[100Δt0000100Δt0000100Δt0000100Rbg[0:]ut[0:3]Δt000010Rbg[1:]ut[0:3]Δt000001Rbg[2:]ut[0:3]Δt0000001]\displaystyle=\left[\begin{array}[]{ccccccc}1&0&0&\Delta t&0&0&0\\ 0&1&0&0&\Delta t&0&0\\ 0&0&1&0&0&\Delta t&0\\ 0&0&0&1&0&0&R^{\prime}_{bg}[0:]u_{t}[0:3]\Delta t\\ 0&0&0&0&1&0&R^{\prime}_{bg}[1:]u_{t}[0:3]\Delta t\\ 0&0&0&0&0&1&R^{\prime}_{bg}[2:]u_{t}[0:3]\Delta t\\ 0&0&0&0&0&0&1\end{array}\right]

We define RbgR^{\prime}_{bg} as xt,ψ\frac{\partial}{\partial x_{t,\psi}}, defined as Equation 71 from Diebel 2006:

Rbg=[cosθsinψsinϕsinθsinψcosϕcosψcosϕsinθsinψ+sinϕcosψcosθcosψsinϕsinθcosψcosϕsinψcosϕsinθcosψ+sinϕsinψ000]\displaystyle R^{\prime}_{bg}=\left[\begin{array}[]{ccc}-\cos\theta\sin\psi&-\sin\phi\sin\theta\sin\psi-\cos\phi\cos\psi&-cos\phi\sin\theta\sin\psi+\sin\phi\cos\psi\\ \cos\theta\cos\psi&\sin\phi\sin\theta\cos\psi-\cos\phi\sin\psi&\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi\\ 0&0&0\end{array}\right]

7.3 Measurement Model

We provide measurement models for the GPS and Magnetometer. We use the IMU as a control input so do not provide a measurement model for it here.

7.3.1 GPS

We assume we get position and velocity from the GPS. We considered using heading from the GPS, but this does not take into account the drone’s orientation, only the direction of travel. Hence we are removing it from the observation.

zt\displaystyle z_{t} =[xyzx˙y˙z˙]\displaystyle=\left[\begin{array}[]{c}x\\ y\\ z\\ \dot{x}\\ \dot{y}\\ \dot{z}\\ \end{array}\right]

Then the measurement model is:

h(xt)=[xt,xxt,yxt,zxt,x˙xt,y˙xt,z˙]\displaystyle h(x_{t})=\left[\begin{array}[]{c}x_{t,x}\\ x_{t,y}\\ x_{t,z}\\ x_{t,\dot{x}}\\ x_{t,\dot{y}}\\ x_{t,\dot{z}}\\ \end{array}\right]

Then the partial derivative is the identity matrix, augmented with a vector of zeros for xt,ϕh(xt)\frac{\partial}{\partial x_{t,\phi}}h(x_{t}):

h(xt)=[100000001000000010000000100000001000000010]\displaystyle h^{\prime}(x_{t})=\left[\begin{array}[]{ccccccc}1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ 0&0&1&0&0&0&0\\ 0&0&0&1&0&0&0\\ 0&0&0&0&1&0&0\\ 0&0&0&0&0&1&0\\ \end{array}\right]

7.3.2 Magnetometer

We assume we get a reading from the magnetometer reporting yaw in the global frame. (This measurement may need to be computed using roll and pitch from the attitude filter and the mag vector.)

zt=[ψ]\displaystyle z_{t}=\left[\begin{array}[]{c}\psi\end{array}\right]
h(xt)=[xt,ψ]\displaystyle h(x_{t})=\left[\begin{array}[]{c}x_{t,\psi}\end{array}\right]

Again since this is linear, the derivative is a matrix of zeros and ones.

h(xt)=[0000001]\displaystyle h^{\prime}(x_{t})=\left[\begin{array}[]{ccccccc}0&0&0&0&0&0&1\end{array}\right]

7.4 Further Information

Based on not tracking acceleration as part of the state, we will use the accelerometer and gyro inputs as control inputs. Note that Erdem and Ercan 2015 describe using the accelerometer inputs in the measurement or control input phases for an EKF over camera and IMU, for all eight combinations. They show that it is always better to fuse both sensors in the measurement stage. We would also then use a multirate Kalman Filter [4, 12] to perform updating, with separate updates for the GPS, magnetometer, and IMU.

However this means that both acceleration and angular velocity appear in the state vector, which we decided not to do in order to simplify the math and implementation. As one example, Ardupilot does it in this way [1]. The advantage is that it does not need to connect as deeply to the control system. Additionally, one could even use different state transition models in the control compared to the EKF. One might make this decision for computational reasons, for example.

Erdem and Ercan 2015 give the math for an EKF with an IMU and camera (on a mobile device), showing the IMU and gyro measurements treated as a prediction and measurement input (all eight combinations). They show that the best perforance is obtained using the IMU as a measurement input. However the difference in pose accuracy estimation is around 1 cm1\text{\,}\mathrm{c}\mathrm{m}. Note that the Ardupilot open source code base [1] makes a different decision, using the IMU as the control/prediction update, and that Bry et al. 2012 describes this decision as a “commonly-used technique.”

We will use the North/East/Down frame where the positive xx gives the distance along the surface of the earth in the direction of north; the yy coordinate gives the distance in the direction of easteast, and zz is altitude, which is negative for distances above the surface of the earth.

We define an intermediate frame, pp, which is a quaternion initialized from the roll, θ\theta, and the pitch, ϕ\phi, from the complementary filter, and no yaw correction. We define the global frame, qq, as a quaternion filled out with the roll, θ\theta, the pitch, ϕ\phi, and the yaw, ψ\psi, from the EKF.

8 Conclusion

We have provided equations for filtering for a quadrotor helicopter, combining information and notation from Thrun et al. 2005, Quan 2017 and other sources.

References

  • ard [Accessed March 5, 2018] Extended kalman filter navigation overview and tuning, Accessed March 5, 2018. http://ardupilot.org/dev/docs/extended-kalman-filter.html#extended-kalman-filter.
  • Bry et al. [2012] Adam Bry, Abraham Bachrach, and Nicholas Roy. State estimation for aggressive flight in gps-denied environments using onboard sensing. In Robotics and Automation (ICRA), 2012 IEEE International Conference on, pages 1–8. IEEE, 2012.
  • Crassidis et al. [2007] John L Crassidis, F Landis Markley, and Yang Cheng. Survey of nonlinear attitude estimation methods. Journal of guidance, control, and dynamics, 30(1):12–28, 2007.
  • Cristi and Tummala [2000] Roberto Cristi and Murali Tummala. Multirate, multiresolution, recursive kalman filter. Signal Processing, 80(9):1945–1958, 2000.
  • Diebel [2006] James Diebel. Representing attitude: Euler angles, unit quaternions, and rotation vectors. Matrix, 58(15-16):1–35, 2006.
  • Erdem and Ercan [2015] Arif Tanju Erdem and Ali Özer Ercan. Fusing inertial sensor data in an extended kalman filter for 3d camera tracking. IEEE Transactions on Image Processing, 24(2):538–548, 2015.
  • Higgins [1975] Walter T Higgins. A comparison of complementary and kalman filtering. IEEE Transactions on Aerospace and Electronic Systems, (3):321–325, 1975.
  • Johansen and Kristiansen [2017] Tor-Aleksander Johansen and Raymond Kristiansen. Quadrotor attitude estimation using adaptive fading multiplicative ekf. In American Control Conference (ACC), 2017, pages 1227–1232. IEEE, 2017.
  • Kandepu et al. [2008] Rambabu Kandepu, Bjarne Foss, and Lars Imsland. Applying the unscented kalman filter for nonlinear state estimation. Journal of process control, 18(7-8):753–768, 2008.
  • Markley [2003] F Landis Markley. Attitude error representations for kalman filtering. Journal of guidance, control, and dynamics, 26(2):311–317, 2003.
  • Nowicki et al. [2015] Michał Nowicki, Jan Wietrzykowski, and Piotr Skrzypczyński. Simplicity or flexibility? complementary filter vs. ekf for orientation estimation on mobile devices. In Cybernetics (CYBCONF), 2015 IEEE 2nd International Conference on, pages 166–171. IEEE, 2015.
  • Quan [2017] Quan Quan. Introduction to multicopter design and control. Springer, 2017.
  • Thrun et al. [2005] Sebastian Thrun, Wolfram Burgard, and Dieter Fox. Probabilistic robotics. MIT press, 2005.
  • Wan and Van Der Merwe [2000] Eric A Wan and Rudolph Van Der Merwe. The unscented kalman filter for nonlinear estimation. In Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000. AS-SPCC. The IEEE 2000, pages 153–158. Ieee, 2000.