Design and implementation of a low-cost IoT-based agroclimatic monitoring system for greenhouses

1.   Introduction

Spacecraft, underwater vehicles, aerial vehicles, and mobile robots require accurate knowledge of their orientation with respect to a known inertial frame. Typically, attitude estimators rely on the measurements of angular velocity and known inertial vectors in the body-fixed frame. Therefore, the rigid body attitude estimation problem using angular velocity and inertial vectors measurements in a body-fixed frame has been widely studied in past research. In practice, angular velocity is measured at a higher rate than inertial vectors. In this work^*, we address the attitude estimation problem given multi-rate measurements of inertially fixed vectors and angular velocity by minimizing the "energy" stored in the state estimation errors.

^*A preliminary arxiv version of this article is available at ^[1]

One of the earliest solutions to attitude determination from vector measurements is the TRIAD algorithm, used to determine the rotation matrix from two linearly independent inertial vector measurements ^[2]. In ^[3], Wahba presented the attitude determination problem as an optimization problem using three or more vector measurements where the cost function is the sum of the squared norms of vector errors. Various methods have been proposed in the literature to solve the Wahba's problem. Davenport was the first to reduce the Wahba's problem to finding the largest eigenvalue and the corresponding eigenvector of the so-called Davenport's K-matrix ^[4]. In a similar approach, Mortari presented the EStimator of the Optimal Quaternion (ESOQ) algorithm in^[5], which provides the closed-form expressions of a $4\times4$ matrix's eigenvalues and then computes the eigenvector associated with the greatest of them, representing the optimal quaternion. The QUEST algorithm of ^[6] determines the attitude that achieves the best-weighted overlap of an arbitrary number of reference vectors. A Singular Value Decomposition (SVD) based method of solving the Wahba's problem was proposed by ^[7]. Markley also devised a Fast Optimal Matrix Algorithm (FOMA) to solve Wahba's problem in ^[8]. A coordinate-free framework of geometric mechanics was used in ^[9] to obtain a solution to Wahba's problem with robustness to measurement noise. ^[10] provides a numerical solution to Wahba's problem.

Attitude estimation methods based on minimizing "energy" stored in state estimation errors can be found in ^{[11,12,13,14]}. Prior research that has designed attitude estimation schemes based on geometric mechanics includes ^{[15,16,17,18]}. Comprehensive surveys of various attitude estimation methods are available in ^[19,20]. However, most of the aforementioned schemes for attitude estimation work only in continuous time or measurement rich environments and neglect the sparsity in the inertial vector measurements. Inertial vector measurements are usually obtained with the help of Sun (and star) sensors or magnetometers which are accurate but suffer from lower sampling rates. ^[21] provides one of the first solutions to rigid body attitude estimation with multi-rate measurements using uncertainty ellipsoids. A recursive method based on the cascade combination of an output predictor and an attitude observer can be found in ^[22]. An attitude estimation scheme on the Special Orthogonal group using intermittent body-frame vector measurements was presented in ^[23]. In ^[13], a filtering scheme in continuous-time is proposed by applying the Lagrange-d'Alembert principle on suitably formulated artificial kinetic and potential energy functions where the authors formulate filter equations assuming that inertial vector measurements and angular velocity measurements are available synchronously and continuously. A discrete-time estimation scheme in the presence of multi-rate measurements is proposed in ^[14]. This work presents the Lyapunov analysis but does not contain a variational interpretation.

The previous works in attitude estimations are presented in a continuous time domain and then discretized for numerical implementation. This voids the theoretical guarantee of asymptotic stability provided by the continuous-time estimator. Furthermore, most of the previous work in attitude estimation do not consider the realistic scenario when the inertial vectors and angular velocities are measured at a different rate. In the current work, we try to overcome these shortcoming by focusing on developing an asymptotically stable optimal geometric discrete-time attitude estimator based on the minimization of "energy" stored in the errors of the state estimators in the presence of multi-rate measurements. The measurements can be corrupted by noise and we do not assume any specific statistics (like normal distribution) on the measurement noise. However, the noise is assumed to be bounded. We represent the attitude as a rotation matrix which precludes potential singularity issues due to local coordinates (such as Euler angles). The multi-rate discrete-time filtering scheme presented here is obtained by applying the discrete-time Lagrange-D'Alembert principle ^[24] on a discrete-time lagrangian followed by a discrete-time Lyapunov analysis using a Lyapunov candidate that depends on the state estimation errors. The filtering scheme provided is asymptotically stable with an almost global region of convergence.

This paper is organized as follows. In the section 2, the attitude estimation problem is formulated as Wahba's optimization problem and then some important properties of the Wahba's cost function are presented. In the section 3, continuous-time rigid body attitude kinematics has been discretized and the propagation model for the measurements in the multi-rate measurement case is presented. Section 4 contains the application of variational mechanics to obtain a filter equation for attitude estimation. The filter equations obtained in the 4 are proven to be asymptotically stable with an almost global domain of convergence by deriving an appropriate dissipation torque using the discrete-time Lyapunov method in the 5. Filter equations are numerically verified with realistic measurements (corrupted by bounded noise) in the 6. Finally, 7 presents the concluding remarks with contributions and future work.

