Robust Nonlinear Tracking Control with Exponential Convergence Using Contraction Metrics and Disturbance Estimation

This paper presents a tracking controller for nonlinear systems with matched uncertainties based on contraction metrics and disturbance estimation that provides exponential convergence guarantees. Within the proposed approach, a disturbance estimator is proposed to estimate the pointwise value of the uncertainties, with a pre-computable estimation error bounds (EEB). The estimated disturbance and the EEB are then incorporated in a robust Riemannian energy condition to compute the control law that guarantees exponential convergence of actual state trajectories to desired ones. Simulation results on aircraft and planar quadrotor systems demonstrate the efficacy of the proposed controller, which yields better tracking performance than existing controllers for both systems.


Introduction
Robotic systems generally have nonlinear dynamics and are subject to model uncertainties and disturbances. Moreover, many robotic systems are underactuated, i.e., having fewer independent control inputs than degrees of freedom, including fixed-wing aircraft, quadrotors and dynamic walking robots. The design of tracking controllers for underactuated robotic systems is a much more challenging problem compared to that for fully-actuated systems. Recently, the concept of a control contraction metric (CCM) was introduced in [1] to synthesize trajectory tracking controllers for general nonlinear systems, including underactuated ones. The CCM extends contraction theory [2] from analysis to constructive control design, while contraction theory is focused on analyzing nonlinear systems in a differential framework by studying the convergence between pairs of state trajectories toward each other. It was shown in [3] that CCM reduces to conventional sliding and energy-based designs for fully-actuated systems. On the other hand, for underactuated systems, compared to prior approaches based on local linearization [4], the CCM approach leads to a convex optimization problem for controller synthesis and generates controllers that stabilize every feasible trajectory in a region, instead of just a single target trajectory that must be known a priori [3].
On the other hand, control design methods to deal with dynamic uncertainties in the deterministic setting can be roughly classified into adaptive and robust approaches. Robust approaches, such as H ∞ control [5], µ synthesis [6] and robust/tube model predictive control (MPC) [7,8], usually consider parametric uncertainties or bounded disturbances and aim to find controllers with performance guarantees for the worst case of such uncertainties. The consideration of worst-case scenarios associated with robust approaches often leads to conservative nominal performance. Disturbance-observer (DOB) based control and related methods such as active disturbance rejection control (ADRC) [9] lump all uncertainties that may include parametric uncertainties, unmodeled dynamics and external disturbances, together as a "total disturbance", estimate it via an observer and then compute control actions to compensate for the estimated disturbance [10] to recover the nominal performance. However, for state-dependent uncertainties, DOB-based control methods usually ignore the dependence of "disturbance" on system states and rely on assumptions on the derivative of the "disturbance" that are difficult to verify for theoretical guarantees [10,11]. Alternatively, adaptive control methods such as model reference adaptive control (MRAC) [12] usually need a parametric structure for the uncertainties, rely on online estimation of the parameters for control law construction and provide asymptotic performance guarantees in most cases. One of the exceptions is L 1 adaptive control [13] that does not need a parameterization of the uncertainties (similar to DOB-based control) and focuses on transient performance guarantees in terms of uniformly bounded error between the ideal and uncertain systems.
Both robust and adaptive control approaches have been explored in the context of CCM-based control in the presence of uncertainties and disturbances. In particular, adaptive control was combined with CCM to handle nonlinear control-affine systems with both parametric [14] and non-parametric uncertainties [15]. The case of bounded disturbances in CCM-based control was addressed by leveraging input-to-state stability analysis [16] or robust CCM [17,18]. CCM for stochastic systems was developed in [19] to minimize the mean squared tracking error in the presence of stochastic disturbances. Closely relevant to this paper, [15] designed an L 1 adaptive controller to augment a baseline CCM-based controller to compensate for matched nonlinear non-parametric uncertainties that can depend on both time and states. The authors of [15] proved that transient tracking performance was guaranteed in the sense that the actual state trajectory exponentially converges to a neighborhood or a tube around the desired one. Compared to [15], our approach relies on a disturbance observer that yields an estimation error bound and robust Riemannian energy condition and ensures that the actual state trajectory exponentially converges to the nominal one.
Statement of Contributions: We present a tracking controller for nonlinear systems subject to matched uncertainties that can depend on both time and states based on contraction metrics and disturbance estimation. Our controller leverages a disturbance estimator to estimate the pointwise value of the uncertainties, with a pre-computable estimation error bound. The estimated disturbance and the estimation error bound are then incorporated into a robust Riemannian energy condition to compute the control law that guarantees exponential convergence of actual state trajectories to nominal ones. We validate the efficacy of our controller on two simulation examples and demonstrate its advantages over existing controllers.
The idea presented in this paper is leveraged in [20] for safe learning of uncertain dynamics using deep neural networks. Compared to [20], this paper is not relevant to learning and allows the uncertainty to be dependent on both time and states, as opposed to the dependence on states only in [20]. Additionally, this paper includes an additional aircraft example for performance illustration and conducts extensive comparisons with existing adaptive approaches in simulations that are not available in [20].
Notations: Let R n , R + and R m×n denote the n-dimensional real vector space, the set of non-negative real numbers and the set of real m by n matrices, respectively. I and 0 denote an identity matrix, and a zero matrix of compatible dimensions, respectively; · denotes the 2-norm of a vector or a matrix. For a vector y, y i denotes its ith element. For a matrix-valued function M : R n → R n×n and a vector y ∈ R n , ∂ y M(x) ∑ n i=1 ∂M(x) ∂x i y i denotes the directional derivative of M(x) along y. For symmetric matrices P and Q, P > Q (P ≥ Q) means P − Q is positive definite (semidefinite). X is the shorthand notation of X + X . Finally, denotes the Minkowski set difference.

