On-The-Fly Control of Unknown Systems: From Side Information to Performance Guarantees through Reachability

We develop data-driven algorithms for reachability analysis and control of systems with a priori unknown nonlinear dynamics. The resulting algorithms not only are suitable for settings with real-time requirements but also provide provable performance guarantees. To this end, they merge noisy data from only a single finite-horizon trajectory and, if available, various forms of side information. Such side information may include knowledge of the regularity of the dynamics, algebraic constraints on the states, monotonicity, or decoupling in the dynamics between the states. Specifically, we develop two algorithms, $\texttt{DaTaReach}$ and $\texttt{DaTaControl}$, to over-approximate the reachable set and design control signals for the system on the fly. $\texttt{DaTaReach}$ constructs a differential inclusion that contains the unknown dynamics. Then, in a discrete-time setting, it over-approximates the reachable set through interval Taylor-based methods applied to systems with dynamics described as differential inclusions. We provide a bound on the time step size that ensures the correctness and termination of $\texttt{DaTaReach}$. $\texttt{DaTaControl}$ enables convex-optimization-based control using the computed over-approximation and the receding-horizon control framework. Besides, $\texttt{DaTaControl}$ achieves near-optimal control and is suitable for real-time control of such systems. We establish a bound on its suboptimality and the number of primitive operations it requires to compute control values. Then, we theoretically show that $\texttt{DaTaControl}$ achieves tighter suboptimality bounds with an increasing amount of data and richer side information. Finally, experiments on a unicycle, quadrotor, and aircraft systems demonstrate the efficacy of both algorithms over existing approaches.


I. INTRODUCTION
Consider a scenario in which significant and unexpected changes in the dynamics of a system occur. The changes in the dynamics are such that the a priori known model cannot be used, and there is a need to learn the new dynamics on the fly. In such a scenario, the system has access to data from only its current trajectory and needs to retain a certain degree of control. This paper considers the problem of data-driven, onthe-fly control of systems with unknown nonlinear dynamics under severely limited data.
We develop data-driven algorithms, DaTaReach and DaTaControl, for the reachability analysis and control of systems with a priori unknown dynamics. DaTaReach exploits the data from a given and and ongoing trajectory of the system to compute an over-approximation of the set of states it This material is based on work partly supported by Air Force Office of Scientific Research (FA9550-19-1-0005) and National Aeronautics and Space Administration (80NSSC19K0209).
F. Djeumou Over-approximations for near-optimal controls Near-optimal control synthesis for t ≥ tN Optimal trajectory Fig. 1. We use limited data and side information for on-the-fly control of systems with unknown dynamics. At each sampling time, we compute an over-approximation of the set of states the system may reach to describe the uncertainty in its unknown trajectory and construct one-step optimal controllers via convex optimization. may reach. DaTaControl incorporates the computed overapproximation into a constrained optimal control problem, which is then solved on the fly. Specifically, DaTaReach and DaTaControl can work with possibly noisy data from only a single finite-horizon trajectory of the system and take advantage of various forms of side information on the dynamics. More specifically, the data include finite (noisy) samples of the states, the derivatives of the states, and the control signals applied. The side information may be a priori knowledge of the regularity of the dynamics, bounds on the vector field locally, monotonicity of the vector field, decoupling in the dynamics among the states, algebraic constraints on the states, or knowledge of parts of the dynamics. For example, such side information may be due to known elementary laws of physics or may be directly extracted from the data, as we illustrate in the numerical simulations.
DaTaReach provides closed-form expressions for overapproximations of the reachable set of unknown dynamical systems. It first uses the available data and the given side information to construct differential inclusions that contain the unknown vector field. Unlike existing work [1] that provides state-independent differential inclusions, the constructed differential inclusions account for dependencies on the states and control signals. Then, DaTaReach provides closedform expressions for over-approximations of the reachable set associated with these data-driven differential inclusions. Specifically, it builds on interval Taylor-based methods [2]- [4] to over-approximate the reachable set of dynamics described by the differential inclusions. The closed-form expressions can incorporate the available side information, and increasing amount of data and richer side information provide tighter over-approximations of the reachable set.
The closed-form expressions enable convex-optimizationbased, on-the-fly control of unknown dynamical systems through the one-step receding-horizon control framework [5]. Specifically, we seek to sequentially minimize a given one-step cost function in a discrete-time setting. We refer to such an optimization problem as a one-step optimal control problem. The one-step cost function, which encodes the desired behavior of the system, has to be optimized in a black-box manner since it is typically a function of the unknown future state, in addition to the known current state and the currently applied control input. DaTaControl computes approximate solutions to the one-step optimal control problem through convex optimization relaxations. The convex relaxations are obtained by replacing the unknown future state with a control-affine linearization of the corresponding over-approximation of the reachable set.
More specifically, we prove that DaTaControl achieves near-real-time and near-optimal control of the unknown system. Indeed, we provide bounds on the suboptimality of the relaxed convex problems' solutions with respect to the optimal solutions if the dynamics were known. Then, we theoretically show that the obtained bounds on the suboptimality become tighter with an increasing amount of data and richer side information. Besides, we prove that DaTaControl is suitable for real-time control of the system through an explicit upper bound on the number of primitive operations required to terminate.
Empirically, through a series of experiments on a unicycle, quadrotor, and aircraft systems, we show that DaTaControl is significantly more computationally efficient, robust to noise, and less suboptimal than the baselines SINDYc [6], DeePC [7], C2Opt [8], and CGP-LCB [9]. Indeed, the experiments show that DaTaControl is three to five orders of magnitude faster than the baselines, achieves the control objective using at least three times fewer system's interactions (less suboptimal), and outperforms the baselines in terms of robustness to noise. Contributions. This paper significantly extends our conference paper [10] by first augmenting DaTaControl and DaTaReach with the ability to handle noise in the data. Indeed, the approach in [10] assumes exact measurements of the state and the state derivative. Then, we provide a bound on the discrete-time step size such that DaTaReach and DaTaControl are guaranteed to be correct and to terminate. The algorithms in [10] assume "sufficiently small" time step to work, without any explicit bound to characterize such small time step. Finally, we prove an explicit bound on the number of primitive operations required by DaTaControl to terminate and provide additional numerical experiments. With respect to state-of-the-art data-driven techniques, this paper proposes near-optimal and real-time data-driven algorithms that work with an extremely scarce amount of data and that can incorporate side information on the underlying dynamics when learning to control nonlinear systems. We provide performance guarantees such as bounds on the suboptimality and computational complexity of the algorithms. At last, through extensive comparisons with state-of-the-art baselines, we empirically demonstrate the efficiency of the proposed approach. Related work. Several approaches for data-driven control combine system identification with model predictive control [6], [8], [11]- [13]. In [11], the authors use Koopman theory to lift the unknown nonlinear dynamics to a higherdimensional space where they perform linear system identifi-cation. SINDYc [6] utilizes a sparse regression over a library of nonlinear functions for nonlinear system identification. DMDc [12] uses the spectral properties of the collected data to obtain approximate linear models. Myopic control [13] uses a finite sequence of perturbations to learn a local linear model, which is then used to optimize a model-driven goodness function that encodes desirable behaviors. DMDc, SINDYc, Myopic control, and approaches based on the Koopman theory require significantly more data than the proposed approach, do not provide real-time guarantees, and, in their current forms, cannot incorporate any of the side information in this paper.
Contextual optimization-based approaches tackle the datadriven control problem via surrogate optimization and skip the system identification step [8], [9]. These approaches iteratively minimize the one-step cost function in a black-box manner using the data. C2Opt [8] exploits the structure in the given problem and utilizes side information (for example, smoothness) and convex optimization to solve this problem. It overcomes the drawbacks of the traditional Gaussian processbased approaches [9], namely high computational costs, expensive hyperparameter tuning, and inability to incorporate side information. Unlike DaTaControl, C2Opt considers limited forms of side information and relies on the knowledge of the gradient of the one-step cost, that may not be accessible.
Recent work [7], [14]- [18] have proposed data-driven control techniques for unknown dynamical systems that bypass the system identification step. These techniques mostly assume linear time-invariant dynamical systems and are extremely performant in such a setting. In contrast, the approach in this paper works with general unknown nonlinear dynamical systems and can incorporate side information on the underlying dynamics. Besides, on a linear version of the aircraft dynamics, we provide a comparison with DeePC [7], which is a traditional baseline using the same behavioral systems theory foundation [19] as the approaches in [14]- [18]. We show that DaTaControl is two orders of magnitude faster than DeePC while being more accurate with a low amount of data and as accurate with a high amount of data.
The authors of [20]- [22] have considered the problem of data-driven estimation of the reachable sets of partially unknown dynamical systems. The approaches in these work rely on system identification using either supervised learning algorithms [22] or Gaussian process-based algorithms [20], [21]. Such approaches are unable to take advantage of the side information, require significantly more data than DaTaReach, and provide only probabilistic guarantees of the correctness of the computed reachable sets. Instead, DaTaReach provides correct over-approximations of the reachable set at the expense of sometimes being conservative.

