State and Unknown Terrain Estimation for Planetary Rovers via Interval Observers

Herein, the problem of state and unknown terrain estimation is considered, where the unknown planetary terrain parameters, e.g., terrain stiffness and ground height, are inferred from how it affects rover motion through vehicle‐terrain interaction. In particular, an alternative framework for terrain estimation based on set‐valued or set‐membership estimation is proposed, where the goal is to find set‐valued estimates (in the form of hyper‐rectangles or intervals) of the states and unknown terrain parameters. For this purpose, a state and model interval observer is designed for partially unknown nonlinear systems with bounded noise. By leveraging a combination of nonlinear bounding/decomposition functions, affine abstractions, and a data‐driven function abstraction method (to overestimate the unknown dynamics model from noisy input–output data), the proposed observer is capable of simultaneously estimating the states and learning the unknown dynamics. Further, a tractable sufficient condition is derived for guaranteeing the stability of the designed observer, i.e., such that the sequence of interval estimate widths are uniformly bounded. When applied to the state and unknown terrain estimation problem, the simulation results indicate that our approach can more reliably find the range of possible terrain parameters when compared with the cubature Kalman filter.


Introduction
A recent survey on slippage estimation and compensation for planetary exploration rovers [1] (also see references therein) reveals that the entrapment and embedding of lunar and planetary exploration vehicles are often due to the lack of understanding of the lunar or planetary terrains. Thus, terrain properties estimation plays a critical role in rover mobility on unknown or uncertain terrains, especially in the context of lunar and planetary exploration. [1] By estimating or identifying terrain parameters, e.g., soil cohesion and internal friction angle, rovers (or off-road vehicles) can reduce the uncertainty in terrain understanding and predict longitudinal and lateral wheel slippage. [2] Moreover, the rover/vehicle control systems, e.g., traction control, stability control, assistive braking, and cruise control, as well as the path planning algorithms can adapt to the current/actual terrain properties [1,3] and as a result, improve their performance and reduce power consumption.