Problem Statement and Preliminaries
Consider a nonlinear control-affine system with uncertaintieṡ where x(t) ∈ X ⊂ R n is the state vector, u(t) ∈ U ⊂ R m is the control input vector, f : R n → R n and B : R n → R m are known and locally Lipschitz continuous functions, and d(t, x) represents the matched model uncertainty that can depend on both time and states. We assume that B(x) has full column rank for any x ∈ X . Suppose X is a compact set that contains the origin, and the control constraint set U is defined as U {u ∈ R m : u ≤ u ≤ū}, where u,ū ∈ R m denote the lower and upper bounds of all control channels, respectively. Furthermore, we make the following assumptions on B(x) and d(t, x).

Assumption 1.
There exist known positive constants L B , L d , l d and b d such that for any x, y ∈ X and t, τ ≥ 0, the following inequalities hold: Remark 1. Assumption 1 indicates that the uncertain function d(t, x) is locally Lipschitz in both t and x with known Lipschitz constants and is uniformly bounded by a known constant in the compact set X .
In fact, given the local Lipschitz constants L d and l d , a uniform bound on d(t, x) in X can always be derived by using Lipschitz continuity properties if the bound on d(t, x * ) for an arbitrary x * in X and any t ≥ 0 is known. For instance, assuming for any x ∈ X and g ≥ 0. In practice, some prior knowledge about the actual system and the uncertainty may be leveraged to obtain a tighter bound than the one based on the Lipschitz continuity explained earlier, which is why we directly make an assumption on the uniform bound. With Assumption 1, we will show (in Section 3.3) that the pointwise value of d(t, x(t)) at any time t can be estimated with a pre-computable estimation error bound.
For the system in (1), assume we have a nominal state and input trajectory, x (·) and u (·), which satisfy the nominal, i.e., uncertainty-free, dynamics: We would like to design a state-feedback controller in the form of so that the actual state trajectory x(·) exponentially converges to the nominal one x (·). Our solution is based on CCM and disturbance estimation. Next, we briefly review CCM for uncertainty-free systems.