2.   Problem formulation and notation

2.1.   Notation and preliminaries

We define the trace inner product on $\mathbb{R}^{m\times n}$ as

The group of orthogonal frame transformations on $\mathbb{R}^3$ is defined by $\mathrm{O}(3) : = \{Q \in \mathbb{R}^{3\times 3} \; | \; \text{det}(Q) = \pm 1\}$. The Special orthogonal group on $\mathbb{R}^3$ is denoted as ${\text{SO(3)}}$ defined as ${\text{SO(3)}} : = \{R \in \mathbb{R}^{3\times 3} \; | \; R^{\text{T}} R = RR^{\text{T}} = I_{3}\}$. The corresponding Lie algebra is denoted as $\mathfrak{so}(3) : = \{ M \in \mathbb{R}^{3\times 3} \; | \; M + M^{\text{T}} = 0\}$. Let $(\cdot)^\times:\mathbb{R}^3 \rightarrow \mathfrak{so}(3) \subset \mathbb{R}^{3\times3}$ be the skew-symmetric matrix cross-product operator and denotes the vector space isomorphism between $\mathbb{R}^3$ and $\mathfrak{so}(3)$:

Further, let $\text{vex}(\cdot):\mathfrak{so}(3)\rightarrow \mathbb{R}^3$ be the inverse of $(\cdot)^\times$. $\exp{(\cdot)}:\mathfrak{so}(3)\rightarrow {\text{SO(3)}}$ is the map defined as

We define $Ad:{\text{SO(3)}}\times\mathfrak{so}(3)\rightarrow \mathfrak{so}(3)$ as

In the rest of the article, the text "Consider the time interval $[t_0, T] \subseteq \mathbb{R}^+$", indicates that the estimation process will be carried out for the time interval $[t_0, T]$ and the time interval is divided into $N$ equal sub-intervals $[t_i, t_{i+1}]$ for $i = 0, 1, \dots, N$ with $t_N = T$ with time step size is denoted as, $h : = t_{i+1} - t_i$; unless specified otherwise. Furthermore, throughout this article we use $(\cdot)^m$ to denote measured quantities. For example, if $\Omega$ is used to denote angular velocity, $\Omega^m$ denotes measured value of the angular velocity.

2.2.   Attitude determination from vector measurements

For the attitude estimation, we consider $k \in \mathbb{N}$ known and linearly independent inertial vectors in the body-fixed frame. Let's denote these vectors in the body-fixed frame by $u_j^m \in \mathbb{R}^3$ for $j = 1, \ldots, k$, where $k \geq 2$. Note, that $k \geq 2$ is necessary for determining the attitude uniquely. When $k = 2$, the cross product of the two measured vectors is used as the third independent measurement. Let $e_j\in \mathbb{R}^3$ be the corresponding known inertial vectors. We denote the true vectors in the body-fixed frame by $u_j : = R^{\text{T}} e_j$, where $R$ is the rotation matrix of the body-fixed frame with respect to the inertial frame. This rotation matrix provides a coordinate-free global and unique description of the attitude of the rigid body. Define the matrix populated by all $k$ measured vectors expressed in the body-fixed frame as column vectors as

and the corresponding inertial frame vectors as

The true body vector matrix is as below:

2.2.1.   Formulation of Wahba's cost function for instantaneous attitude determination from vector measurements

The optimal attitude determination problem using a set of vector measurements is finding an estimated rotation matrix $\hat{R} \in {\text{SO(3)}}$, such that the weighted sum of squared norms of the vector errors

is minimized. This attitude determination problem is known as Wahba's problem and consists of minimizing

with the respect to $\hat{R} \in {\text{SO(3)}}$, and the weights $w_j > 0$ for all $j \in \{1, 2, \ldots, k\}$. we can express (2.9) as

where $U^m$ is given by (2.5), $E$ is given by (2.6), and $W = \text{diag}(w_i)$ is the positive definite diagonal matrix of the weight factors for the measured directions. $W$ in (2.10) can be generalized to be any positive definite matrix. It is to be noted that, having different weights for different vector measurements is helpful as in general we can have some sensors that provide more accurate measurements than others. Therefore, it makes sense to provide more weight to those accurate measurements during attitude estimation.

2.2.2.   Properties of Wahba's cost function in the absence of measurements errors