Literature Review
To complement efforts on terrain estimation in mobile robotics from exteroceptive sensors, e.g., light detection and ranging (LIDAR), vision and radar, [4][5][6] several recent works have considered terrain estimation using proprioceptive sensors based on vehicle-terrain (or more precisely, wheel-terrain) interactions. [2,3,[7][8][9] In fact, a related problem of wheel-road friction coefficient estimation has been extensively studied in the automotive field and those approaches are already being commercialized. Nonetheless, these methods are often unsuitable for off-road driving and uncertain planetary terrains. [3] To address this issue, several wheel-terrain interaction approaches for off-road and planetary terrain estimation have been proposed in previous studies [2,3,7,8] that are based on linear regression, discrimination analyses, and Kalman filtering. Specifically, methods based on state and parameter estimation via Kalman filtering have been shown to provide optimal solution to estimation problems in the automotive and mobile robotics fields when the system dynamics are linear and the noise is Gaussian-distributed, [3,10] whereas extended Kalman filters (EKF) and unscented Kalman filters (UKF) have been used when the models are nonlinear for estimating vehicle states and parameters such as side-slip angles. [11][12][13] More recently, methods based on cubature Kalman filter (CKF) [14] have also been proposed to further improve the terrain estimation performance. [2] However, these approaches often hinge on the idea that the disturbance inputs caused by terrain irregularity can be modeled as stochastic uncertainties with known statistics, e.g., as a zeromean Gaussian noise.
In contrast, when bounds on the environmental uncertainties or sensor errors are available (instead of distributional information), set-valued and interval observers have been designed for DOI: 10.1002/aisy.202100040 Herein, the problem of state and unknown terrain estimation is considered, where the unknown planetary terrain parameters, e.g., terrain stiffness and ground height, are inferred from how it affects rover motion through vehicleterrain interaction. In particular, an alternative framework for terrain estimation based on set-valued or set-membership estimation is proposed, where the goal is to find set-valued estimates (in the form of hyper-rectangles or intervals) of the states and unknown terrain parameters. For this purpose, a state and model interval observer is designed for partially unknown nonlinear systems with bounded noise. By leveraging a combination of nonlinear bounding/decomposition functions, affine abstractions, and a data-driven function abstraction method (to overestimate the unknown dynamics model from noisy input-output data), the proposed observer is capable of simultaneously estimating the states and learning the unknown dynamics. Further, a tractable sufficient condition is derived for guaranteeing the stability of the designed observer, i.e., such that the sequence of interval estimate widths are uniformly bounded. When applied to the state and unknown terrain estimation problem, the simulation results indicate that our approach can more reliably find the range of possible terrain parameters when compared with the cubature Kalman filter. multiple classes of systems, including linear time-invariant (LTI), [15] linear parameter-varying (LPV), [16,17] Metzler and/or partial linearizable, [18,19] cooperative, [18,20] Lipschitz nonlinear, [21] monotone nonlinear, [22,23] and uncertain nonlinear [24] systems. Moreover, the set-valued observer design framework was recently extended to additionally estimate arbitrary unknown disturbance inputs for LTI, [25] LPV, [26] switched linear, [27] and nonlinear [28,29] systems with bounded-norm noise.
The aforementioned state and parameter estimation approaches often assume that a mathematical dynamic model of the system of interest is known. Hence, when the system model is not exactly known, the unknown dynamics needs to be estimated or learned from data, where learning methods such as support vector machines (SVMs), Bayesian networks and deep learning have been used to model and learn terrain properties. [1,30,31] In particular, a learning/identification technique known as nonlinear set membership prediction, Lipschitz interpolation or kinky inference [32][33][34][35][36] that can capture worst-case generalization errors holds great promise for being integrated into the set-valued estimation framework that we consider in this article for terrain estimation.

Contribution
This article presents a novel approach for state and unknown terrain estimation based on wheel-terrain interaction using the framework of a nonlinear interval observer/estimator for partially known dynamical systems with bounded noise. Specifically, the proposed approach bridges between setmembership state estimation approaches, e.g., in the studies by Mazenc and Bernard, Rassi et al., Efimov et al., Khajenejad and Yong,[15,18,21,26,27] and data-driven function approximation methods, e.g., in the studies by Milanese and Novara, Canale et al., Zabinsky et al.,Beliakov,Calliess,[32][33][34][35][36] to design intervalvalued observers for nonlinear dynamical systems with bounded noise and partially known dynamics, where the state and observation vector fields belong to a fairly general class of nonlinear functions and the partially known system dynamics contains an unknown function, which can represent dynamic terrain uncertainty.
Similar to stochastic filtering, our recursive interval observer approach consists of a state propagation step to construct framers of system states and unknown parameters (i.e., upper and lower bounds that contain/sandwich the true states and parameters) that are consistent with the dynamic model by leveraging techniques from nonlinear decomposition/bounding functions [29,[37][38][39] and affine abstractions, [40] as well as a measurement update step to tighten the framers based on sensor measurements. In addition, as the system dynamics is only partially known, we adopt the data-driven abstraction approach based on the study by Jin et al. [41] to recursively over-approximate the unknown dynamics function from noisy observation data and interval estimates from the update step. Furthermore, we provide a sufficient condition for stability of our observer, which if satisfied, guarantees that the sequence of interval estimate widths are uniformly bounded. Finally, the effectiveness of the proposed observer is demonstrated on the problem of state and unknown terrain estimation, and our approach is observed to outperform methods based on CKF. Note that a preliminary version of the interval observer design appeared in a conference proceeding. [42] The article is organized as follows. Section 2 describes the wheel-terrain interaction model and formulates a general interval observer design problem for discrete-time partially known dynamical systems with bounded noise. Section 3 designs a state and model interval observer (SMIO) and proves the correctness and stability of our design. Simulation results are presented in Section 4 to validate our algorithm and Section 5 concludes the article.

Preliminaries and Problem Formulation
In this section, we first present the wheel-terrain interaction model we will be using for state and unknown terrain estimation. Then, after presenting some background mathematical ideas, we formally formulate the SMIO design problem whose solution will later be applied for terrain estimation in the simulation section.

Wheel-Terrain Interaction Modeling
Similar to the study by Reina et al., [2] the wheel-terrain interaction model we consider in this article is based on the standard quarter-car model, which consists of two mass bodies, i.e., the vehicle sprung mass m s and unsprung mass m ns , with lumped stiffness and viscous friction parameters given by k and c, as shown in Figure 1. The vertical displacement of m s and of m ns are represented by y 1 and y 2 , respectively. Furthermore, we model the wheel-terrain (or tire-soil) interaction as a spring with a combined stiffness k tot that represents the equivalent stiffness corresponding to soil deformability k s and tire stiffness k t , given by where the combined stiffness k tot is unknown, while the parameters m s , m ns , k, and c are known. The equations of motion for the quarter-car model are given by the following 2 Þ þ kðy 1 À y 2 Þ ¼ 0 m nsÿ2 þ cðy : 2 À y : 1 Þ þ kðy 2 À y 1 Þ þ k tot ðy 2 À hÞ ¼ 0 (2) where h is the terrain elevation profile, whose unknown variation is assumed to satisfy an unknown functionf , as follows to capture the variation of the terrain height, where v is the forward velocity that is assumed to be known or measured. Moreover, as k tot is unknown but constant, we adopt a common approach in state and parameter estimation to describe the constant k tot as a state with zero dynamics as follows Then, defining the state vector as we find the resulting equations of motion for the quarter-car and wheel-terrain interaction model as Due to the presence of disturbance/noise signals and unknown dynamics, the states of the system in (8) cannot be always measured directly, hence need to be estimated. Consequently, the problem of state and unknown terrain estimation can be considered as a specific instance of a more general problem of designing interval observers for partially unknown bounded-error nonlinear systems in the form of (8), which will be described Section 3.

Preliminary Material
Before formulating the problem, we first introduce some notations, definitions, and related results that will be useful throughout the article.
Notation. ℝ n denotes the n-dimensional Euclidean space and ℝ þþ positive real numbers. For vectors v, w ∈ ℝ n and a matrix M ∈ ℝ pÂq , kvk ≜ ffiffiffiffiffiffiffi v ⊤ v p and kMk denote their (induced) 2-norm, and v ≤ w is an elementwise inequality. Moreover, the transpose, Moore-Penrose pseudoinverse, ði, jÞth element and rank of M are given by M ⊤ , M † , M i,j and rkðMÞ, whereas M ðr∶sÞ is a submatrix of M, consisting of its rth through sth rows, and its row support is r ¼ rowsuppðMÞ ∈ ℝ p , where r i ¼ 0 if the ith row of M is zero and r i ¼ 1 otherwise, ∀i ∈ f1 : : : We call M a non-negative matrix, i.e., M ≥ 0, if M i,j ≥ 0, ∀i ∈ f1 : : : pg, ∀j ∈ f1 : : : qg.
Definition 1 (Interval, Maximal and Minimal Elements, Interval Width). An (multidimensional) interval ℐ ⊂ ℝ n is the set of all real vectors x ∈ ℝ n that satisfies s ≤ x ≤ s, where s, s, and ks À sk are called minimal vector, maximal vector, and width of ℐ, respectively.
As a corollary, if A is non-negative, then Ax ≤ Ax ≤ Ax.

Definition 3 (Mixed-Monotone Mappings and Decomposition Functions). [37, Definition 4] A mapping f
Note that the decomposition function of a vector field is not unique and a specific one is given in the study by Yang et al. [37,Theorem 2]: If a vector field q ¼ ½q ⊤ 1 : : : q ⊤ n ⊤ ∶ X ⊆ ℝ n ! ℝ m is differentiable and its partial derivatives are bounded with known bounds, i.e.,  13)]. Consequently, for x ¼ ½x 1 : : : x j : : : x n ⊤ , y ¼ ½y 1 : : : y j : : : y n ⊤ , we have In contrast, when the precise lower and upper bounds, a i,j , b i,j , of the partial derivatives are not known or are hard to compute, we can obtain upper and lower approximations of the bounds using Proposition 3 with the slopes set to zero, or by leveraging interval arithmetic. [43]