II. PRELIMINARIES
We denote an interval by [a, b] = {x ∈ R|a ≤ x ≤ b} for some a, b ∈ R such that a ≤ b, the set {i, . . . , j} by N [i,j] for i, j ∈ N with i ≤ j , the 2-norm by || · || 2 , the k th component of a vector x and the (k, j) component of a matrix X by (x) k and (X) k,j , respectively, and the Lipschitz constant of , referred as f is C k , with k ≥ 0 if f is continuous on X ⊆ R n and all the partial derivatives of order 1, . . . , k exist and are continuous on X , and f is piecewise-C k with k ≥ 0 if there exists a partition of X such that f is C k on each set in the partition. Comparisons (e.g., ≥) between matrices or vectors are conducted elementwise.

A. Interval Analysis
We denote the set of intervals on R by IR = {A = [A, A] | A, A ∈ R, A ≤ A}, the set of n-dimensional interval vectors by IR n , and the set of n × m-dimensional interval matrices by IR n×m . We carry forward the definitions [23] of arithmetic operations, set inclusion, and intersections of intervals to interval vectors and matrices by applying them componentwise. We denote the absolute value of an interval A by |A| = max{|A|, |A|}, the infinity norm of A ∈ IR n by A ∞ = sup i∈N [1,n] |(A) i |, the Cartesian product between intervals by ⊗, i.e., A⊗B = [[A, A], [B, B]] ∈ IR 2 for any two intervals A, B ∈ IR. We define A n as the Cartesian product of any interval A ∈ IR with itself n times. We use the term interval to specify an interval vector or interval matrix when it is clear from the context.
Thus, given an interval A, f (A) is an interval that overapproximates the range of values taken by f over A. We denote interval-valued vector functions via bold lowercase and interval-valued matrix functions via bold uppercase symbols.