This is an observer design that, like any other observer, can filter measurement noise from realistic measurements. As observer designs are based on perfect (deterministic) measurements for which their stability properties apply, it is useful to confirm these properties through numerical simulations. We have $U^m = U = R^{\text{T}} E$ in the absence of measurement errors or noise. Let $Q = R\hat{R}^{\text{T}} \in {\text{SO(3)}}$ denote the attitude estimation error. The following lemmas from ^[13] stated here without proof give the structure and characterization of critical points of the Wahba's cost function.

Lemma 2.1. Let rank($E$) = 3 and the singular value decomposition of $E$ be given by

and Diag$^+(n_1, n_2)$ is the vector space of $n_1\times n_2$ matrices with positive entries along the main diagonal and all the other components zero. Let $\sigma_1, \sigma_2, \sigma_3$ denote the main diagonal entries of $\Sigma_E$. Further, Let $W$ from $(2.10)$ be given by

and the first three diagonal entries of $W_0$ are given by

Then, $K = EWE^{{T}}$ is positive definite and

is its eigen-decomposition. Moreover, if $d_i \neq d_j$ for $i \neq j$ and $i, j \in \{1, 2, 3\}$ then $\langle\, I-Q, K\rangle$ is a Morse function whose set of critical points given as the solution of $S_K(Q) : = {{vex}}\left(KQ - Q^{{T}} K\right) = 0$:

and $a_i$ is the $i^{th}$ column vector of the identity matrix $I \in {\text{SO(3)}}$.

Lemma 2.2. Let $K = EWE^{{T}}$ have the properties given by lemma 2.1. Then the map ${\text{SO(3)}} \ni Q \mapsto \langle\, I-Q, K\rangle \in \mathbb{R}$ with critical points given by $(2.15)$ has a global minimum at the identity $I\in {\text{SO(3)}}$, a global maximum and two hyperbolic saddle points whose indices depend on the distinct eigenvalues $d_1, d_2,$ and $d_3$ of $K$.

3.   Discretization of attitude kinematics

3.1.   Case-1: The sampling rates are integer related

Consider the time interval $[t_0, T] \subseteq \mathbb{R}^+$. Let the true angular velocity in the body-fixed frame be denoted by $\Omega \in \mathbb{R}^3$. The true and measured angular velocities at the time instant $t_i$ will be denoted by $\Omega_i$ and $\Omega_i^m$ respectively. Further, let $U_i$ and $U_{i}^m$ denote the matrix formed by true and measured inertial vectors in the body-fixed frame at the time instant $t_i$ respectively. The assumption is that angular velocity measurements and inertial vectors measurements in the body-fixed frame are coming at a different but constant rate. In general coarse rate gyros have much higher sampling rate than that of a coarse attitude sensor. Therefore, in a realistic scenario, angular velocities are measured at a higher rate than the inertial vector measurements in the body-fixed frame. Therefore, we assume that the measurements of angular velocity ($\Omega^m$) are available after each time interval $h$ say, $\Omega^m_0, \Omega^m_1, \dots, \Omega^m_N$ while, inertial vector measurements in the body-fixed frame are available after time interval $nh, n \in N$ say, $U^m_0, U^m_n, U^m_{2n}, \dots$.

We have, $U = R^{\text{T}} E$. Therefore, at time instants $t_i$ and $t_{i+1}$, the following relations will hold true respectively; $U_i = R_i^{\text{T}} E_i, \; U_{i+1} = R_{i+1}^{\text{T}} E_{i+1}$. Here, $R_i$ and $R_{i+1}$ are the rotation matrices from body-fixed frame to inertial frame at time instants $t_i$ and $t_{i+1}$ respectively. $E_i = E_{i+1} = E$ are the corresponding known vectors expressed in the inertial frame. Note that the vectors are fixed in the inertial frame and do not change with the time. The continuous time attitude kinematics are,

We discretize the kinematics in (3.1) as follows,

where, $\exp{(\cdot)}:\mathfrak{so}(3)\mapsto {\text{SO(3)}}$ is the map defined as,

Using (2.7) and the discretization from (3.2) we get,

For the instants of time when inertial vector measurements in the body-fixed frame are not available we will use (3.4) to obtain the missing values of $U_i^m$. This implies that for the time instants $(n-1)h < t_i < nh, n \in \mathbb{N}$, by employing the propagation scheme in (3.4), we propagate direction vector measurements between the instants at which they are measured, using the angular velocity measurements that are obtained at a faster rate. We now formalize the aforementioned inertial vector measurement model as below,

Note that in the absence of measurements errors, we have $\Omega_i^m = \Omega_i, \, \forall i \in \{0, 1, \ldots, N\}$. Also, $U_i^m = U_i$ for the time instants when inertial vector measurements are available. Now, at time instant $t_0$, we have $\tilde{U}^m_{0} = U^m_0 = U_0$ and $\Omega^m_0 = \Omega_0$. Using (3.5) at time instant $t_1$, noting that $\Omega^m_1 = \Omega_1$, we get $\tilde{U}^m_{1} = \exp{\left(-\frac{h}{2}(\Omega_{0}+\Omega_1)^\times\right)}{U}_{0}$. Comparing it with (3.4), we have $\tilde{U}^m_{1} = U_1$. Using the relation from (2.7) we have $\tilde{U}^m_{1} = R_1^{\text{T}} E_1$. Similarly, combining (3.4), and (3.5), and using the relation in (2.7) we get the following relation for all $i \in \{0, 1, \ldots, N\}$ in the absence of measurement errors:

3.2.   Case-2: The sampling rates are time-varying

Consider the time interval $[t_0, T] \subseteq \mathbb{R}^+$. Let the true angular velocity in the body-fixed frame be denoted by $\Omega \in \mathbb{R}^3$. The true and measured angular velocities at the time instant $t_i$ will be denoted by $\Omega_i$ and $\Omega_i^m$ respectively. It is assumed that a measurement $U^m_0$ is available at time instant $t_0$. Now the next measurement will be available after a time interval $nh$ where $n$ is an integer variable taking values between $n_1$ and $n_2$ which are user provided. Note that since $n$ is varying, that sampling rate of inertial vector sensor is not fixed but varying over time. Considering the discrete kinematics of $R$ and $U$ from (3.2) and (3.4) we propose the following discrete-time vector measurement model,

Note that now we have obtained values of $\Omega^m_i$ and $\tilde{U}^m_i$ for each time instant $t_i$ and similar to the previous arguments, it can be proven that in the absence of measurement errors, we have $\Omega^m_i = \Omega_i$ and $\tilde{U}_i^m = R_i^{\text{T}} E_i$.

3.3.   Case-3: The sampling rates are non-integer related

Again consider the time interval $[t_0, T] \subseteq \mathbb{R}^+$. Let the angular velocity measurements $\Omega^m \in \mathbb{R}^3$ be available after each time interval of $h$. Now suppose the inertial vector measurements are available after a time interval $rh$ where $r = n_1/n_2$ is a rational number with $n_1, n_2 \in \mathbb{N}$. Let $n$ be the greatest common factor of $n_1$ and $n_2$. Let $h^\prime = nh/n_2$. Now, consider the time interval $[t_0, T] \subseteq \mathbb{R}^+$ divided into $N^\prime$ equal sub-intervals $[t^\prime_i, t^\prime_{i+1}]$ for $i = 0, 1, \dots, N^\prime$ with $t^\prime_N = T$ and let $t^\prime_{i+1} - t^\prime_i = h^\prime$. Note that a measurement $\Omega^m$ will be available after a time interval of $n_2h^\prime/n$ and a measurement $U^m$ will be available after time interval $rh = \frac{n_1}{n_2}\times\frac{n_2}{n}h^\prime = n_1h^\prime/n$.

We define the measurement model for $\Omega^m$ as

for all $i = 0, 1, \dots, N^\prime$. Now we can define the vector measurement model to be,

Note that the division of the time interval $[t_0, T]$ with the time step size of $h^\prime$ is finer than that of with the time step size of $h$ with the relation as $h = \frac{n_2}{n}h^\prime$. Therefore for each $i \in \{0, 1, \ldots, N\}$ we have $j = \frac{n_2}{n}i \in \{0, 1, \ldots, n^\prime\}$ such that $t_i = t_j^\prime$. Redefining $\Omega^m_i : = \tilde{\Omega}^m_{\frac{n_2}{n}i}$, and $\tilde{U}^m_i = \tilde{U}^m_{\frac{n_2}{n}i}$ for all $i \in \{0, 1, \ldots, N\}$. Note that now we have obtained values of $\Omega^m_i$ and $\tilde{U}^m_i$ for each time instant $t_i$ and similar to previous arguments, it can be proven that in the absence of measurement errors, we have $\Omega^m_i = \Omega_i$ and $\tilde{U}_i^m = R_i^{\text{T}} E_i$.

Observe that in any of the aforementioned three cases, the measurement models provide values $\Omega^m_i$ and $\tilde{U}^m_i$ for all the discrete-time instants $t_i, \; i \in \{0, 1, \ldots, N\}$ in the time interval $[t_0, T]$.

4.   Discrete-time optimal attitude estimator based on Lagrange-d'Alembert's principle

The value of the Wahba's cost function at each instant encapsulates the error in the attitude estimation. We can consider the Wahba's cost function as an artificial potential energy-like term. Therefore using (2.10) we have

where $\tilde{U}_i^m$ is according to the inertial vector propagation model presented in section 3. The term encapsulating the "energy" in the angular velocity estimation error is denoted by the map $\mathcal{T}^v$ defined as

where $m > 0$ is a scalar and $\hat{\Omega}_i$ stands for the value of estimated angular velocity at the discrete time instant $t_i$. We can write (4.2) in terms of the angular velocity estimation error $\omega_i : = \Omega^m_i - \hat{\Omega}_i$ as