Problem Formulation
Consider a partially unknown nonlinear discrete-time system in the following form (cf. (8) and (9) for the terrain estimation example) where x k ∈ X ⊂ R n is the state vector at time k ∈ ℕ, u k ∈ U ⊂ R m is a known input vector, y k ∈ R l is the measurement vector and d k ∈ D ⊂ R p is a dynamic input vector whose dynamics is governed by an unknown vector field ϕð⋅Þ (Note that if the vector field ϕð⋅Þ is partially known (i.e., consists of the sum of a known componentφð⋅Þ and an unknown componentφð⋅Þ), we can simply consider d kþ1 Àφð⋅Þ as the output data for the model learning procedure to learn a model of the (completely) unknown functionφð⋅Þ). It might be worth emphasizing that d k can be considered as a dynamic exogenous input, whose dynamics/changes can be described by the unknown function ϕð⋅Þ, which is an unknown mapping from the current values of state, x k , known input u k , exogenous input d k , and disturbance w k to the exogenous input value d kþ1 at the next time step.
Moreover, we refer to z k ≜ ½x ⊤ k d ⊤ k ⊤ as the augmented state. The process noise w k ∈ R n w and the measurement noise v k ∈ R l are assumed to be bounded, with w ≤ w k ≤ w and v ≤ v k ≤ v, where w, w and v, v are the known lower and upper bounds of the process and measurement noise signals, respectively. We also assume that lower and upper bounds, z 0 and z 0 , for the initial augmented state ∀j ∈ f1 : : : pg is known to be Lipschitz continuous. For simplicity and without loss of generality, we assume that the Lipschitz constant L ϕ j is known; otherwise, we can estimate the Lipschitz constants with any desired precision using the approach in the study by Singh et al. [41,Equation (18) and Proposition 3], which is briefly recapped in Section 3.5. Moreover, we assume the following: k ⊤ and the known inputs u k , ∀k ∈ f0 : : : ∞g, respectively.
The observer design problem for a general partially unknown nonlinear discrete-time system in the form of (8) can be cast as follows: Problem 1. Given a partially known nonlinear discrete-time system (8) with bounded noise signals and unknown dynamics ϕð⋅, ⋅ , ⋅ , ⋅Þ, design a stable observer that simultaneously finds bounded intervals of compatible augmented states, k ⊤ , and learns an unknown dynamics model for ϕð⋅, ⋅ , ⋅ , ⋅Þ.
the quarter-car and wheel-terrain interaction models described in Section 2 for simultaneously estimating system states and unknown terrain parameters, where the combined soil-tire stiffness and the terrain elevation variation function are unknown.