B. Over-Approximations of the Reachable Set
Consider a nonlinear dynamical system, where the state x : R + → X is a continuous-time signal evolving in X ∈ IR n , the control u ∈ U is a signal of time evolving in the control set U ∈ IR m with U = {v : R + → U | v is piecewise-C Du } for D u ≥ 0, and the function h : Definition 1 (TRAJECTORY OF THE SYSTEM). Given an initial state x i = x(t i ) at time t i and a control signal u ∈ U, a trajectory of (4) is a continuous function of time We are interested in the set of states reachable by trajectories of the system when the initial state and the control signal are uncertain quantities. We call such a set the reachable set of the system, and we define it as follows.
Definition 2 (REACHABLE SET). Given a set I i ⊆ X of states at time t i and a set V ⊆ U of control signals, the reachable set of the dynamics (4) at time t ≥ t i is given by Given a set V ⊆ U and a set I 0 ⊆ X of states at time t 0 , we compute over-approximations of R(t, I 0 , V) at time t ≥ t 0 using interval Taylor-based methods [3], [23], [24]. Specifically, we consider a time grid where u (l) is the l th derivative of u and h [d] is a function of higher order derivatives of h. Next, we start with R + 0 = I 0 , and compute the over-approximations Here, the set S i ⊆ X is an a priori rough enclosure of R(·, and is a solution of Since it is a hard problem to compute exactly the range , we instead compute the a priori rough enclosure S i recursively by solving where h is an interval extension of the known function h. Note that a solution S i of (8) is also a solution of (7) by definition of the interval extension h. The existence of the a priori rough enclosure S i is guaranteed for a sufficiently small time step size t i+1 − t i . Hence, picking adequately the time step size is primordial in order to find S i in a finite number of iterations when solving the fixed-point equation (8).

C. One-Step Optimal Control
Constrained receding-horizon control [5] is a form of constrained control in which, at each sampling time, we obtain the current control input by solving a finite-horizon openloop optimal control problem, using the current state of the system as the initial state. Consider a system with dynamics in the form (4), the desired behavior of the system is generally encoded as a one-step cost function c : X × U × X → R that assigns preferences over the tuple of the current state, the current control input, and the corresponding state at the next time step. Given a constant time step size ∆t ≥ 0, we define the next state x p+1 = x(t p+1 ; x p , u p ) as the value of the trajectory of the system at time t p+1 = t p + ∆t, given the initial state x p ∈ X at time t p and the constant control signal u p ∈ U. Therefore, the finite horizon open-loop optimal control problem at a state x i ∈ X at time t i is given by where N ≥ 1 is the planning horizon of the control problem. Note that even if the cost function c is convex, the optimization problem (9) is generally nonconvex since the next state x p+1 is possibly highly nonlinear in the control variable u p due to the nonlinear dynamics function h.
In the setting where the dynamics of the system are unknown, we cannot compute its trajectories, and c(x p , u p , x p+1 ) is therefore unknown. Hence, the optimization problem (9) is similar to the problem of optimizing an a priori unknown function. In such a setting, the inability to accurately predict the next states of the system motivates using the short planning horizon N = 1. Thus, at each time t i , we solve the one-step optimal control problem [8], [25], [26] given by where x i is the state of the system at t i , x i+1 is the unknown next state at time t i+1 = t i + ∆t, and ∆t is the constant time step size. The reduction to a one-step optimal control problem greatly improves the tractability of the problem of controlling the unknown dynamical system [25].

III. PROBLEM STATEMENT
In this paper, we consider control-affine nonlinear systems, where f : X → R n and G : X → R n×m are unknown vectorvalued and matrix-valued functions, respectively, for X ∈ IR n . Note that even though we consider control-affine dynamics, in the general case, we can construct a control-affine model of the system locally and apply the results of the paper. Assumption 1 is common in the frameworks of reachability analysis and receding-horizon control. Since the domain X ∈ IR n is bounded, f and G are globally Lipschitz-continuous on X . We exploit such global Lipschitz continuity on X by assuming known upper bounds on the Lipschitz constants. That is, we have access to L f ∈ R n + and L G ∈ R n×m + with (L f ) k = L f k and (L G ) k,l = L G k,l as known upper bounds on the Lipschitz constants of (f ) k and (G) k,l for all k ∈ N [1,n] and l ∈ N [1,m] . We emphasize that the Lipschitz bounds can be directly estimated from data at the expense of weakening some of the guarantees in this paper. Let denote a single finite-horizon trajectory containing N noisy samples of the exact state x i = x(t i ), the derivativeẋ i =ẋ(t i ) of the state, and the control signal u i = u(t i ) from a trajectory of the system. We build on the widely-used bounded noise assumption and consider that the following bounds hold: for all t ∈ R + and for some vector values η,η ∈ R n + . Here the absolute value and the comparison are conducted elementwise.
Given the current noisy measurements of the statex N +1 = x(t N +1 ) and the trajectory T N , we first seek to overapproximate the reachable set of the system.

Problem 1 (REACHABLE SET OVER-APPROXIMATION).
Given a set V ⊆ U of admissible control signals, a time step size ∆t > 0, and a maximum number T > N of time steps, compute an over-approximation of the reachable set at time Next, we seek to control the system by computing approximate solutions to the one-step optimal control problem (10).

Problem 2 (APPROXIMATE ONE-STEP OPTIMAL CONTROL).
Given a single finite-horizon trajectory T i−1 for some i ∈ N, i > N , the current measurementx i at time t i , compute, at sampling time t i , approximate solutions of the one-step optimal control problem (10) and characterize the suboptimality of such solutions with respect to the solutions of the optimal control problem if the dynamics were known.
Example 2 (UNICYCLE SYSTEM). Consider, as a running example, the unicycle system with dynamics given bẏ where the components of the state x = [p x , p y , θ] represent, respectively, the position in the x plane, the y plane, and the heading of the unicycle. The components of the control u = [v, ω] represent the speed and the turning rate, respectively. We consider the constraint set U = [−3, 3] × [−π, π]. In the control-affine form (11), we have f = 0, (G) 1,1 = cos(θ), (G) 2,1 = sin(θ), (G) 3,2 = 1, and (G) k,l = 0 otherwise. We assume the dynamics (13) are unknown and are given the loose Lipschitz bounds L f = [0.01, 0.01, 0.01], L G1,1 = L G2,1 = 1.1, L G3,2 = 0.1, and L G k,l = 0 otherwise. These bounds encode the knowledge about the influence of the current state and control signals on the unicycle dynamics. For example, since p x and p y are unchanged when no velocity is applied, we can deduce that (G) 1,2 = (G) 2,2 = 0 and therefore L G1,2 and L G2,2 . The same reasoning allows to say that f = 0, but we do not use such knowledge. Furthermore, we consider the knowledge that the vector field does not depend on its positions. That is, f (x) = f (θ) and G(x) = G(θ).
IV. OVER-APPROXIMATIONS OF THE REACHABLE SET OF UNKNOWN DYNAMICAL SYSTEMS We develop DaTaReach to address Problem 1. It constructs a differential inclusion that contains the unknown vec-tor field. Then, it utilizes the interval Taylor-based method to over-approximate the reachable set of dynamics described by the constructed differential inclusion. We also provide a bound on the time step size that guarantees that an a priori rough enclosure exists. Henceforth, the bound enables to ensure the correctness and termination of DaTaReach.

A. Differential Inclusion based on Lipschitz Continuity
First, to aid the construction of the differential inclusion, we over-approximate f and G at each data point of T N .
for all k ∈ N [1,n] , are the smallest intervals enclosing f (x i ) and G(x i ), respectively, given only the data point, F i , and G i .
The proof for Lemma 1 is derived exactly as in our previous work [10] for noiseless data, where we replace the unknown exact derivativeẋ i with the uncertain quantitiesx i + [−η,η] and propagate the result with interval arithmetic. Intuitively, we use the to remove from F i and G i some values of f (x i ) and G(x i ) that do not satisfy the constraint. Note that, in practise, |ẋ i − x i | ≤η holds since measurements of the state derivative are usually obtained from derivatives of the noisy measurements of the statex i . We also provide an illustration of Lemma 1 in Appendix A for the unicycle system.
Next, we provide functions to globally over-approximate f and G, via uncertain knowledge of f and G at the data points.
and the bounds L f and L G , the functions f : X → IR n and G : X → IR n×m , given for all k ∈ N [1,n] and l ∈ N [1,m] by We provide a proof for Lemma 2 in our previous work [10], where we use interval arithmetic with the definition of the upper bounds on the Lipschitz constants to over-approximate the unknown functions f and G. Note that interval extensions of f and G can be obtained by replacing occurences of · 2 Algorithm 1 Optimal over-approximation of the values of f and G at each data point of a finite-horizon trajectory. Input: Single trajectory T N , sufficiently large M > 0, upper bounds on the Lipschitz constants, noise boundη.
Execute lines 3-6 with E N instead of E i−1 on line 4 9: while E N is not invariant 10: return E N with its interval extension given in Example 1. In the rest of the paper, when the input of f and G are intervals, we use such interval extensions.
Finally, we develop Algorithm 1 that utilizes the trajectory T N to compute the set E N required in Lemma 2 to overapproximate f and G. Then, we construct a differential inclusion containing the unknown vector field.
Theorem 1 (DIFFERENTIAL INCLUSION). Given a trajectory T N , the bounds L f and L G , and the noise boundη, the dynamics (11) are contained in the differential inclusioṅ where the functions f : X → IR n and G : X → IR n×m are obtained by (15) with E N being the output of Algorithm 1.
A proof of Theorem 1 is provided in our previous work [10], where we show that the output E N of Algorithm 1 is such that for all ( Therefore, Lemma 2 enables to conclude. Figure 2 shows that the differential inclusion holds on the unicycle system using a randomly generated trajectory T 15 . Time evolution ofṗx and its over-approximation (f (x) + G(x)u) 1 for the unicycle system of Example 2 in a noiseless setting. We generate the trajectory corresponding toẋ(t) using a randomly generated piecewiseconstant input u(t) ∈ U for t ≤ t 15 = 1.5s, and u(t) = [1, cos(6(t − t 15 ))] for t ∈ [t 15 , 4].
Remark 1 (PERSISTENT EXCITATION). The quality of the differential inclusion of Theorem 1 depends on how much information on f and G can be obtained from T N . Thus, if the trajectory is not diverse, it is likely impossible to obtain tight differential inclusions. For example, consider one of the corner cases where u i = 0 for every data point in T N . In this case, it is impossible to retrieve any information on G. When the goal is to control the system, control values can be synthesized to diversify the trajectory and obtain tight differential inclusions that help future control decisions. For example, excitation-based control (i.e., control inputs that are either zero or nonzero only on a single axis) can be performed to learn the dynamics.

B. Bound on the Time
Step Size for the Existence of an a Priori Rough Enclosure In this section, we consider the problem of finding an a priori rough enclosure S i solution of the fixed-point equation (1) represents the range of a function over an interval domain, Specifically, we give a bound on the time step size ∆t such that S i exists and provide an analytic expression for S i . Recall that S i and ∆t enable to compute an over-approximation of the reachable set when using interval Taylor-based methods.
Theorem 2 (EXPLICIT A PRIORI ROUGH ENCLOSURE). Assume that the time step size ∆t satisfies Then, the set S i given by is an a priori rough enclosure solution of the fixed-point equation (17). The functions f and G in (19) can be any known interval extensions of the unknown functions f and G.
A proof of Theorem 2 is based on several lemmas and is provided in Appendix A. Intuitively, we use the bounds on the Lipschitz constants to provide explicit bounds on the rate of changes of solutions of the differential equation (11), which provide intuition on how to pick a set S i . However, the condition involves β i which is also dependent on ∆t, as v([t i , t i + ∆t]) depends on ∆t. The following corollary gives an explicit bound on ∆t based on the constraint set U. Corollary 1. Under the notation of Theorem 2, assume that where Then, S i given by (19) is a solution of (17).
Proof. This is immediate by observing that 1 Hence, Theorem 2 applies.

C. Interval Taylor-Based Method for Differential Inclusions
We compute an over-approximation of the reachable set of the dynamics described by the differential inclusion (16). The obtained over-approximation is also an over-approximation of the reachable set of the unknown system. Theorem 3 provides a closed-form expression for such over-approximating sets.
Given a trajectory T N , a set V ⊆ U of piecewise-C Du control signals for D u ≥ 1, the bounds L f and L G , the set R + i ∈ IR n of states at time t i , assume that ∆t > 0 satisfies (18). Then, an over-approximation R + i+1 of the reachable set at t i + ∆t of dynamics described by the differential inclusioṅ where ) denote the intervals that contain the range of all the control signals in V and their first derivative on [t i , t i + ∆t] and at t i . The interval matrices J f ∈ IR n×n and J G ∈ IR n×m×n are interval extensions of the Jacobian of f and G, given by and the set S i can be either obtained by (19) of Theorem 2 or by solving the fixed-point equation A proof of Theorem 3 is provided in our previous work [10], where we use a Taylor expansion (6) of order D = 2 with an over-approximation of the unknown Jacobian of f and G obtained from the Lipschitz bounds. Besides, Theorem 2 ensures the existence of an a priori rough enclosure S i .

Remark 2. Theorem 3 implicitly relies on the knowledge of the intervals
given solely by the definition of the set V. For such intervals to exist, all the controls u ∈ V must be at least C 1 on the subdomain [t i , t i + ∆t]. Hence, the assumption that the controls in V are piecewise-C Du with D u ≥ 1. However, when V is such that D u = 0, we propose Corollary 2 to over-approximate the reachable set of the system. We will focus on Theorem 3 in the remainder of the paper since, in practise, the control signals of interest are piecewise-constants.
Compute v(t i ), V i , and V Get tighter f (R + i ) and G(R + i ) if any side information 6: Compute S i via (19) or (24), f (R + i ), and G(R + i ) 7: Compute f (S i ) and G(S i ) via Theorem 1 8: Get tighter f (S i ) and G(S i ) if any side information 9: Compute J f and J G via (22) and (23) 10: Get tighter J f and J G if any side information 11: Compute R + i+1 via (21) 12: end for 13: Corollary 2. Under the notation of Theorem 3, assume that D u = 0. Then, an over-approximation R + i+1 of the reachable set of the dynamics described by the differential inclusionẋ ∈ f (x) + G(x)u with u ∈ V is given by

D. Side Information
As motivated in the introduction, prior knowledge on the underlying unknown dynamics might be available. However, Theorem 3, as it is, does not incorporate any side information. We now show how to incorporate a variety of side information to obtain tighter over-approximations of the reachable sets.
Side information 1 (VECTOR FIELD BOUNDS). We are given R f A ∈ IR n and R G A ∈ IR n×m as supersets of the range of f and G, respectively, over a given set A ⊆ X .
Given a set S ⊆ A, tight extensions of f and G over S can be obtained by the update Side information 2 (GRADIENT BOUNDS). We are given bounds on the gradient of some components of f and G. Such side information may include the monotonicity of f and G.
These bounds can be used in equations (22) and (23) to provide tight interval extensions J f and J G . For example, if the function (f ) k is known to be non-decreasing with respect to the variable x p on a set A ⊆ X , then we obtain a tighter R + i+1 by the update This side information expresses properties such as decoupling between states. For example, if the state (x(t)) p does not directly affect (ẋ(t)) k for some p, k ∈ N [1,n] under any control signal in U, we can obtain a tighter over-approximation of the reachable set by setting to zero the intervals (J f ) k,p and (J G ) k,l,p for all l ∈ N [1,m] . For the unicycle system of Example 2, since f and G depends only on the heading θ, the Jacobian terms and (J G ) 3,2,2 must all be set to zero.
Side information 3 (ALGEBRAIC CONSTRAINTS). We are given a differentiable map g : R n → R q satisfying the constraint g(x) ≥ 0 over the state x of the unknown system. Without loss of generality, we consider that g : R n → R. This side information provides bounds on f , G, J f , and J G locally. Specifically, if we derivate the function t → g(x(t)), then we obtain a new constraint w : R n ×R n×m ×R m ×R n → R such that w(f (x), G(x), u, x) ≥ 0. If we compute the partial derivatives of w with respect to all the states x i , then we obtain a new constraint z : R n×n × R n×m×n × R m × R n → R n such that z( ∂f ∂x (x), ∂G ∂x (x), u, x) ≥ 0. The obtained constraints w and z can be incorporated in the computation of R + i+1 through the framework of contractors [27], [28]. Specifically, given a set A of states, a set V of control values, and the overapproximating sets f (A), G(A), J f , and J G , the methods based on contractor programming [27] compute the intervals In other words, it contracts f (A), G(A), J f , and J G by pruning out some values that do not satisfy the constraints imposed by w and z. In the computation of R + i+1 , the set A is either R + i or S i , and the set V is V i . Side information 4 (PARTIAL DYNAMICS KNOWLEDGE). We are given terms of some components of the vector field of (11). That is, we consider that the dynamics are in the forṁ where the functions f kn and G kn are known while f ukn and G ukn are unknown. Such side information is usually obtained from the application of elementary laws of physics.
The new functions f and G are given by f = f kn + f ukn and G = G kn + G ukn , where the functions f kn and f ukn are interval extensions of known functions f kn and G kn , respectively, and the functions f ukn and G ukn are obtained by Theorem 1 using L f ukn , L G ukn , and the new trajectory Furthermore, the new Jacobian terms in the computation of R + i+1 are given by J f = ∂fkn ∂x (S i ) + J f ukn and (J G ) k,l,p = ∂(Gkn)k,l ∂xp (S i ) + (J G ukn ) k,l,p , respectively.
DaTaReach (Algorithm 2) addresses Problem 1 by combining Theorem 3 and the available side information. We demonstrate the value of additional side information on the unicycle system. Figure 3 shows that, as expected, the overapproximation becomes tighter with richer side information.

V. CONTROL SYNTHESIS FOR UNKNOWN DYNAMICS
We develop DaTaControl to address Problem 2. It computes approximate solutions to the one-step optimal control problem (10) using the over-approximation of the reachable set. Specifically, it replaces the unknown next state x i+1 in problem (10) with the corresponding over-approximation of the reachable set R + i+1 , which may be computed using Theorem 3. Unfortunately, R + i+1 is a nonconvex function of the control signal u, the decision variable. Therefore, it constructs an over-approximation of R + i+1 , in Theorem 4, that is convex in the given constant control signal u applied in the time interval [t i , t i+1 ]. The new over-approximation provides tractable approximations of the optimization problem (10).

A. Approximate One-Step Optimal Control
First, Theorem 4 provides a control-affine linearization of the over-approximation R + i+1 given by Theorem 3. Theorem 4. Under the notation of Theorem 3 and given a constant control signal u : t → u i applied between t i and t i+1 = t i + ∆t with u i ∈ U, the reachable set R i+1 satisfies where the intervals A − i , A + i , and B i are given by where (J T G ) k,p,l = (J G ) k,l,p for k, p ∈ N [1,n] and l ∈ N [1,m] , and the a priori rough enclosure S i is either obtained by (19) with v([t i , t i+1 ]) replaced by U or is a solution of A proof of Theorem 4 is provided in our previous work [10], where we linearize the quadratic term in u of the closedform expression (21). Such linearization enables to have a control-affine over-approximation at time t i + ∆t. Intuitively, the linearization comes from over-approximating the quadratic term (J G u i )(G(S i )u i ) of Theorem 3 by linear terms (J G U)(G(S i )u i ) and (J G u i )(G(S i )U). Thus, we obtain two different over-approximations characterized by A + and A − , which are finally intersected to provide a tighter set.
Next, we use the control-affine linearization of Theorem 4 to propose two convex optimization problems which solutions approximate the one-step optimal control problem (10). The first problem, called optimistic control problem, is given by where the goal is to minimize the best possible cost value over all possible state x i+1 in the over-approximation (26). The second problem, called idealistic control problem, is an idealistic approximation of (10) given by where the goal is to minimize the cost associated to a specific trajectory x i+1 = b ide i +A ide i u i in the over-approximation (11), idealistically considered as the unknown next state evolution. Specifically, we choose x i+1 according to the middle point between a point Remark 3. Another possible approximation of the one-step optimal control problem (10) is a pessimistic control problem where, in contrast to the optimistic control problem, the goal is to minimize the worst possible one-step cost value over all possible state x i+1 in the over-approximation of the reachable set. Such a pessimistic control problem can be formulated as a semidefinite program using tools from robust optimization.
Finally, we provide, in Theorem 5, bounds on the suboptimality of the optimistic and idealistic control with respect to the optimal control when the cost function c is quadratic.
Assumption 2 (QUADRATIC ONE-STEP COST). We consider that the one-step cost function c is a convex quadratic function, where q ∈ R n , r ∈ R m , Q = Q T ∈ R n×n , R = R T ∈ R m×m , and S ∈ R n×m .
Theorem 5 (SUBOPTIMALITY BOUND). Let c i , c opt i , and c ide i be the optimal cost of the one-step optimal control problem (10), the optimistic control problem (29), and the idealistic control problem (30), respectively, at the sampling time t i . The following inequality holds Algorithm 3 DaTaControl: Approximate one-step optimal control solution at t i ≥ t N Input: Single trajectory T i−1 of length i > N , ime step size ∆t satisfying (20), Lipschitz bounds, noise bounds η ≥ 0 andη ≥ 0, one-step cost c, current statex i =x(t i ).
Optional: Any of the side information 1-4. Output: Constant controlû i ∈ U to apply between t i and t i + ∆t that approximates a solution of (10).
We provide a proof for Theorem 5 in our previous work [10] where we use the local Lipschitz property of the cost function c to characterize its variations over the over-approximation set (26).

Remark 4.
The main term of the suboptimality bound in Theorem 5 is directly related to the width of the overapproximation of the reachable set at time t i + ∆t. Thus, as the over-approximation of the reachable set becomes tighter, the gap with the unknown optimal cost decreases. Specifically, it is straightforward to see that the more data and richer side information are used in the computation of B i , A + i , and A − i , the tighter the widths w(B i ), w(A + i ), and w(A − i ) become.

B. Real-Time Approximate Control Synthesis
In this section, we investigate efficient algorithms to compute solutions of the optimistic control problem (29) and the idealistic control problem (30) when the one-step cost function c satisfies Assumption 2. We provide a bound on the number of elementary operations required by DaTaControl to synthesize near-optimal control values. Such a bound makes DaTaControl a good candidate for real-time control of systems with a priori unknown dynamics.

Lemma 3 (APPROXIMATE QUADRATIC PROBLEMS).
Without loss of generality, consider that either (U) l ≥ 0 or (U) l ≤ 0 for a given l ∈ N [1,m] . The optimistic control problem (29) can be reformulated as the quadratic programming problem minimize ui∈U ,xi+1∈X and the idealistic control problem (30) as the convex, boxconstrained quadratic programming problem where Q i , p i , and q i are given by the matrices A l+ i and A s+ i are such that, for all k ∈ N [1,n] and l ∈ N [1,m]  Remark 5. Lemma 3 assumes that the components of U are either nonnengative or nonpositive in order to simplify the notation and provide a simple expression for the optimistic control problem. In the general case, U can be partitioned into subsets where each component is either nonnegative or nonpositive. Then, we solve the optimistic control problem on each partition and return the control value from the partition with the minimum cost. It is important to note that the idealistic control problem (35) does not require any partitioning of U and can therefore be solved as a one-shot quadratic problem. Besides, if m i ≤ m is the number of components of U having an indefinite sign, then solving the optimistic control problem (29) requires to solve 2 mi quadratic programming problems similar to (34). As a consequence, the number of problems drastically reduces when the intervals (U) l are known to be nonnegative or nonpositive a priori.
We seek for algorithms to compute solutions of (34) and (35) while providing bounds on the number of elementary operations to obtain such solutions. The number of elementary operations, referred as flops, is obtained by counting the number of operations +, −, * , /, max, min, and (·) executed by the algorithms. When the cost function c is assumed to be strongly convex, i.e., either Q is positive definite or R is positive definite, the authors of [29] develop an algorithm [29, Algorithm 1] based on applying the fast gradient-projection method on the dual of the optimization problem. The algorithm can be directly used to solve both the optimistic and idealistic control problem. The rate of convergence of such an algorithm is proportional to 1/ √ , where is the desired accuracy for optimality. Furthermore, the authors provide a maximum number of flops, cubic polynomial in both n and m, required by the algorithm to terminate. We consider the general case where c might not be strongly convex and develop an algorithm inspired from [30] to solve the idealistic control problem with a linear rate of convergence to a solution, i.e. proportional to ln 1/ . We deduce a maximum number of flops, quadratic in both n and m and linear in the length of the Algorithm 4 AdaRES: Accelerated proximal gradient with adaptive restart for optimization problems satisfying the hypothesis of Lemma 4 Input: Desired accuracy for the optimality, a point y 0 ∈ R p , smoothness parameter L, an initial estimation µ 0 of µ P (y 0 ), and a function Π P for projection onto P. Output:ŷ ∈ P such that l(ŷ) − l ≤ 1: s ← −1, t s ← 0, y s,ts ← y 0 , y s+1,0 ← Π P (y s,ts ) 2: do 3: for k ∈ N [0,...,Ks−1] do 9: x k+1 ← Π P (z k − 1 L ∇l(z k )) 10: end for 14: y s,t+1 ← z Ks , t ← t + 1 15: while trajectory T p , required by DaTaControl to compute nearoptimal control values when the idealistic formulation is used. THEOREM 2). Consider the convex optimization problem inf y∈P l(y), where P ⊆ R p is a convex set and l : R p → R is a differentiable convex function. Let l * ∈ R be the optimal cost, P the set of minimizers of the problem, and assume that there exists L > 0 and µ P (y 0 ) > 0 for any y 0 ∈ R p such that l satisfies the smoothness and local error bound assumption given respectively by for all y, w ∈ R p and for all z ∈ P such that l(z) ≤ l(y 0 ). Then, for any µ 0 > 0 and > 0, Algorithm 4 returns a solutionŷ ∈ P such that l(ŷ) − l ≤ with a total number of iterationsN satisfyinĝ µ P (y 0 ) ln 2(l(y 0 ) − l µ P (y 0 ) , otherwise. (41) Remark 6. Lemma 4 guarantees a linear rate of convergence of AdaRES without knowing explicitly the parameter µ P (y 0 ).
Specifically, AdaRES achieves this rate of convergence by being able to check at iteration s whether the current estimate µ s of µ P (y 0 ) is too large and finally by refining such estimate. Therefore, AdaRES can also be used to estimate a lower bound on µ P (y 0 ) and thus providing the guarantees (41).
Note that the number of iterationsN represents also the total number of calls to the projection Π P (·) onto P.
We exploit Lemma 4 to provide a maximum number of flops required by DaTaControl to terminate.
Theorem 6 (NUMBER OF FLOPS BY DATACO N T R O L). At each sampling time t i ≥ t N , given a trajectory T i of length i ≥ N , if DaTaControl based on the idealistic control problem (30) uses Algorithm 5 with µ 0 ≤ µ U (v 0 ) and without any side information 1-4, then DaTaControl takes at most N ide flops to compute near-optimal control values with N ide roughy given by + 2i(8n + 9)(m + 1) + 24n 2 m + 8n 2 + 8m 2 where c max = |c(X , U, X )| is the maximum value in absolute value taken by the cost function over the range X and U. Furthermore, if any of the side information 1-4 is provided, then only a maximum of R(n, m) flops will be added to N ide where R is a polynomial of maximum degree 2.
We provide in Appendix A a proof of Theorem 6 where we show that (30) satisfies the assumption of Lemma 4 in order to bound the number of flops required by DaTaControl.
Remark 7. Even though the optimistic control problem (29) satisfies the assumptions of Lemma 4, computing the projection Π P is as hard as solving directly (29). Hence, we choose to focus on the idealistic control problem. Note that µ 0 can be empirically estimated offline by solving the idealistic control problem on many instances of the same problem.

Time step
Cost function Fig. 4. DaTaControl outperforms state-of-the-art data-driven approaches in terms of computation and optimality of the computed control when applied to a scenario with no noise in the measurements. Opt. traj. corresponds to the one-step optimal control with the known dynamics, and the solid black line on the far right figure represents the optimality threshold for reaching the goal. We emphasize that the extremely low computation time of DeePC at the beginning is due to randomly generated control inputs until the satisfaction of the persistent excitation condition.
and SINDYc [6]. These algorithms also solve the one-step optimal control problem (10) approximately. CGP-LCB assumes that the unknown one-step approximate cost function x i , u i )) to minimize is described by a Gaussian process, while C2Opt assumes that C has Lipschitz continuous gradients. On the other hand, SINDYc uses the limited data to perform sparse identification of the dynamics over a library of nonlinear functions, which then permits an approximate solution to (10) via numerical solvers. DeePC directly computes control values by solving an N -step model predictive control problem without the system identification step. The experiments on the unicycle system, a quadrotor, and an aircraft system demonstrate that DaTaControl outperforms these algorithms both on the computation time at each sampling time, the robustness to noise, and the suboptimality of the control. We provide the code for reproducibility of the results of this paper at the link: https://github.com/wuwushrek/datacontrolreach.git.