The discrete time Lagrangian can be written as:

We use variational mechanics ^[25,26] approach to obtain an optimal estimation scheme for the constructed Lagrangian. In the absence of a dissipative term, the variational approach would result in (a generalization to the Lie group of) the Euler-Lagrange equations that we obtain in the context of optimal control. The generalization would be an Euler-Poincare equation on SO(3) ^[27]. Therefore, if the estimation process is started at time $t_0$, then the discrete-time action functional corresponding to the discrete-time Lagrangian (4.4) over the time interval $[t_0, T]$ can be expressed as

4.1.   Discrete-time attitude state estimation based on the discrete-time Lagrange-d'Alembert principle

Consider attitude state estimation in discrete-time in the presence of multirate measurements with noise and initial state estimate errors. Applying the discrete-time Lagrange-d'Alembert principle ^[24] from variational mechanics to the action functional $\mathfrak{s}_d(\mathscr{L}(\hat{R}_i, \tilde{U}_i^m, \omega_i, \omega_{i+1}))$ given by (4.5), in the presence of a dissipation term on $\omega_i = \Omega_i^m - \hat{\Omega}_i$, leads to the following attitude and angular velocity filtering scheme.

Proposition 1. Consider the time-interval $[t_0, T]$. Consider the multi-rate measurement models from Section 3 such that we have values $\Omega^m_i$ and $\tilde{U}^m_i$ for all the discrete-time instants $t_i, \; i \in \{0, 1, \ldots, N\}$ in the time interval $[t_0, T]$. Let the $W_i$ be chosen such that $K_i = E_iW_iE_i^{{T}}$ satisfies eigen-decomposition condition $(2.14)$ of lemma 2.1. Also, let $\tau_{D_i} \in \mathbb{R}^3$ denote the value of the dissipation torque at the time instant $t_i$. A discrete-time optimal filter obtained by applying the discrete-time Lagrange-d'Alembert principle would be as follows:

where $S_{L_{i}}(\hat{R}_i) = {{vex}}(L_i^T\hat{R}_i - \hat{R}_i^{{T}} L_i) \in \mathbb{R}^3$, $L_i = E_iW_i(\tilde{U}^m_i)^{{T}}$, and $(\hat{R}_0, \hat{\Omega}_0) \in {\text{SO(3)}}\times\mathbb{R}^{3\times3}$ are initial estimated states.

Proof. Consider a first variation in the discrete attitude estimate as

where $\Sigma_i \in \mathbb{R}^3$ represents a variation for the discrete attitude estimate. For fixed end-point variations, we have $\Sigma_0 = \Sigma_N = 0$. A first order approximation is to assume that $\hat{\Omega}^\times$ and $\delta\hat{\Omega}^\times$ commute. Taking the first variation of the discrete-time attitude kinematics according to the first equation of (4.6) and comparing with (4.7) we get

(4.8) can be rearranged to

(4.9) can be equivalently written as an equation in $\mathbb{R}^3$ as follows:

Note that, $\omega_i = \Omega_i^m - \hat{\Omega}_i$ gives us

Consider artificial potential energy term as expressed in (4.1). Taking its first variation with the respect to estimated attitude $\hat{R}$, we get

Consider the first variation in the artificial kinetic energy $\mathcal{T}^v(\omega_i, \omega_{i+1})$ as in (4.3) with the respect to the angular velocity estimation error,

Using the results in (4.10) and (4.11) we get

The first variation of the discrete-time action sum in (4.5) using (4.12) and (4.14) can be written as

Applying the discrete-time Lagrange-d'Alembert principle to the attitude motion, we obtain

For $0 \leq i < N$, (4.16) leads to

which in turn leads to the second filter equation.

Remark 1. Please note that, we introduced a virtual dissipation torque in the Proposition 1 and the performance of the estimator will depend on the value of dissipation torque. In the next section, we determine the value of dissipation torque via discrete-Lyapunov analysis so that the resulting estimator is almost globally asymptotically stable.

5.   Stability of the filter using the discrete Lyapunov Approach

For the Lyapunov stability of the filter equations, we need to construct a suitable Lyapunov candidate function. We use the Wahba's cost function expressed in (4.1) as the artificial potential energy which encapsulates the error in the estimation of attitude. A new term encapsulating the "energy" in the angular velocity estimation error can be constructed as the map $\mathcal{T}^l:\mathbb{R}^3\times\mathbb{R}^3\rightarrow \mathbb{R}$ defined as

where $m > 0$ is a scalar same as before. Further, (5.1) can be written in terms of angular velocity estimation error, $\omega_i : = \Omega^m_i - \hat{\Omega}_i$ as follows:

In the absence of measurement errors, we have $\tilde{U}_i^m = R_i^{\text{T}} E_i$. Therefore we can we can write (4.1) in terms of state estimation error $Q_i = R_i\hat{R}^{\text{T}}_i$ as