Control Contraction Metrics (CCMs)
We first introduce some notations related to Riemannian geometry, most of which are from [1]. A Riemannian metric on R n is a symmetric positive-definite matrix function M(x), smooth in x, which defines a "local Euclidean" structure for any two tangent vectors δ 1 and δ 2 through the inner product δ 1 , δ 2 x δ 1 M(x)δ 2 and the norm δ 1 , δ 2 x . A metric is called uniformly bounded if a 1 I ≤ M(x) ≤ a 2 I holds ∀x and for some scalars a 2 ≥ a 1 > 0. Let Γ(a, b) be the set of smooth paths connecting two points a and b in R n , where each c ∈ Γ(a, b) is a piecewise smooth mapping, c : [0, 1] → R n , satisfying c(0) = a, c(1) = b. We use the notation c(s), s ∈ [0, 1], and c s (s) ∂c ∂s . Given a metric M(x) and a curve c(s), we define the Riemannian energy of c(s) as E(c) 1 0 c s M(c(s))c s (s)ds. The Riemannian energy between a and b is defined as E(a, b) inf c∈Γ(a,b) E(c).
Contraction theory [2] draws conclusions on the convergence between pairs of state trajectories toward each other by studying the evolution of the distance between any two infinitesimally close neighbouring trajectories. CCM generalizes contraction analysis to the controlled dynamics setting in which the analysis jointly searches for a controller and a metric that describes the contraction properties of the resulting closed-loop system. Following [1,14], we now briefly review CCMs by considering the nominal, i.e., uncertaintyfree, system:ẋ where x(t) ∈ R n and u(t) ∈ R m . The differential form of (7) is given byδ . Consider a function V(x, δ x ) = δ x M(x)δ x for some positive definite metric M(x), which can be viewed as the Riemannian squared differential length at point x. Differentiating and imposing that the squared length decreases exponentially with rate 2λ, one obtainṡ We first recall some basic results related to CCM. (7) is said to be universally exponentially stabilizable if, for any feasible desired trajectory x (t) and u (t), a feedback controller can be constructed that for any initial condition x(0), a unique solution to (7) exists and satisfies

Definition 1 ([1]). The system
where λ and R are the convergence rate and overshoot, respectively, independent of the initial conditions.

Lemma 1 ([1]
). If there exists a uniformly bounded metric M(x), i.e., α 1 I ≤ M(x) ≤ α 2 I for some positive constants α 1 and α 2 , such that for all x and δ x = 0 satisfying δ x MB = 0, hold, then the system (7) is universally exponentially stabilizable in the sense of Definition 1 via continuous feedback defined almost everywhere, and everywhere in the neighborhood of the target trajectory with the convergence rate λ and overshoot R = α 2 The condition (9) ensures that the dynamics orthogonal to the input are contracting, i.e., (8) holds in the presence of δ x MB = 0 and is often termed as the strong CCM condition [1]. In particular, the condition (9b) can be satisfied by enforcing that each column of B(x) forms a killing vector field for the metric M(x), i.e., M ∂b i ∂x +∂ b i M = 0 for all i = 1, . . . , m. The CCM condition (9) can be transformed into a convex constructive condition for the metric M(x) by a change of variables. Let W(x) = M −1 (x) (commonly referred to as the dual metric), and B ⊥ (x) be a matrix whose columns span the null space of the input matrix B (i.e., B ⊥ B = 0). Then, condition (9) can be cast as convex constructive conditions for W(x): The existence of a contraction metric M(x) is sufficient for stabilizability via Lemma 1. What remains is constructing a feedback controller that achieves the universal exponential stabilizability (UES). As mentioned in [1,16], one way to derive the controller is to interprete the Riemann energy, E(x (t), x(t)), as an incremental control Lyapunov function and use it to construct a min-norm controller that renders for any time ṫ Specifically, at any time t > 0, given the metric M(x) and a desired/actual state pair (x (t), x(t)), a minimum-energy path, i.e., a geodesic, γ(·, t) connecting these two states (i.e., γ(0, t) = x (t) and γ(1, t) = x(t)), can be computed (e.g., using the pseudospectral method in [21] to solve a nonlinear programming problem). Consequently, the Riemannian energy of the geodesic is defined as ∂γ ∂s , can be calculated. As noted in [16], from the formula for the first variation of energy [22] Therefore, (11) can be rewritten as . Therefore, the control signal with a minimum norm for u(t) − u (t) can then be obtained by solving the following quadratic programming (QP) problem: at each time t, which is guaranteed to be feasible under condition (9) [1]. The minimization problem (13) is often termed as the pointwise min-norm control problem and has an analytic solution [23]. The above discussions can be summarized in the following theorem. The proof is trivial by following Lemma 1 and the subsequent discussions and is thus omitted.
. Given a nominal system (7), assume that there exists a uniformly bounded metric W(x) that satisfies (10) for all x ∈ R n . Then, the control law constructed by solving (13) with M(x) = W −1 (x), universally exponentially stabilizes the system (7) in the sense of Definition 1, where R = α 2 α 1 with α 1 and α 2 being two positive constants satisfying