A. The Unicycle System
We consider the problem of driving the unicycle of Example 2 to the origin. We encode this control objective using the one-step cost function c(x i , u i , x i+1 ) = 0.5 x i+1 2 2 . We chose the time step size ∆t = 0.1s and generated a random initial trajectory T 10 starting from the state x 0 = [−2, −2.5, π/2] and such that the unicycle goes away from the target. Note that the time step size ∆t = 0.1s is within the limit enforced by (20) as 1/( √ nβ inf ) ≈ 0.123. Hence, we used the a priori rough enclosure given by Theorem 2.
We applied DaTaControl with the side information discussed in Example 2. We use DaTaControl with Algorithm 5 as the baseline to solve the idealistic control problem. The weights w + and w − were chosen randomly, and µ 0 = 1.0. We used the default parameters in GPyOpt [31] when implementing CGP-LCB. We chose the Lipschitz constant for the gradient L = 10 and trade-off hyperparameter α = 1/2 for C2Opt [8]. For SINDYc, we considered monomials (up to degree 6), sines and cosines of the state, and the products of these functions with the velocity v and the turning rate ω as the library functions. To perform sparse identification, we swept the regularization parameter [6] λ ∈ {10 p : p ∈ N [−6,5] } and rounded-down the coefficients smaller than 10 −3 to zero. For DeePC, we use the data-driven model predictive control scheme (8) in [7] with the slack variable and the regularization in the cost to account for the nonlinearities in the unicycle dynamics. After extensive tuning, we use the parameters λ y = 10 4 , λ g = 1, T ini = 8, and the horizon control N = 10. We note that finding working hyperparameters for nonlinear dynamics was a hard task and the hyperparameters did not work for different initializations of the unicycle. We relaxed the target state to the 0.1-sublevel set of the one-step cost function as the stopping criteria for the algorithms. Figure 4 shows the trajectories and the evolutions of the cost functions for the different algorithms under noiseless measurements. DaTaControl performs significantly better than GPyOpt, SINDYc, DeePC, and C2Opt. It reaches the origin in fewer time steps and significantly lower computation time. Instead, SINDYc failed as it gets stuck close to the origin and is unable to further improve its one-step cost function.
Additionally, we compare DaTaControl with the datadriven baselines under the same conditions as the previous experiment but with different levels of noise in the measurements of the state and the state derivatives. Figure 5 shows that DaTaControl has improved robustness to noise.