The weights $W_i$'s are chosen such that $K_i$ is always positive definite with distinct eigenvalues according to lemma 2.1.

Theorem 5.1. Consider the time-interval $[t_0, T]$. Consider the multi-rate measurement models from Section 3 such that we have values $\Omega^m_i$ and $\tilde{U}^m_i$ for all the discrete-time instants $t_i, \; i \in \{0, 1, \ldots, N\}$ in the time interval $[t_0, T]$. Then the estimation scheme in Proposition 1 with the following value of the dissipation torque:

leads to the estimation scheme

where $S_{L_{i}}(\hat{R}_i) = {{vex}}(L_i^{{T}}\hat{R}_i - \hat{R}_i^{{T}} L_i) \in \mathbb{R}^3$, $L_i = E_iW_i(\tilde{U}^m_i)^{{T}}$, $l > 0$, $l \neq m$ and $k_p > 0$, which is asymptotically stable at the estimation error state $(Q, \omega) : = (I, 0)$ ($Q_i = R_i\hat{R}^{\text{T}}_i$) in the absence of measurement errors. Further, the domain of attraction of $(I, 0)$ is a dense open subset of ${\text{SO(3)}}\times\mathbb{R}^3$.

Proof. Using the third equation from (5.5) we have

We choose the following discrete-time Lyapunov candidate:

where $k_p > 0$ is a constant.

The stability of the attitude and angular velocity estimation error can be shown by analyzing $\Delta V_i = k_p\Delta\mathcal{U}_i + \Delta\mathcal{T}^l_i$.

Assuming $K_i$ to be constant and letting $K = K_i = K_{i+1}$ we obtain

where $\Delta Q_i = Q_{i+1} - Q_i$. Therefore,

Considering the first order expansion of $\exp{\left(\frac{h}{2}(\hat{\omega}_{i+1} + \hat{\omega}_i)^\times\right)}$ as

we have

It has to be noted that approximation in (5.10) is accurate for small values of $h$ and may affect the stability results for very high values of $h$. In the absence of measurement errors, we have $\tilde{U}_i^m = R_i^{\text{T}} E_i$. Therefore, it follows that

and noting that $L_i = E_iW_i(\tilde{U}_i^m)^{\text{T}}$, we get

Similarly we can compute the change in the kinetic energy as follows:

The change in the value of the candidate Lyapunov function can be computed as,

Similarly, we obtain

Substituting the value of $\omega_{i+2}$ from the filtering scheme presented in Proposition 1 we get

Now, for $\Delta V$ to be negative definite for all $i$ we require

where $l > 0, l \neq m$, and $\Delta V_{i+1}$ simplifies to

Substituting $\omega_{i+2}$ from the third equation presented in Proposition 1 into (5.18),

which further simplifies to,

which upon simple manipulations yields (5.4). It can be seen that after substituting (5.4) into Proposition 1, we obtain

(5.22) can also be rewritten as

in terms of $\omega_i, \omega_{i+1}$ and $S_{L_i}(\hat{R}_i)$. From (5.19), $\Delta V_i$ can be written as

We employ the discrete-time La-Salle invariance principle from ^[28] considering our domain (${\text{SO(3)}}\times\mathbb{R}^3$) to be a subset of $\mathbb{R}^{12}$. We use Theorem 6.3 and Theorem 7.9 from Chapter-1 of ^[28]. For this we first compute $\mathscr{E} : = \{(Q_i, \omega_i) \in {\text{SO(3)}}\times\mathbb{R}^3 | \Delta V_i(Q_i, \omega_i) = 0\} = \{(Q_i, \omega_i) \in {\text{SO(3)}}\times\mathbb{R}^3 \; | \; \omega_{i+1} + \omega_i = 0\}$. From (5.6), $\omega_{i+1} + \omega_i = 0$ implies that

Also, from (5.13) we have $\Delta\mathcal{U} = 0$ whenever $\omega_{i+1} + \omega_i = 0$. This implies that the potential function, which is a Morse function according to lemma 2.1, is not changing and therefore has converged to one of its stationary points. Stationary points of the Morse function $\langle I-Q, K \rangle$ are characterized as the solutions of

Now, $L_i = E_iW_i(\tilde{U}_i^m)^{\text{T}} = E_iW_i(R_i^{\text{T}} E_i)^{\text{T}} = (E_iW_iE_i^{\text{T}})R_i = KR_i$, which further gives us

Multiplying (5.27) from the right hand side by $\hat{R}_i^{\text{T}}$ and from the left hand side by $\hat{R}_i$,

At the critical points from (5.26), the right side of the above expression vanishes. Therefore, as $\hat{R}_i$ is an orthogonal matrix, the following holds true at the critical points:

Substituting this information in (5.23) yields,

Now if, $\omega_{i+1} + \omega_i = 0$, we have