Remark 2.
According to Definition 1 and Theorem 1, under the conditions of Theorem 1, given any feasible trajectory (x ((·), u (·)) of (7), a controller can always be constructed to ensure that the actual state trajectory x(·) exponentially converges to x (·).

Robust Trajectory Tracking Using CCM and Disturbance Estimation
In Section 2, we have shown that existence of a CCM for a nominal (i.e., uncertaintyfree) system can be used to construct a feedback control law to guarantee the universal exponential stabilizability (UES) of the system. In this section, we present a controller based on CCM and disturbance estimation to ensure the UES of the uncertain system (1), whose architecture is depicted in Figure 1.

Uncertain System QP with Robust Riemannian Energy Condition
Disturbance Estimator

CCMs for the Actual System
To apply the contraction method to design a controller to guarantee the UES of the uncertain system (1), we need to first search a valid CCM for it. Following Section 2, we can derive the counterparts of the strong CCM condition (9) or (10). Due to the particular structure with (1) attributed to the matched uncertainty assumption, we have the following lemma. A similar observation has been made in [14] for the case of matched parametric uncertainties. The proof is straightforward and thus omitted. One can refer to [14] for more details.

Lemma 2.
The strong (dual) CCM condition for the uncertain system (1) is the same as the strong (dual) CCM condition, i.e., (9) and (10), for the nominal system.

Remark 3.
As a result of Lemma 2, a metric M(x) (dual metric W(x)) satisfying the condition (9) and (10) for the nominal system (7) is always a CCM (dual CCM) for the true system (1).
x) ∈ D for any t ≥ 0 and x ∈ X . As mentioned in Section 2, given a CCM and a desired trajectory x (t) and u (t) for a nominal system, a control law can be constructed to ensure exponential convergence of the actual state trajectory x(t) to the desired state trajectory x (t). In practice, we have access to only the nominal dynamics (5) instead of the true dynamics to plan a trajectory x (t) and u (t). The following lemma gives the condition when x (t), planned using the nominal dynamics (5), is also a feasible state trajectory for the true system. Lemma 3. Given a desired trajectory x (t) and u (t) satisfying the nominal dynamics (5) with then x (t) is also a feasible state trajectory for the true system (1).
x (t)) ∈ D, which is due to x (t) ∈ X and Assumption 1, we haveū (t) ∈ U . By comparing the dynamics in (1) and (5), we conclude that x (t) andū (t) satisfy the true dynamics (1) and thus are a feasible state and input trajectory for the true system.
Lemma 3 provides a way to verify whether a trajectory planned using the nominal dynamics is a feasible trajectory for the true system in the presence of actuator limits. In the absence of such limits, any feasible trajectory for the learned dynamics is also a feasible trajectory for the true dynamics due to the particular structure of (1) associated with the matched uncertainty assumption.

Robust Riemannian Energy Condition
Section 2 shows that, given a nominal system and a CCM for such a system, a control law can be constructed via solving a QP problem (13) with a condition to constrain the decreasing rate of the Riemannian energy, i.e., condition (12). When considering the uncertain dynamics in (1), the condition (12) becomes (5). Several observations follow immediately. First, it is clear that (15) is not implementable due to its dependence on the true uncertainty d(x(t)) throughẋ(t). Second, if we could have access to the pointwise value of d(x(t)) at each time t, (15) will become implementable even when we do not know the exact functional representation of d(x). Third, if we could estimate the pointwise value of d(x(t)) at each time t with a bound to quantify the estimation error, then we could derive a robust condition for (15). Specifically, assume d(x(t)) is estimated asd(t) at each time t with a uniform estimation error bound (EEB) δ, i.e., d (t) − d(x(t)) ≤ δ, ∀t ≥ 0. Then, we could immediately get the following sufficient condition for (15): Moreover, since M(x) satisfies the CCM condition (9), u(t) that satisfies (16) is guaranteed to exist for any t ≥ 0, regardless of the size of δ, if the input constraint set U is sufficiently large. We term condition (16) the robust Riemannian energy (RRE) condition.

Disturbance Estimation with a Computable EEB
We now introduce a disturbance estimation scheme to estimate the pointwise value of the uncertainty d(x) with a pre-computable EEB, which can be systematically improved by tuning a parameter in the estimation law. The estimation scheme is based on the piecewiseconstant estimation (PWCE) law in [24], which was originally from [25]. The PWCE law consists of two elements, namely a state predictor and a piecewise-constant update law. The state predictor is defined as: is the prediction error, and a is an arbitrary positive constant. The estimation,σ(t), is updated in a piecewise-constant way: where T is the estimation sampling time, and i = 0, 1, 2, · · · . Finally, the pointwise value of d(x(t)) at time t is estimated asd where B † (x(t)) is the pseudoinverse of B(x(t)). The following lemma establishes the EEB associated with the estimation scheme in (18) and (19). The proof is similar to that in [24]. For completeness, it is given in Appendix A.