B. Damaged Aircraft System
In this section, we consider the problem of controlling an aircraft with nonlinearities introduced as damage in the pitch.
The model above is based on the linearized dynamics of an undamaged Boeing 747 in landing configuration at sea level [32], [33], where the damages are introduced via nonlinear factors in the dynamics [13]. The state of the aircraft is given by x = [w l , w v , q, θ, h] and the control input by u = [δ e , δ t ]. w l is the deviation from the steady-state longitudinal speed, w v the downwards-pointing vertical speed, q the pitch rate, θ the pitch angle and h the altitude. The control inputs δ e and δ t are the deviations from the elevator and thrust values.
The introduced damages emulate damages to the horizontal stabilizer of the aircraft, which result in the pitch rate behaving more erratically. We seek to control the pitch angle to a setpoint θ sp = 5 while assuming that the dynamics are unknown. We encode this control objective using the cost function c(x i , u i , x i+1 ) = 0.5((x i+1 ) 4 − θ sp ) 2 . We chose the time step size ∆t = 0.01s and generated a random initial trajectory T 10 starting from the state x 0 = [0, 0, 0, 0, 100.0]. We used DaTaControl, CGP-LCB with the same hyperparameters as for the unicycle dynamics. For SINDYc, we considered monomials (up to degree 2), sines and cosines of the state. We did not investigate C2Opt due to the inability to compute the gradient of C and the fact that it is a strong assumption to assume measurements of such gradient. For DaTaControl, we consider the elementary side informatioṅ θ = q. Furthermore, we consider that G(x) = G(q, w v ) and f (x) = f (w l , w v , q, θ) as the longitudinal speed and the altitude do not directly affect the derivative of the states. We considered the loose Lipschitz bounds L f = [0.4, 3, 4, 1, 3], L G3,2 = 0.5, and L G k,l = 0.01 otherwise.
We first consider a linear version of the aircraft dynamics where we remove the nonlinear terms sin(x 1 ), cos(x 1 ), and sin(x 2 ) introduced for the damage. The objective here is to show that DaTaControl is more performant on linear systems than state-of-the art approaches based on behavioral systems theory foundation, e.g., DeePC. For linear systems, DeePC does not use slack variables and considers the parameters λ y = 0, λ g = 0 and T ini = 8 as instructed in [7]. Figure 6 shows the near-optimality of DaTaControl, while CGP-LCB fails to stabilize the pitch in a noise-free scenario.
Then, we consider the nonlinear aircraft dynamics above with bounded noise on the measurements on the state and the state derivatives. In this scenario, we use the parameters λ y = 10 3 , λ g = 10 to account for nonlinearities when using DeePC. Figure 7 shows that DeePC is less performant on the nonlinear dynamics, while SINDYc and DeePC fail when there is noise in the measurements. We did not compare with CGP-LCB due to its extremely high computation time for solving the control problem at each time step.