This leads to the conclusion that the set of estimation errors, $\mathscr{E} = \{(Q_i, \omega_i) \in {\text{SO(3)}}\times\mathbb{R}^3 \; | \; Q_i \in C_Q, \omega_i = 0\}$, is the largest invariant set for the estimation error dynamics, and we obtain $\mathscr{M} = \mathscr{E} = \{(Q_i, \omega_i) \in {\text{SO(3)}}\times\mathbb{R}^3 \; | \; Q_i \in C_Q, \omega_i = 0\}$. Therefore, we obtain the positive limit set as the set,

In the absence of measurement errors, all the solutions of this filter converge asymptotically to the set $\mathscr{I}$. More specifically, the attitude estimation error converges to the set of critical points of $\langle\, I-Q, K\rangle$. The unique global minimum of this function is at $(Q, \omega) = (I, 0)$ from lemma 2.2. Therefore, $(Q, \omega) = (I, 0)$ is locally asymptotically stable. The remainder of this proof is similar to the last part of the proof of stability of the variational attitude estimator in ^[13].

Consider the set,

which consists of all the stationary states that the estimation errors may converge to, besides the desired estimation error state $(I, 0)$. Note that all states in the stable manifold of a stationary state in $\mathscr{C}$ will converge to this stationary state. From the properties of the critical points $Q^i \in C_Q \backslash (I)$ of $\mathit{\Phi}(\langle\, K, I-Q\rangle)$ given in lemma 2.2. we see that the stationary points in $\mathscr{I}\backslash(I, 0) = \{(Q^i, 0) : Q^i \in C_Q \backslash (I)\}$ have stable manifolds whose dimensions depend on the index of $Q^i$. Since the angular velocity estimate error $\omega$ converges globally to the zero vector, the dimension of the stable manifold $\mathcal{M}_i^S$ of $(Q^i, 0) \in {\text{SO(3)}}\times\mathbb{R}^3$ is

therefore, the stable manifolds of $(Q, \omega) = (Q^i, 0)$ are three-dimensional, four dimensional, or five-dimensional, depending on the index of $Q^i \in C_Q \backslash (I)$ according to (5.34). Moreover, the value of the Lyapunov function $V(Q_i, \omega_i)$ is non decreasing (increasing when $(Q_i, \omega_i) \notin \mathscr{I}$) for trajectories on these manifolds when going backwards in time. This implies that the metric distance between error states $(Q, \omega)$ along these trajectories on the stable manifolds $\mathcal{M}_i^S$ grows with the time separation between these states, and this property does not depend on the choice of the metric on ${\text{SO(3)}} \times \mathbb{R}^3$. Therefore, these stable manifolds are embedded (closed) sub-manifolds of ${\text{SO(3)}} \times \mathbb{R}^3$ and so is their union. Clearly, all states starting in the complement of this union, converge to the stable equilibrium $(Q, \omega) = (I, 0)$; therefore the domain of attraction of this equilibrium is,

which is a dense open subset of ${\text{SO(3)}}\times\mathbb{R}^3$.

6.   Numerical simulations

This section presents numerical simulation results of the discrete-time estimator presented in the section 5. The rigid body is assumed to have an initial attitude and angular velocity given by,

$W$ is selected based on the measured set of inertial vectors $E$ at each instant such that it satisfies lemma 2.1. Initially estimated states have the following initial estimation errors:

It has been assumed that there are at most 9 inertially known directions that are being measured by the sensors attached to the rigid body. The number of observed direction can vary randomly between 2 to 9 at each time instant. In the case where the number of observed directions is 2, the cross product of the two measurements is used as the third measurement. The standard rigid body dynamics are used to produce true states of the rigid body by applying sinusoidal forces. The observed directions in the body-fixed frame are simulated with the help of the aforementioned true states. The true quantities are disturbed by bounded, random noise with zero mean to simulate realistic measurements. Based on coarse attitude sensors like sun sensors and magnetometers, a random noise bounded in magnitude by $2.4^\circ$ is added to the matrix $U = R^{\text{T}} E$ to generate measured $U^m$. Similarly, a random noise bounded in magnitude by $0.97^\circ/s$, which is close to real noise levels of coarse rate gyros, is added $\Omega$ to generate measured $\Omega^m$.

6.1.   Case-1 simulation results

The estimator is simulated over a time interval of $T$ = 60s, with a step-size of $h = 0.01 s$. The inertial scalar gain is $m = 1.5$ and the dissipation term is chosen to be $l = 0.3$. The difference of sampling rate between measurements of angular velocity and measurements inertial vectors in body-fixed frame is taken to be $n = 10$. Furthermore, the value of gain $k_p$ is chosen to be $k_p = 1$. The principle angle $\phi$ of the rigid body's attitude estimation error $Q$ is shown in the Figure 1. Components of estimation error $\omega$ in the rigid body's angular velocity are shown in Figure 2.