Lemma 4.
Given the dynamics (1) subject to Assumption 1, and the estimation law in (18) and (19), if x ∈ X and u ∈ U for any t ≥ 0, the estimation error can be bounded as where φ max x∈X ,u∈U with constants L B , L d and b d from Assumption 1, and φ defined in (23). Moreover, lim T→0 δ(t, T) = 0, for any t ≥ T.

Proof. See Appendix A.
Remark 4. Lemma 4 implies that theoretically, for t ≥ T, the disturbance estimation after a single sampling interval can be made arbitrarily accurate by reducing T, which further indicates that the conservatism with the RRE condition can be arbitrarily reduced after a sampling interval.
In practice, the value of T is subject to the limitations related to computational hardware and sensor noise. Additionally, using a very small T tends to introduce high frequency components in the control loop, potentially harming the robustness of the closed-loop system, e.g., against time delay. This is similar to the use of a high adaptation rate in model reference adaptive control schemes as discussed in [13]. Therefore, one should avoid the use of a very small T for the sake of robustness unless a low-pass filter is used to filter the estimated disturbance before fed into (16), as suggested by the L 1 adaptive control theory [13].

Remark 5.
The estimation in [0, T) cannot be arbitrarily accurate. This is because the estimation in [0, T) depends onx(0) according to (19). Considering thatx(0) is purely determined by the initial state of the system, x(0), and the initial state of the predictor,x(0), it does not contain any information of the uncertainty. Since T is usually very small in practice, lack of a tight estimation error bound for the interval [0, T) will not cause an issue from a practical point of view. Additionally, the estimation of φ defined in (23) could be quite conservative. Further considering the frequent use of Lipschitz continuity and inequalities related to matrix/vector norms in deriving the constant α(T), α(T) can be overly conservative. Therefore, for practical implementation, one should leverage some empirical study, e.g., performing simulations under a few user-selected functions of d(t, x) and determining a bound for δ(t, T). In our experiments, we found the theoretical bound δ(t, T) computed according to (21) was usually at least 10 and could be 10 4 times more conservative.

Exponentially Convergent Trajectory Tracking
Based on the review of contraction control in Section 2 and the discussions in Sections 3.2 and 3.3, the control law can be obtained by solving the following QP problem at each time t: whereẋ(t) = f (x) + B(x)(k +d(t)), according to (17), depends ond(t), which is from the disturbance estimation law defined by (18) to (20), δ(t, T) as defined in (21), anḋ x (t) = f (x ) + B(x )u as defined in (5). Similar to (13), problem (24) is a pointwise min-norm control problem and has an analytic solution [23]. Specifically, denot- (25) can be written as φ 0 (t, x , x) + φ 1 (x , x)(k − u (t)) ≤ 0, and the solution for (24) is given by To move forward with analysis, we need to verify that when x(t), x (t) ∈ X , the control signal u(t) resulting from solving the QP problem (24) satisfies u(t) ∈ U . Deriving verifiable conditions to ensure this set bound is outside the scope of this paper and will be addressed as future work. We are now ready to state the main result of the paper. Theorem 2. Given an uncertain system represented by (1) satisfying Assumption 1, assume that there exists a metric W(x) such that for all x ∈ X , (10) holds and α 1 I ≤ M(x) = W −1 (x) ≤ α 2 I holds for positive constants α 1 and α 2 . Furthermore, suppose that a nominal trajectory (x (·), u (·)) planned using the nominal dynamics (5) and the initial actual states x(0) satisfy (14) and for any t ≥ 0. Then, if u(t) from solving (24) satisfies u(t) ∈ U for any t ≥ 0, the control law constructed by solving (24) ensures x(t) ∈ X for any t ≥ 0, and furthermore, universally exponentially stabilizes the uncertain system (1) in the sense of Definition 1 with R = α 2 α 1 , i.e., Proof. We use contradiction to show x(t) ∈ X for all t ≥ 0. Assume this is not true. According to (27), x(0) ∈ X . Since x(t) is continuous, there must exist a time τ such that Now let us consider the system evolution in [0, τ − ]. Since u(t) ∈ U by assumption and x(t) ∈ X for any t in [0, τ − ], the EEB in (21) holds in [0, τ − ]. As a result, the control law obtained from solving (24) ensures satisfaction of the RRE condition (16) and thus satisfaction of the Riemannian energy condition (15) for the uncertain system (1), and thereby universally exponentially stabilizes the uncertain system (1) in [0, τ − ], in the sense of Definition 1 with R = α 2 α 1 , according to Theorem 1. On the other hand, satisfaction of (14) implies that x (t) is a feasible state trajectory for the uncertain system (1) according to Lemma 3. Further considering Theorem 1, we have (27), the preceding inequality indicates that x(t) remains in the interior of X for t in [0, τ − ]. This, together with the continuity of x(t), immediately implies x(τ) ∈ X , which contradicts (29). Therefore, we conclude that x(t) ∈ X for all t ≥ 0. From the development of the proof, it is clear that with the control law given by the solution of (24), the UES of the closed-loop system in the sense of Definition 1 with R = α 2 α 1 for all t ≥ 0 is achieved, which is mathematically represented by (28). The proof is complete.