Why Set-Valued Observers?
As mentioned in Section 1, many state and parameter estimation approaches via a Kalman filtering framework have been proposed for terrain estimation based on wheel-terrain interaction, e.g., the studies by Reina and coworkers. [2,3,10] Implicit in this Kalman filtering framework are the assumptions that the disturbance signals caused by terrain irregularity can be modeled as stochastic uncertainties with known statistics (e.g., as a zeromean Gaussian noise, so that the available statistical information on the model uncertainty can be leveraged in the design process) and that the system dynamics is fully known. By contrast, this article considers the setting where these assumptions are not available or justified. Instead, we consider a distribution-free framework that treats disturbances/noise as nondeterministic and bounded signals, and we assume that the system dynamics is only partially known. Under this setting (with a lack of additional statistical assumptions), Kalman filtering approaches are unfortunately not directly applicable. In contrast, set-valued and interval observers are viable and appropriate approaches when only bounds of the disturbance/noise signals are available and also when only distribution-free bounds on the generalization errors of the model learning methods can be derived.

Recursive Interval Observer
In this section, we briefly recap a three-step recursive interval observer that was introduced in our previous work [42] to solve the state and unknown terrain estimation problem described in Section 2.1 that is modeled as a partially known system (8). The main idea is to combine a data-driven model learning approach to learn the unknown function ϕð⋅Þ and a model-based interval observer approach to estimate the augmented states, z k ≜ ½x ⊤ k d ⊤ k ⊤ (consisting of the state and the exogenous input). Three main constituents are combined to design the observer structure: a state propagation (SP), a measurement update (MU) step, and a model learning (ML) step. Via the state propagation step, the interval estimate for the augmented states is propagated for one time step through the nonlinear state equation and the estimated model of the unknown dynamics function obtained in previous time step. Then, the update step iteratively updates compatible intervals of the augmented states, given new measurements and the nonlinear observation function, and finally, the upper and lower framer functions (abstractions) for the unknown dynamics function are estimated in the model learning step. Mathematically speaking, the three observer steps can be described in the following form (with where ℱ p and ℱ u are to-be-designed interval-valued mappings and ℱ l a to-be-constructed function over-approximation procedure (abstraction model), whereas ℐ z p k and ℐ z k are the intervals of compatible propagated and estimated augmented states, respectively. Moreover, fϕ k ð⋅Þ, ϕ k ð⋅Þg is a data-driven abstraction or over-approximation model for the unknown function ϕð⋅Þ, at time step k, i.e.
with D ϕ being the domain of ϕð⋅Þ and ζ k ≜ ½z ⊤ k u ⊤ k w ⊤ k ⊤ . To avoid the computational complexity of optimal observers, [44] while taking the advantages of intervals, [17] the following form of interval estimates in the propagation and update steps is considered where the estimation would be equivalent to finding the maximal and minimal values of ℐ z p k and ℐ z k , i.e., z p k , z p k , z k , z k . Further, given the sequence of interval estimates up to the current time, at the model learning step, we aim to apply the data-driven function abstraction/over-approximation approach developed in our previous work [41] to update and refine the learned/estimated model of the unknown dynamics function ϕð⋅Þ at the current time step.
In particular, defining the augmented state our proposed interval observer at each time step k ∈ ℕ is given as follows: State Propagation Measurement Update x Model Learning ϕ k,j ðζ k Þ ¼ min t∈f0, : : : , TÀ1g ϕ k,j ðζ k Þ ¼ max t∈f0, : : : , TÀ1g d kÀt,j À L ϕ j kζ k Àζ kÀt k À ε j kÀt (24) www.advancedsciencenews.com www.advintellsyst.com where j ∈ f1 : : : pg, fζ kÀt ¼ 1 2 ðζ kÀt þ ζ kÀt Þg k t¼0 , and fd kÀt , d kÀt g k t¼0 are the augmented input-output data set. At each time step k, the estimated framers gathered from the initial to the current time step construct the augmented data set, which is used in the model learning step to recursively derive over-approximations of the unknown function ϕð⋅Þ, i.e., fϕ k ð.Þ, ϕ k ð.Þg by applying [41,Theorem 1]. In addition, the propagated state framers at time step k are computed as Moreover, the sequences of updated framers fz u i,k , z u i,k g ∞ i¼1 are iteratively computed as follows where kÀt , ∀q ∈ ff , ϕg, J ∈ fA, Wg, i ∈ f1, : : : , ∞g, j ∈ f1, : : : , pg, are to-be-designed observer parameters, matrix gains of appropriate dimensions at time k and iteration i (given in Theorem 1 and Appendix), whereas f d ð., ., ., .Þ is the bounding function (based on (10)), with the purpose of achieving desirable observer properties.
Note that the measurement update step is done iteratively, because the tightness of the upper and lower bounding functions for the observation function g (cf. Propositions 2 and 3) depends on the a priori interval ℬ, (see proof of Theorem 2 for more explanation). Hence, if tighter updated intervals are obtained starting from the compatible intervals from the propagation step, we can use them as the new ℬ to obtain better abstraction/ bounding functions for g, which in turn may lead to even tighter updated intervals. Repeating this process results in a sequence of monotonically tighter updated intervals, that is convergent by the monotone convergence theorem, and its limit is chosen as the final interval estimate at time k. It is worth mentioning that in practice, the iterations can be terminated when the improvement is less than a user-specified stopping/threshold criterion, without jeopardizing the correctness and stability properties we will prove in the next sections.
Further, in the model learning step with the history of obtained compatible intervals up to the current time, given f½z s , z s g k s¼0 as the noisy input data and the compatible interval of unknown inputs, ½d k , d k , as the noisy output data, by leveraging our previous result in the study by Jin et al. [ is recursively constructed for the unknown function ϕð⋅Þ, that by construction, satisfies (45). In other words, it is guaranteed that our model estimation is correct (i.e., is guaranteed to frame/bracket the true function) and becomes more precise with time (cf. Lemma 1).

Correctness of the Observer
The objective of this section is to guarantee the framer property [19] of the proposed SMIO observer by designing the observer gains. In other words, we desire to make sure that the observer returns correct interval estimates, in the sense that starting from the initial interval z 0 ≤ z 0 ≤ z 0 , the true augmented states of the dynamic system (8) are guaranteed to be within the estimated intervals, given by (18)-(28). We call fz k , z k g ∞ k¼0 an augmented state framer sequence for system (8), if the observer is correct. Prior to derive our result on correctness of the observer, we state a modified version of our previous result in the study by Singh et al. [40,Theorem 1], in a unified manner that enables us to derive parallel global and local affine bounding functions for our known f ð⋅Þ, gð⋅Þ and unknown ϕð⋅Þ vector fields.
Proposition 3 (Parallel Affine Abstractions). Let the entire space be defined as X and suppose that Assumption 2 holds. Consider the vector fields qð.Þ, qð.Þ∶X ⊂ ℝ n 0 ! ℝ m 0 , where ∀ζ ∈ X, qðζÞ ≤ qðζÞ, along with the following Linear Program (LP) where ℬ is an interval with ζ, ζ and V ℬ being its maximal, minimal and set of vertices, respectively, 1 m ∈ ℝ m is a vector of ones, σ q is given in the study by Singh et al. A q ζ þ e q ≤ qðζÞ ≤ qðζÞ ≤ A q ζ þ e q , ∀ζ ∈ X.
( 3 2 ) Using the aforementioned proposition, we first solve (30) on the entire space X, i.e., with ℬ ¼ X (where the constraint (31) is trivially satisfied and is thus redundant) and obtain a tuple of ðθ q , A q , e q , e q Þ that satisfies (32), i.e., we construct a global affine abstraction model for the pair of functions qð.Þ, qð.Þ on the entire space X.
Next, given the (global) tuple ðA q , e q , e q Þ computed as described earlier, we solve (30) on ℬ subject to (31) to obtain a tuple of local parallel affine abstraction matrices for the pair of functions fqð⋅Þ, qð⋅Þg on the interval ℬ, satisfying the following: ∀ζ ∈ ℬ www.advancedsciencenews.com www.advintellsyst.com Now, equipped with all the required tools, we state our first main result on the framer property of the SMIO observer.
Theorem 1 (Correctness of the Observer). Consider the system (8) with its augmented state defined as z ≜ ½x ⊤ d ⊤ ⊤ , along with the SMIO observer in (18)- (24). Suppose that Assumptions 1-2 hold, f d ð⋅Þ is a decomposition function of f ð⋅Þ and observer gains and parameters are designed as given in (6.1). Then, the SMIO observer estimates are correct, i.e., the sequences of intervals fz k , z k g ∞ k¼0 are framers of the augmented state sequence of system (8) that satisfy z k ≤ z k ≤ z k for all k.
Proof. We will prove this by induction. For the base case, by assumption, z 0 ≤ z 0 ≤ z 0 holds. Now, for the induction step, suppose that z kÀ1 ≤ z kÀ1 ≤ z kÀ1 . Then, Propositions 1-3 as well as (8), (20)- (25) and the study by Singh et al. [41,Theorem 1] imply that z p k ≤ z k ≤ z p k . Given this, iteratively obtaining upper and lower abstraction matrices for the observation function gð.Þ based on Proposition 3 and applying Proposition 1, we have where α i,k , α i,k are given in (29) and A g i,k is a solution of the LP in (30), i.e., the parallel abstraction slope for function gð⋅Þ at iteration i in the corresponding compatible interval ½z u iÀ1,k , z u iÀ1,k . Then, multiplying (34) by A g † i,k , Proposition 1 and using the fact that z u iÀ1,k , z u iÀ1,k are framers for the augmented state z k at time k and, [45] we obtain z u i,k ≤ z k ≤ z u i,k , with z u i,k , z u i,k given in (27). Now, note that by construction, the sequences of updated upper and lower framers, fz u i,k g ∞ i¼0 and fz u i,k g ∞ i¼0 with z u 0,k ¼ z p k and z u 0,k ¼ z p k , are monotonically decreasing and increasing, respectively, and hence are convergent by the monotone convergence theorem. Consequently, their limits z k , z k are the tightest possible framers, i.e., ∀i ∈ f1 : : : ∞g i,k ≤ : : : ≤ z u i,k ≤ : : : ≤ z u 0,k ≤ : : : ≤ z u i,k ≤ : : : where z k , z k are the returned updated augmented state framers by the observer. This completes the proof. Next, through the following lemma, we show that the abstraction model of the unknown dynamics function becomes tighter (i.e., more precise) over time given correct interval estimates. Hence, our model estimate of the unknown dynamics becomes more accurate as time increases. Lemma 1. Consider the system (8) and the SMIO observer in (18)- (28) and suppose that all the assumptions in Theorem 1 hold. Then, the following sequence of inequalities holds ϕ 0 ðζ 0 Þ ≤ : : : ≤ ϕ k ðζ k Þ ≤ : : that is, the unknown input model estimations/abstractions are correct and become more precise or tighter with time.
Proof. It directly follows from the study by Singh et al. [41,Theorem 1] and Theorem 1 that the model estimates are correct, i.e., ∀k ∈ f0 : : : ∞g∶ϕ k ðζ k Þ ≤ ϕðζ k Þ ≤ ϕ k ðζ k Þ. Moreover, considering the data-driven abstraction procedure in the model learning step, note that by construction, the data set used at time step k is a subset of the one used at time k þ 1. Hence, from the study by Singh et al. [41,Proposition 2] the abstraction model satisfies monotonicity, i.e., (45) holds.

Observer Stability
This section investigates the stability of the designed observer, which will be first formally defined as follows Definition 4 (Observer Stability). The observer SMIO (18)-(24) is stable, if the sequence of interval widths fkΔ z kÀ1 k ≜ kz kÀ1 À z kÀ1 kg ∞ k¼1 is uniformly bounded, and consequently, the sequence of estimation errors fk˜z kÀ1 k ≜ maxðkz kÀ1 À z kÀ1 k, kz kÀ1 À z kÀ1 kÞ is also uniformly bounded.
Remark 1. Note that it is straightforward to observe that the above notion of (global) stability implies that the estimation error system is uniformly bounded-input bounded-state (UBIBS), which is a widely used stability notion for distribution-free approaches. [46][47][48][49][50][51] In particular, a dynamic system is UBIBS, if bounded initial states x 0 and bounded (disturbance/noise) inputs u produce uniformly bounded trajectories [46, Section 3.2], i.e., there exist two К-functions (A function σ∶ℝ þ ! ℝ þ is a К-function if it is continuous, strictly increasing and σð0Þ ¼ 0 σ 1 and σ 2 ) such that Next, a useful property for the decomposition function given in (10) will be derived, which will be helpful in obtaining sufficient conditions for the observer stability.
We are now ready to state our next main result on the SMIO observer stability in the following theorem.
Theorem 2 (Observer Stability). Consider the system (8) along with the SMIO observer in (18)- (28). Let D m be the set of all diagonal matrices in ℝ mÂm with their diagonal arguments being 0 or 1. Suppose that all the assumptions in Theorem 1 hold and the decomposition function f d is constructed using (10). Then, the observer is stable if there exist D 1 ∈ D nþp , D 2 ∈ D l , D 3 ∈ D n that satisfy D 1,i,i ¼ 0 if rðiÞ ¼ 1, i.e., if there exist (40) such that (10). Proof. Note that, we aim to obtain sufficient stability conditions that can be checked a priori instead of for each time step k. In contrast, for the implementation of the update step, we iteratively find new local parallel abstraction slopes A g i,k by iteratively solving the LP (35) for g on the intervals obtained in the previous iteration, (29)), with additional constraints given in (36) in the optimization problems, which guarantees that the iteratively updated local intervals obtained using the local abstraction slopes are inside the global interval, i.e.
z u k ≤ z u 0,k ≤ : : : ≤ z u i,k ≤ : : and ðD Ã 1 , D Ã 2 , D Ã 3 Þ is a solution of the following problem Consequently, the sequence of interval widths fkΔ z k kg ∞ k¼1 is uniformly upper bounded by a convergent sequence as Proof. The proof is straightforward by applying Proposition 1, computing (47) iteratively, using the fact that by Theorem 2, AðD 1 , D 2 , D 3 Þ is a stable matrix for the tuple of ðD 1 , D 2 , D 3 Þ that is a solution of (41), and from triangle inequality.

Estimation of Lipschitz Constant
In previous sections, the Lipschitz constants are assumed to be given. In the case when the constants are not known, they can be estimated from the noisy sampled data set D ¼ fðs j ,ỹ jþ1 Þj j ¼ n y , · · · , N À 1g as followŝ whereỹ k ands k are the augmented measured/estimated data input and data output vectors of unknown function ϕð⋅Þ at time step k. This can be considered as an extension of the lazy update rule in the study by Calliess [36,Section 4.3.2] to the case where both the input and output data, i.e.,s j andỹ jþ1 for all j, are corrupted by bounded noise. The expression in (62) simply follows from the definition of Lipschitz continuity and the use of triangle inequality.
To guarantees the accuracy of L ðiÞ p that is critical for the results in the previous sections, we proceed to find some confidence that we obtain the right estimate with high probability. To achieve this goal, we leverage a classical result on probably approximately correct (PAC) learning for linear separators, which is summarized as follows Definition 5 (Linear Separators introduction). [53] For Γ ⊂ ℝ Â ℝ, a linear separator is a pair ða, bÞ ∈ ℝ 2 such that ∀ðx, yÞ ∈ Γ∶ x ≤ ay þ b (63) Proposition 4 (PAC Learning introduction). [53] Let ε, δ ∈ ℝ þ . If number of sampling points is N ≥ 1 ε ln 1 δ , where the sample points Γ are drawn from a distribution , then, with probability greater than 1 À δ, a linear separator ða, bÞ has an error of less than ε, where the error of a pair ða, bÞ is defined as .
From Definition 5, it is easy to verify that ourL ðiÞ P estimate in (62) is a special case of the linear separator in (63) with b ¼ 0. Thus, the estimatedL ðiÞ P using (62) is guaranteed to be close to the true Lipschitz constant of the original unknown function with high probability if we have sufficient data.
However, from the nature of the Equation (62), the estimated value of the Lipschitz constantL p tends to be smaller than the true value of the Lipschitz constant L p , especially when the available data is limited, which may lead to the violation of the framer property of our proposed interval observer. To overcome this potential concern, we adopt a heuristic approach to estimate the Lipschitz constant. First, we initialize the algorithm with a prior sampled data set D with size n 0 , a large constant L p , and a desired data size N 0 which can be determined based on Proposition 4. Then, we expand the data set D with newly observed data and estimateL p using (36) for each time step. The Lipschitz constant is selected by the following equation In Equation (64), if the size of the current data set n 0 satisfies the desired data size N 0 , the estimated Lipschitz constantL p will be kept, otherwise, the large constant L p is used for compensation. Our proposed SMIO with the Lipschitz constant estimation heuristic is summarized in Algorithm 1. ðz k , z k Þ ¼ ðz u ∞,k , z u ∞,k Þ; ℐ z k ¼ fz ∈ R n ∶z k ≤ z ≤ z k g;

Simulation Results
In this section, we apply the proposed SMIO in Section 3 to the quarter-car and wheel-terrain interaction models described in Section 2 for simultaneously estimating system states and unknown terrain parameters, where the combined soil-tire stiffness and the terrain elevation variation function are unknown.
In particular, we compare the estimation performance of the proposed SMIO is compared with the performance of CKFs, which was the estimator of choice in the study by Reina et al., [2] under the same operating conditions and some additional statistical assumptions to enable the use of the CKF. The function trackingCKF in MATLAB is used to implement the CKF algorithm and interested readers are referred to the literature, e.g., the studies by Reina et al. and Arasaratnam and Haykin, [2,14] for more details. In addition, we verify that stability condition of the interval observers in Theorem 2 is indeed satisfied in this example. The parameters of a typical off-road vehicle (representing a planetary rover) that are used in the simulations are shown in Table 1. The sampling time is chosen to be δt ¼ 0.01s and the true terrain stiffness k s is selected to be k s ¼ 651.1 kN, thus, the true unknown k tot ¼ k s k t k s þk t ¼ 137.93 kN m À1 . Further, the true unknown terrain elevation variation function is chosen to be h : ¼f ðh, v, tÞ ¼ sinðhÞ þ 3 sinðvtÞ with v ¼ 5 m s À1 , which is treated as a Gaussian-distributed noise in the CKF.
We initialize the simulations with the initial unknown function model obtained with 1000 sample points and choose the compensation constant L p to be 10. The initial augmented state bound of z is given by Moreover, the process and measurement noise signals are assumed to be uniformly distributed with the following bounds: jw 1,k j ≤ 0.1897, jw 2,k j ≤ 0.1897, jw 3,k j ≤ 0.0949, jv 1,k j ≤ 2.1213, jv 2,k j ≤ 2.1213, jv 3,k j ≤ 0.001, jv 4,k j ≤ 0.001, jv 5,k j ≤ 0.001, jv 6,k j ≤ 0.001, and jv 7,k j ≤ 0.001.
To compare our method with the CKF approach, we relate the noise bounds to the variance matrix by setting jw i,k j=jv i,k j ¼ 3σ w i,k =3σ v i,k , with a starting guess for the CKF as follows As shown in Figure 2-7, the state estimates obtained by the proposed SMIO performed better than the CKF in terms of the estimation errors, and the difference is particularly notable for the estimate of the combined tire-soil stiffness where the percentage error is 6.5% for SMIO and 34% for the CKF approaches, respectively. Moreover, SMIO can correctly identify bounds/ framers for the states, whereas CKF does not provide such guarantees and appears to be unstable. The poor performance of the CKF method can be a result of treating the unknown dynamics function and uniformly distributed noise as Gaussian distributions. In Figure 8, the percentage error for the terrain stiffness estimate is 41.1% for SMIO and 73.7% for CKF. The error is relatively large because the nonlinear function that outputs the k s amplifies the estimation error of k tot .
Further, to verify that the correctness of our approach, we first obtain finite-valued upper and lower bounds (horizontal abstractions) for the partial derivatives of f ð⋅Þ using Proposition 3 with abstraction slopes set to zero, as follows    www.advancedsciencenews.com www.advintellsyst.com are correct. This is shown in Figure 2-9, where the true states and unknown inputs as well as interval estimates are depicted. In addition, solving the optimization problem in Proposition 3 for the global abstraction matrices, we obtain A ϕ ¼ ½1.14 À0.40 (72) Figure 7. Actual value, estimates, and estimation errors of k tot when using our proposed SMIO versus the CKF.  www.advancedsciencenews.com www.advintellsyst.com Next, from the study by Yang et al. [37, (10)-(13)]), we obtain C f ¼ ½0 5Â6 when using (10). Consequently, (47) is satisfied and so, the sufficient condition in Theorem 2 holds. Thus, as expected, we obtain uniformly bounded and convergent interval estimate errors when applying our observer design, as can be seen in Figure 10, where at each time step, the actual error sequence is upper bounded by the interval widths, which converge to steady-state values. Further, Figure 9 (right) shows the framer intervals of the learned/estimated unknown dynamics model (depicted by the "kinky" red and blue meshes) that frame the actual unknown dynamics function ϕð⋅Þ, as well as the global abstraction that is computed via Proposition 3 at the initial step. In summary, the simulation example demonstrates that the proposed SMIO is effective for state and unknown terrain estimation and holds great promise for enhancing the mobility of lunar and planetary vehicles.

Conclusion
In this article, we consider the design of a novel terrain estimation approach based on the framework of an SMIO. Specifically, an interval observer was introduced for partially unknown nonlinear systems with bounded noise. By leveraging a combination of nonlinear bounding or decomposition functions, affine abstractions, and a data-driven function abstraction method (to overestimate the unknown dynamics model from noisy inputoutput data), the proposed observer is capable of estimating the augmented states and learns the unknown dynamics simultaneously. In addition, we derived a tractable sufficient condition for the stability of the designed observer, i.e., for the uniform boundedness of the sequence of interval estimate widths. Finally, the effectiveness of the proposed interval observer is demonstrated using a terrain estimation problem, where we observed that the estimation performance is better when compared with the use of a CKF.
Future research on this subject includes the use of more realistic lunar or planetary rover models, e.g., the rocker-bogie system, and the fusion of the estimated parameters with estimates from exteroceptive sensors to further enhance the estimation of lunar and planetary terrain properties.