6.2.   Case-2 simulation results

The estimator is simulated over a time interval of $T$ = 60s, with a step-size of $h = 0.01 s$. The inertial scalar gain is $m = 1.5$ and the dissipation term is chosen to be $l = 0.3$. The difference of sampling rate between measurements of angular velocity and measurements inertial vectors in body-fixed frame is to varying randomly but bounded between 10 and 30. Furthermore, the value of gain $k_p$ is chosen to be $k_p = 1$. The principle angle $\phi$ of the rigid body's attitude estimation error $Q$ is shown in the Figure 3. Components of estimation error $\omega$ in the rigid body's angular velocity are shown in Figure 4.

6.3.   Case-3 simulation results

The estimator is simulated over a time interval of $T$ = 60s, with a step-size of $h = 0.008 s$. The inertial scalar gain is $m = 2.5$ and the dissipation term is chosen to be $l = 0.5$. The assumption is that the angular velocity measurements are available after each time interval of 0.008s and inertial vector measurements are available after each time interval of 0.05s. Furthermore, the value of gain $k_p$ is chosen to be $k_p = 10$. The principle angle $\phi$ of the rigid body's attitude estimation error $Q$ is shown in the Figure 5. Components of estimation error $\omega$ in the rigid body's angular velocity are shown in Figure 6.

6.4.   Performance of the estimator in the presence of large initial estimaition errors

We have discussed the performance of the estimator in various scenarios but with the fixed values of initial estimation errors. It remains to see whether the estimator perform good in case of large initial estimation errors. For this part, we restrict our attention to the Case-2 i.e., when the difference of sampling rate between measurement of angular velocity and measurement of inertial vectors is time varying and bounded between 10 and 30. We check the performance of the estimator for three different values of initial estimation errors in angular velocity and attitude. The results of the same are plotted in Figure 7.

On comparing Figures 7a and 7b, when the initial estimation error in principle angle of attitude estimation is increased to about 2.5 times and estimation error in angular velocity is increased to 10 times, convergence takes about 10s more time. However, on comparing Figure 7a and 7c, it can be see that if the initial attitude estimation error is fixed then change in initial angular velocity estimation error(100 times the original values) has no effect on the convergence time. Therefore, it can be concluded that increasing the initial attitude estimation error can increase the convergence time but increase in initial angular velocity estimation error has no effect on the convergence time.

6.5.   Effect of gain on convergence time

In this section, we discuss how convergence time can be decreased by increasing the value of gain. Same as before, we restrict our attention to the Case-2 i.e., when the difference of sampling rate between measurement of angular velocity and measurement of inertial vectors is time varying and bounded between 10 and 30. We consider very high errors in initial estimation errors and plot the performance of the estimator for three different values of gain $k_p$ in Figure 8. As it can be seen, increasing the value of gain significantly improves the convergence time.

7.   Conclusion

We develop a geometric attitude and angular velocity estimation scheme using the discrete-time Lagrange-D'Alembert principle followed by discrete-time Lyapunov stability analysis in the presence of multi-rate measurements. The attitude determination problem from two or more vector measurements in the body-fixed frame is formulated as Wahba's optimization problem. To overcome the multi-rate challenge, a discrete-time model for attitude kinematics is used to propagate the inertial vector measurements forward in time. The filtering scheme is obtained with the aid of appropriate discrete-time Lagrangian and Lyapunov functions consisting of Wahba's cost function as an artificial potential term and a kinetic energy-like term that is quadratic in the angular velocity estimation error. As it can be observed, the Lyapunov function is not constructed from the same artificial potential and kinetic energy terms that are used for constructing the Lagrangian. There are mainly two reasons behind this; 1) the filtering scheme obtained by applying the discrete Lagrange-d'Alembert principle is implicit in nature and therefore it can increase computational load and runtime making it difficult to use for real-time applications. Therefore, an explicit filtering scheme is more desirable, 2) we also need the filtering scheme to be asymptotically stable. A Lyapunov function, different from system energy constructed appropriately, helps us meet both the requirements. The explicit filtering scheme obtained after the Lyapunov analysis was proven to be asymptotically stable in the absence of measurement noise and the domain of convergence is proven to be almost global. Numerical simulations were carried out with realistic measurement data corrupted by bounded noise. Numerical simulations verified that the estimated states converge to a bounded neighborhood of $(I, 0)$. Furthermore, the rate of convergence of the estimated states to the real state can be controlled by choosing appropriate gains. Future endeavors are towards obtaining a discrete-time optimal attitude estimator in the presence of multi-rate measurements when there a constant or slowly time-varying bias in the measurements of angular velocity while also obtaining a bound on the state estimation errors when there is measurement noise in the inertial vector measurements and the angular velocity measurements.

Conflict of interest

All authors hereby declare that there are no conflicts of interest in this paper.