Discussion
Theorem 2 essentially states that under certain assumptions, the proposed controller guarantees exponential convergence of the actual state trajectory x(t) to a desired one x (t).
With the exponential guarantee, if the actual trajectory meets the desired trajectory at certain time τ, then these two trajectories will stay together afterward. While the exponential convergence guarantee is stronger than the performance guarantees provided by existing adaptive CCM-based approaches [14,15] that deal with similar settings (i.e., matched uncertainties), the proposed method requires the knowledge of the Lipschitz bound of the uncertainty d(t, x) and the input matrix function B(x) to be in a compact set known a priori (see Assumption 1), and the actual control inputs to stay in a compact set known a priori, which cannot be verified at this moment due to the lack of a bound on the control inputs. These requirements are not needed in [14,15].
The approach here is related to the robust control Lyapunov-based approaches [23] which provide robust stabilization around an equilibrium point (as opposed to a trajectory considered in this paper) in the presence of uncertainties.

Remark 6.
The exponential convergence guarantee stated in Theorem 2 is based on a continuoustime implementation of the controller. In practice, a controller is normally implemented on a digital processor or controller with a fixed sampling time. As a result, the property of exponential convergence may be slightly violated.
Computational cost: As can be seen from Sections Section 2 and 3.2-3.4, computation of the control signal at each time t includes three steps: (i) updating the estimated disturbancê d(t) via (18) to (20), (ii) computing the geodesic γ(·, t) connecting the actual and nominal states (see the discussion below (11)), and (iii) computing the control signal u(t) via (26). The computation costs of steps (i) and (iii) are quite low as they only involve integration and algebraic calculation. In comparison, step (ii) has a relatively high computational cost as it necessitates solving a nonlinear programming (NLP) problem. However, since the NLP problem does not involve dynamic constraints, it is much easier to solve than a nonlinear model predictive control (MPC) problem [21]. Following [21], such a problem can be efficiently solved by applying a pseudospectral method.

Simulation Results
In this section, we illustrate the performance of our proposed tracking controller based on the RRE condition and disturbance estimation, denoted as DE-CCM, using aircraft and planar quadrotor examples. For both examples, we perform comparisons of DE-CCM with standard CCM controllers that ignore the uncertainties and adaptive CCM (Ad-CCM) controllers considering parametric uncertainties designed using the approach in [14]. All the computations and simulations were performed in Matlab R2021b.