C. Quadrotor System
Consider a quadrotor with control-affine dynamics [34]  We assume that the dynamics of the quadrotor are unknown and consider the problems of controlling v x to a setpoint v sp x = 5 and p y to a setpoint p sp y = 8. We encode these control objectives, respectively, using the cost functions We chose the time step size ∆t = 0.01s and generated a random initial trajectory T 10 starting from the state x 0 = [0, 0, 5, 0, 0, 0]. We consider the same hyperparameters for Algorithm 5 as in the unicycle case. We applied DaTaControl with the side informationṗ x = v x ,ṗ y = v y , andφ = ω. Such extra knowledge is obtained from elementary laws of physics. Furthermore, DaTaControl considered the loose Lipschitz bounds L f2 = L f4 = 0.3, L f6 = 0.9, L G6,1 = L G6,2 = 0.01, L G2,1 = L G2,2 = L G4,1 = L G4,2 = L G6,1 = L G6,2 = 0.9, and uses the side information G(x) = G(φ) and f (x) = f (v x , v y , w). We used the default parameters of GPyOpt when implementing CGP-LCB. For SINDYc, we considered monomials (up to degree 1), sines and cosines of the state, and the products of these functions with the T 1 and T 2 as the library functions. We did not investigate C2Opt due to the inability to compute the gradient of C. Figure 8 shows the near-optimality of DaTaControl while CGP-LCB and SINDYc fail to reach the setpoints. Furthermore, the figure demonstrates that DaTaControl can achieve near-optimal, near-real-time control of the vertical position and horizontal speed. We justify the suboptimality of Opt. traj. by the fact that the one-step optimal control problem is highly nonlinear, which makes possible the synthesis of local optimum solutions by the numerical solvers.