Longitudinal Dynamics of an Aircraft
We first implement our method on the simplified pitch dynamics of an aircraft borrowed from [26] where θ, α and q are the pitch angle (in rad), angle of attack (in rad), and pitch rate (in rad/s). HereL(α) andM(α) are aerodynamic lift and moment, respectively. Using the flat plat theory, these two aerodynamic terms are approximated byL(α) = 0. We set the convergence rate λ to 1. By gridding the set of α and evaluating the constraints (9) in those grid points, we found a CCM W(x) as a quadratic function of α with the SPOT toolbox [27] (to formulate the convex optimization problem) and Mosek solver [28].
Additionally, the constants α 1 and α 2 in (28) such that α 1 I ≤ M(x) = W −1 (x) ≤ α 2 I for all x ∈ X were found to be α 1 = 0.1 and α 2 = 396.5. We planned a nominal trajectory (x (·), u (·)) using OptimTraj [29,30], to drive the system from the initial states [0, 0, 0] T to the terminal states [180 • , 0, 0] T , while minimizing the task completion time (T a ) and energy consumption characterized by the cost function J = T a 0 u(t) 2 dt + 5T a . For simulation, OPTI [31] and Matlab fmincon solvers were used to solve the geodesic optimization problem (see Section 2). The initial states of the actual system were chosen to be [5 • , 5 • , 0] T , slightly deviated from that planned ones to better illustrate the tracking performance. We implemented our proposed DE-CCM, Ad-CCM from [14] and a standard CCM which neglects all the uncertainty. For Ad-CCM design, the adaptation gain was chosen to be diag([10 3 , 10 3 ]) to achieve a relatively good tracking result, while further increasing it did not help much with the tracking performance. The design procedure for Ad-CCM in [14] requires a parametric structure for the uncertainty, which is given by where ∆(t, x) is the known base function and θ is the unknown parameter vector to be estimated by the adaptive law proposed in [14]. The control signals under all three controllers were updated at 200 Hz. It is easy to notice from Assumption 1 that L B = 0 and l d = 0 since the input matrix is constant, and the uncertainty is time-invariant. We can also verify that the disturbance is bounded by b d = 2.12 and has a Lipschitz constant L d = 3.80. By gridding the space X and making use of the control input bound, the system derivative can also be bounded by a constant φ = 2.11. According to (21), if we want to achieve an EEB δ(t, T) = 0.05 for all t ≥ T, the maximum value for the estimation sample time T is 7.76 × 10 −4 s. However, as mentioned in Remark 5, the way to compute the EEB is quite conservative. In the simulations we found that T s = 0.005 was more than enough to ensure the EEB and therefore used T s = 0.005 for implementing the DE-CCM controller.
As shown in Figures 2 and 3, due to ignoring the uncertainties, CCM yielded a large tracking error between 2 and 6 s. The state trajectories under Ad-CCM had some oscillations, which lasted roughly up to 8 s. All three states yielded by DE-CCM achieve good tracking performance without large deviations from the planned trajectories, unlike the performance yielded by Ad-CCM and CCM. From Figure 3 we notice that the tracking error represented by x − x under DE-CCM monotonically decreases and achieves the smallest steady-state error. The small non-zero tracking error at the end under DE-CCM, which is inconsistent with the performance guarantee in (25), is due to the limited control update frequency, while the performance guarantee in Lemma 4 holds under continuous update of the control signal, i.e., corresponding to an infinitely high update frequency. Table 1 shows the mean squared error (MSE) for state trajectory tracking defined by where N is the number of data points, under DE-CCM, Ad-CCM and CCM. We observe that DE-CCM outperforms CCM and Ad-CCM in terms of MSE by 54% and 2%, respectively.    From Figure 4, we observe that the input of DE-CCM is smoother than Ad-CCM. The small oscillations between 2 s to 8 s in DE-CCM input are due to the finite tolerance in the geodesic optimization. Decreasing the tolerance and the sample time will reduce the oscillations but request more iterations (and thus more time) to compute the control signal at each time step.

Planar Quadrotor
A planar quadrotor system is borrowed from [16]. The state vector is defined as where p x and p z are the positions in x and z directions, respectively, v x and v z are the slip velocity (lateral) and the velocity along the thrust axis in the body frame of the vehicle, φ is the angle between the x direction of the body frame and the x direction of the inertia frame. The input vector u = [u 1 , u 2 ] contains the thrust force produced by each of the two propellers. The dynamics of the vehicle are given bẏ where m and J denote the mass and moment of inertia about the out-of-plane axis, and l is the distance between each of the propellers and the vehicle center, and d(t, x) denotes the unknown disturbances exerted on the propellers. The parameters were set as m = 0.486 kg, J = 0.00383 Kg m 2 , and l = 0.25 m. The uncertainty d(t, x) was set to be d(t, x) = 0.15(v 2 x + v 2 z )[−1 + 0.3 sin(2t); −1 + 0.3 cos(2t)]. We imposed the following constraints: x ∈ X [0, 10] When searching for CCM, we parameterized the CCM W by φ and v x and imposed the constraint W ≥ 0.01I. The convergence rate λ was chosen to be 0.8. More details about synthesizing the CCM can be found in [17]. For estimating the disturbance using (18) to (20), we set a = 10. It is easy to verify that L d = 1.23, l d = 0.64, b d = 1.38, and L B = 0 (due to the fact that B is constant) satisfy (2). By gridding the space X , the constant φ in (23) can be determined as φ = 584.6. According to (21), if we want to achieve an EEB δ(t, T) = 0.1 for all t ≥ T, then the estimation sampling time needs to satisfy T ≤ 6.42 × 10 −7 s. However, as noted in Remark 5, the EEB computed according to (21) could be overly conservative.
In the simulations, we found that the estimation sampling time of 0.002 s was more than enough to ensure the desired EEB and therefore simply set T = 0.002 s.
We consider the task of navigation from (2, 0) to (8,8) while avoiding three obstacles depicted by black circles in Figure 5. A nominal trajectory (x (·), u (·)) was generated using OptimTraj [29] to minimize the cost J = T a 0 u(t) 2 dt + 5T a , where T a is the arrival time. OPTI [31] and Matlab fmincon solvers were used to solve the geodesic optimization problem (see Section 2). The actual start point was set to be (0, 0), which was different from the planned start point, to reveal the trajectory convergence pattern. For comparison, we also designed a standard CCM controller by completely ignoring the uncertainty and two adaptive CCM (Ad-CCM) controllers following the approach in [14]. To apply the approach in [14] which can only handle parametric uncertainties, we parameterized the uncertainty as where ∆(t, x) is the basis function that is assumed to be known, and θ is th vector of unknown parameters. With the parametric structure (33), we designed two adaptive CCM controllers using Γ = 10 and Γ = 100, respectively, where Γ denotes the adaptive gain. Figure 5 shows the planned and actual trajectories under the CCM, Ad-CCM, and our proposed controller based on the RRE condition and disturbance estimation, denoted as DE-CCM, while Figures 6 and 7 show the control inputs and Riemannian energy. One can see that the actual trajectories yielded by the CCM controller deviated quite a lot from the planned ones and collided with one obstacle. On the other hand, the actual trajectories yielded by the DE-CCM controller converged to the desired trajectory as expected and almost overlapped with it afterward. In fact, the slight deviations of actual trajectories from the desired ones under the DE-CCM controller were due to the finite step size associated with the ODE solver used for the simulations (see Remark 6). Table 2 shows the MSE for state trajectory tracking defined in (32). We observe that DE-CCM outperforms CCM and Ad-CCM in terms of MSE by 46% and 14%, respectively.   From Figure 6, one can see that the magnitude of E(x , x) under the RD-CCM controller decreased exponentially, and the magnitude was bounded by the curve E(x (0), x(0))e −2λt from above except at the very end when the energy is close to zero. In comparison, Ad-CCM with Γ = 100 yielded similar tracking performance to DE-CCM, while the tracking performance of Ad-CCM with Γ = 10 was relatively worse and characterized by larger oscillations. Additionally, from Figure 7, one can see that the control inputs generated by both of the Ad-CCM controllers have high-frequency oscillations before 3 s, which is undesired for practical deployment. Finally, Figure 8 shows the actual and estimated disturbances as well as the estimation error. One can see that the estimated disturbance is quite close to the actual one for both channels, and the EEB of 0.1 is respected throughout the simulation.

Conclusions
This paper presents a robust trajectory tracking controller with exponential convergence for uncertain nonlinear systems based on control contraction metrics (CCM) and disturbance estimation. The controller uses a disturbance estimator to estimate the pointwise value of the uncertainty with a pre-computable estimation error bound (EEB). The estimated disturbance and the EEB are then incorporated into a robust Riemannian energy condition, which guarantees exponential convergence of actual trajectories to desired ones. The efficacy of the proposed controller is validated in simulations. In particular, the proposed controller outperforms an existing adaptive CCM controller in terms of tracking performance by 2% for the aircraft example and 14% for the planar quadrotor example, while not needing to know the basis functions to parameterize the uncertainties that are needed by the adaptive CCM controller.
This paper considers only matched uncertainties, which are added to the system through the same channels as control inputs. In the future, we would like to address unmatched uncertainties that widely exist in practical systems. Additionally, we would like to experimentally validate the proposed controller on real hardware.

Conflicts of Interest:
The authors declare no conflict of interest.