VII. CONCLUSION
We developed two data-driven algorithms, DaTaReach and DaTaControl, for on-the-fly over-approximation of the reachable set and constrained near-optimal, real-time control of systems with unknown dynamics. These algorithms are suitable for scenarios with severely limited noisy data and can take advantage of various side information on the underlying a priori unknown dynamics. The numerical experiments demonstrate the efficacy of the algorithms over existing approaches both in terms of the suboptimality of the control, robustness to noise, and the computation time to synthesize control values.

APPENDIX A
In this appendix, we provide the proofs for all lemmas and theorems presented in this paper.

PROOF OF THEOREM 2
First, we provide a bound on the maximum variation of the function h given upper bounds on the Lipschitz constants of f and G.
Lemma 5 (MAXIMUM VARIATION OF h). Let L f and L G be upper bounds on the Lipschitz constants of f and G. For all x 1 , x 2 ∈ X and u 1 ∈ U, we have that Proof. This is a direct application of the definition of the Lipschitz bounds of f k and G k,l . Specifically, we have that for all k ∈ N [1,n] . The last inequality is obtained by the definitions of L f k and L G k,l . Hence, we obtain (43).
with γ(s) ≥ 0 for all s ∈ [t 0 , T ]. Then, for all t ∈ [t 0 , T ], Next, we combine the results of Lemma 5 and Lemma 6 to provide a bound on the variations of points in a trajectory. (11)). Let x be a continuous-time signal satisfying (11) for a given control signal u ∈ V. For all t ∈ [t i , t i + ∆t], we have that

Lemma 7 (VARIATION OF TRAJECTORIES OF
where the parameters β i ∈ R + and α i ∈ R n are given by Proof. For all t ∈ T i = [t i , t i + ∆t], we have that We obtain (50) since x is a solution of (11). The passage from (50) to (51) results from applying the triangle inequality. We use Lemma 5 with x 1 = x(s), x 2 = x(t i ), and u 1 = u(s) to obtain (52). We obtain (53) by showing that the inequalities β(u(s)) ≤ β i and h(x(t i ), u(s)) 2 ≤ α i 2 hold for all s ∈ T i . Specifically, the former inequality holds by upperbounding |(u(s)) l | in β(u(s)) with |(v([t i , t i + ∆t])) l | ≥ sup s∈Ti |(u(s)) l |. The latter inequality holds since The inequality (53) satisfies the condition (46) of Lemma 6 with Ψ(t) = x(t)−x(t i ) 2 , α(t) = α i 2 (t−t i ), and γ(t) = β i . Thus, we have that By integration by parts, we have that Finally, by combining (55) and (54), we obtain (47).
Finally, we provide the proof of Theorem 2 below. The expression (19) is derived by scaling adequately the bound (47) from Lemma 7. Specifically, we seek for an a priori rough enclosure S i such that where µ > 0 is a parameter to find in order for S i to satisfy the fixed-point equation (17) and α i is given by We over-approximate the set R(h, S i × v([t i , t i + ∆t])) as a function of ∆t, β ti,ti+∆t , and h = f + Gu. For all x i ∈ R + i , u i ∈ v([t i , t i + ∆t]), and s i ∈ S i , we have that The inequality (57) comes from Lemma 5 with β(u i ) upperbounded by β i and the definition of S i (56) provides an upper bound on s i − x i that yields (58). Additionally, the inequality (58) implies that Hence, S i from (56) solves the fixed-point equation (17) if For notation brevity, let c 1 = α i (e βi∆t − 1). Observe that [0, ∆t][A, A] = [min(0, A), max(0, A)]. We use the observation to find µ > 0 such that the inclusion (59) holds. That is, the inequalities hold for all k ∈ N [1,n] . Therefore, for a step size ∆t satisfying (18), µ given by satisfies the above inequalities. Thus, the inclusion (59) holds and the set S i is solution of the fixed-point equation (17). By replacing µ in S i given by (56), we obtain (19). We can deduce that A + i u i = A s+ i u i and A + i u i = A l+ i u i . Similarly, it is easy to prove that A − i u i = A s− i u i and A − i u i = A l− i u i . As a consequence, the optimistic control problem (29) can be immediately reformulated as the convex quadratic programming problem (34).
Additionally, we have that Hence we obtain the cost function of the idealistic control problem (35). Therefore, it remains to prove that the matrix Q i is positive semidefinite. Let prove that u T i Q i u i ≥ 0 for all u i ∈ R m . Since by assumption, the cost function is convex, then the matrix implied in the quadratic term is positive semidefinite. Hence, we can deduce that for all u i ∈ R m . Expanding the expression above immediately provides that u T i Q i u i ≥ 0 for all u i ∈ R m . PROOF OF THEOREM 6 First, using the notation of Lemma 4, we have that l(y) = 0.5y T Q i y + q T i y + p i and P = U for the idealistic control problem. Thus, l satisfies condition (39) with any parameter L ≥ λ max (Q i ) where λ max (Q i ) is the maximum eigen value of Q i . Recall that since we have x − y 2 2 , then Lemma 3.4 in [36] implies that the condition (39) holds. Since λ max (Q i ) ≤ √ m Q i ∞ , we pick L = √ m Q i ∞ . The local error bound assumption holds due to Theorem 2.1 in [37] where the authors prove that for a function h satisfying h(x) = g(Ax)+c T x on a convex polyhedra P with g strongly convex, there exists µ P (x 0 ) for any x 0 such that the condition (40) holds. As a consequence, since we assume that µ 0 ≤ µ P (x 0 ), Lemma 4 provides that the number of iterationsN required by AdaRES is upper bounded as followŝ As a consequence, since the computation of Π P (·) is simple and requires at most 3m comparisons, and ∇l requires 2m 2 + m flops, we can deduce that the maximum number of flops required by AdaRES to terminate is roughly given by N max (2m 2 + 16m + 12). Such a number is obtained by considering the worst case where each line of Algorithm 4 is executed N max times even though in reality, only lines 9-12 of AdaRES will be executed maximum N max times.
Finally, since computing f and G requires N (8n + 9) and N (8n+9)m respectively, we can deduce that, without any side information 1-4, lines 1-9 of DaTaControl require at most 2N (8n + 9)(m + 1) + 16n 2 m + 8n 2 + 2m 2 + 6nm + 24n + 8. Specifically, it is easy to check that S i using Theorem 2 can be computed using 6mn + 6n + 4 flops. The matrices B i , A + i and A − i can be computed using 4n 2 + 6n flops for B i and 8n 2 m + 3nm + m 2 flops for each A + i and A − i . Hence, since computing A ide i and B ide requires 8n + 8nm flops, we can compute Q i and q i using 6n 2 m + 6m 2 and (2n 2 + 4n + 2)m flops, respectively. Thus, after adding the number of flops for AdaRES, lines 1-9 of DaTaControl, A ide i and B ide , and Q i and q i , we obtain the bound N ide . ILLUSTRATION FOR LEMMA 1 Consider the unicycle system of Example 2, and assume that (G) 1,2 is unknown. Given the data point (x i ,ẋ i , u i ) where Therefore, Lemma 1 provides the contractions (C F ) 1 ⊂ (F i ) 1 , (C Gi ) 1,1 ⊂ (G i ) 1,1 , and (C Gi ) 1,2 ⊂ (G i ) 1,2 .