Anas Makdesi
email: anas.makdesi@l2s.centralesupelec.fr

Senior Member, IEEE Antoine Girard
email: toine.girard@l2s.centralesupelec.fr

Laurent Fribourg
email: fribourg@lsv.ens-cachan.fr

Data-Driven Models of Monotone Systems

Keywords: Monotone maps, Monotone systems, Data-driven models, Data-driven abstraction, Symbolic control

In this paper, we consider the problem of computing from data guaranteed set-valued overapproximations of unknown monotone functions with additive disturbances. We provide a characterization of a simulating map that provably contains all monotone functions that are consistent with the data. This map is also minimal in the sense that any set-valued map containing all consistent monotone functions would also include the map we are proposing. We show that this minimal simulating map is interval-valued and admits a simple construction on a finite partition induced by the data. As the complexity of the partition increases with the amount of data, we also consider the problem of computing minimal interval-valued simulating maps defined on partitions that are fixed a priori. We present an efficient algorithm for their computation. We then use those data-driven over-approximations to build models for partially unknown systems where the unknown part is monotone. The resulting models are used to construct finite-state symbolic abstractions, paving the way for discrete controller synthesis methods to be applied. We extend our approach to handle systems with bounded derivatives and introduce an algorithm to calculate the bounds on those derivatives and on the disturbances from the data. We present several numerical experiments to test the performance of the introduced method and show that the data-driven abstractions are suitable for controller synthesis purposes.

I. INTRODUCTION

The recent advancements in data acquisition tools, coupled with the ability to efficiently deal with the ever-increasing amount of data harnessed by those tools, gave rise to many new data-driven approaches and techniques in control theory. Those approaches have the benefit of being able to tackle unknown or hard-to-model systems, which justifies the increasing interest in them. While some methods rely on finding the controller directly from the data [12, and references therein], others depend on finding data-driven models, which then can be used to find the controllers. A survey of both methods can be found in [START_REF] Brunton | Data-driven science and engineering: Machine learning, dynamical systems, and control[END_REF]. Learning the dynamics to find systems' models leverages the ability to use well-established, wellstudied approaches to find controllers enforcing the desired behaviors. Moreover, data-driven models can sometimes come with formal guarantees regarding the approximation of the real system. Such guaranteed models can be categorized as stochastic or robust [11, and references therein]. The present work follows the latter approaches, where the aim is to find bounds on the function representing the system (nonparametric approach) or bounds on some unknown parameters of the function representing the system (parametric approach).

The set membership approach offers an example of how we can find those robust models. In [START_REF] Novara | Set membership identification of nonlinear systems[END_REF], the set membership approach is used for the identification of nonlinear systems. Set membership approaches are used in [START_REF] Canale | Nonlinear model predictive control from data: a set membership approach[END_REF] to synthesize nonlinear model predictive controllers. In contrast to [START_REF] Canale | Nonlinear model predictive control from data: a set membership approach[END_REF], the introduced approach does not assume anything about the stability of the unknown systems. Also, the number of points that can be handled is greater than the number addressed in the set membership approach. Another example of robust models can be found in [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF], where interval observers for partially unknown mixed-monotone nonlinear systems are found, assuming the unknown part is Lipschitz continuous. Lipschitz dynamical systems with known bounds on the Lipschitz constants are also studied in [START_REF] Sadraddini | Formal guarantees in data-driven model identification and control synthesis[END_REF] to reach piecewise affine set-valued models. In contrast to [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF] and [START_REF] Sadraddini | Formal guarantees in data-driven model identification and control synthesis[END_REF], we do not consider the case of Lipschitz unknown functions. Instead, we start by studying monotone functions. Relying on the monotonicity assumption, we developed an efficient method to handle a large number of data points. Monotone functions can be found in many applications, such as adaptive cruise control [START_REF] Ivanova | Lazy controller synthesis for monotone transition systems and directed safety specifications[END_REF] or power networks [START_REF] Zonetti | Decentralized monotonicity-based voltage control of dc microgrids with zip loads[END_REF]. In [START_REF] Milani Fard | Fast and flexible monotonic functions with ensembles of lattices[END_REF], models that enforce a notion of monotonicity are learned following the argument that, for many machine learning problems, some inputs relate to the output monotonically, such as house pricing. Data-driven approaches are used in [START_REF] Kawano | Data-driven model reduction of monotone systems by nonlinear DC gains[END_REF] for monotone systems reduction.

In this work, given a set of data generated from an unknown monotone function with an additive bounded disturbance, we find an over-approximating set-valued map, which we call a simulating map, guaranteed to include any monotone function capable of generating the data. Moreover, we look for the tightest map with this guarantee. The resulting model (the minimal simulating map) can be naturally defined as an interval-valued map using a partition of the input space, which only depends on the sampled data, and thus it is non-parametric. We also consider a parametric case where a partition of the input space is given a priori. We find the tightest interval-valued overapproximation on this partition. Although fixing the partition introduces more conservatism to the resulting models, those models can be calculated and handled more efficiently.

Then, we use our over-approximation method to learn models of discrete-time dynamical systems with an unknown monotone part. Those models are then used to find finitestate symbolic models of the systems. Symbolic control is a computational approach to controller synthesis where discrete abstractions (symbolic models) of continuous-state dynamical systems are calculated. Then, discrete controller synthesis techniques [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF], [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF] are used to find correct-bydesign controllers; those controllers are capable of ensuring complex specifications (safety and reachability [START_REF] Girard | Controller synthesis for safety and reachability via approximate bisimulation[END_REF], behaviors described by automata [START_REF] Pola | Control of cyber-physical-systems with logic specifications: a formal methods approach[END_REF], or temporal logic formulas [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF]). In the case of an unknown monotone system, we show that working with the finite-state symbolic models does not add any conservatism compared to working with the tightest datadriven over-approximation models.

Although symbolic control is usually presented as a modeldriven technique, some data-driven approaches have emerged recently. For instance, approaches that only require the ability to sample the system dynamics on a given grid of states and inputs are introduced in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], [START_REF] Meyer | Safety control with performance guarantees of cooperative systems using compositional abstractions[END_REF], or [START_REF] Zamani | Symbolic models for nonlinear control systems without stability assumptions[END_REF]. In [START_REF] Devonport | Symbolic abstractions from data: a PAC learning approach[END_REF], datadriven abstractions within the PAC (probably approximately correct) statistical framework are found. PAC guarantees are also used in [START_REF] Coppola | Data-driven abstractions for verification of deterministic systems[END_REF] to build data-driven abstraction through the use of the scenario approach. The scenario approach is also used in [START_REF] Lavaei | Data-driven synthesis of symbolic abstractions with guaranteed confidence[END_REF] and [START_REF] Kazemi | Data-driven abstraction-based control synthesis[END_REF] to offer probabilistic guarantees on the datadriven built abstractions. The main difference between all the previously mentioned work and the approach introduced here is that, in contrast to our approach, they all give probabilistic guarantees. Data-driven control design with regular language specifications is studied in [START_REF] Pola | Datadriven controller synthesis for abstract systems with regular language specifications[END_REF] for plants described as abstract systems. In comparison to [START_REF] Pola | Datadriven controller synthesis for abstract systems with regular language specifications[END_REF], our approach can handle bounded disturbances and is able to satisfy the specifications robustly, not to a predefined desired accuracy. Unknown parts of nonlinear systems are modeled using Gaussian processes; then, those models are used in building symbolic abstractions in [START_REF] Hashimoto | Learning-based symbolic abstractions for nonlinear control systems[END_REF].

Some results in this paper appeared in preliminary form in the conference papers [START_REF] Makdesi | Efficient data-driven abstraction of monotone systems with disturbances[END_REF], [START_REF] Makdesi | Data-driven abstraction of monotone systems[END_REF], where we focused on the case of pure monotone systems. The present paper brings a globally improved presentation of our approach with a redesigned concept of minimal simulating maps and complete proofs of theoretical results (whereas only sketches of proofs, in the undisturbed case, were provided in our previous works). In addition, we introduce a new efficient algorithm to compute our set-valued over-approximations. We also generalize the approach to the case of general systems with an unknown monotone part. Using ideas inspired by those introduced to study mixed-monotone systems in [START_REF] Yang | On sufficient conditions for mixed monotonicity[END_REF], we extend our approach to study systems with bounded derivatives. Finally, we use the scenario approach, developed for robust control design in [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF], to estimate upper and lower bounds on the partial derivatives of the unknown function generating the data. Finally, we validate our results with more extensive numerical experiments.

This paper is organized as follows. In Section II, we introduce the problem of over-approximating monotone maps, give the necessary definitions, and then compute the optimal solution, and optimal solution given a fixed partition a priori.

In the case of fixed partition, we provide an efficient algorithm to calculate the over-approximation. In Section III, we use the introduced over-approximation to build data-driven models for systems with unknown monotone parts. Then, we use those data-driven models to find a finite-state symbolic representation of the system. We also show how we can study systems with bounded derivatives in this section. In Section IV, we provide a way to find the upper and lower bounds on the partial derivatives of a function using data generated by this function. Finally, Section V is dedicated to presenting some numerical examples.

Notations: R = [-∞, +∞] is the set of extended real numbers. We use bold lowercase letters to represent vectors, e.g. z ∈ R n ; subscripts are used to differentiate between multiple vectors z i , whereas normal lowercase letters with superscripts z i are used to denote the i th component of a vector z. Given two vectors z 1 , z 2 ∈ R n , we define the partial

order ⪯ on R n to be z 1 ⪯ z 2 if and only if z i 1 ≤ z i 2 for all i = 1, . . . , n. We use indifferently z 1 ⪯ z 2 and z 2 ⪰ z 1 . [z 1 , z 2 ] = {z ∈ R n | z 1 ⪯ z ⪯ z 2 } defines a closed interval of R n .
We define max(z 1 , z 2 ), or min(z 1 , z 2 ), to be the vector z whose components are

z i = max(z i 1 , z i 2 ), or z i = min(z i 1 , z i 2 ) respectively. Given a set Z ⊆ R n , int Z and cl Z denote the
interior and the closure of the set Z, inf Z and sup Z denote the infimum and the supremum of Z, i.e. the greatest lower and least upper bounds of Z relative to partial order ⪯ on R n . Given a collection of sets Z k ⊆ R n , for k in some index set, we will use the following

convention k∈∅ Z k = R n . For z ∈ R n and Z ⊆ R n , we define z + Z = {z + z ′ | z ′ ∈ Z}. A relation R ⊆ Z × Y is identified with the set-valued map R : Z ⇒ Y where R(z) = {y ∈ Y | (z, y) ∈ R}. We define the converse relation R -1 = {(y, z) ∈ Y × Z| (z, y) ∈ R}.
A single-valued map f : Z → Y can be identified to the deterministic set-valued map F : Z ⇒ Y where F (z) = {f (z)} for all z ∈ Z. In this paper, to avoid confusion, we will refer to set-valued maps as maps and to single-valued maps as functions. A map

F : Z ⇒ Y where Y ⊆ R m is called an interval-valued map if for all z ∈ Z, there exists y 1 , y 2 ∈ Y such that F (z) = [y 1 , y 2 ]. An over-approximation of a map F : Z ⇒ Y is a map F : Z ⇒ Y such that F (z) ⊆ F (z), for all z ∈ Z. A property is true for almost all z ∈ Z if it holds for all z ∈ Z \ Z 0 where Z 0 ⊆ Z is of measure zero. Then given two maps F 1 , F 2 : Z ⇒ Y , we say that F 1 = F 2 almost everywhere (a.e.) if F 1 (z) = F 2 (z)
for almost all z ∈ Z.

II. LEARNING MONOTONE MAPS FROM DATA

In this section, we consider the problem of learning monotone maps from data. For a class of monotone maps, we provide an approach to compute over-approximations of the maps that are as tight as possible, given the available data.

As we are dealing with set-valued maps, it is adequate to properly define the notion of monotonicity for such maps.

Definition 1: The map An illustration of the notion of monotone maps is shown in Figure 1. Let us remark that in the case of a function f : Z → Y , Definition 1 coincides with the usual definition of monotone functions, that is

F : Z ⇒ Y , with Z ⊆ R n and Y ⊆ R m , is monotone if for all z, z ′ ∈ Z with z ⪯ z ′ , ∀y ∈ F (z), ∃y ′ ∈ F (z ′ ), y ⪯ y ′ , and
∀y ′ ∈ F (z ′ ), ∃y ∈ F (z), y ⪯ y ′ .
∀z, z ′ ∈ Z, z ⪯ z ′ =⇒ f (z) ⪯ f (z ′ ).
This paper considers a class of maps given by monotone functions with additive bounded disturbances. Formally, let us consider a map

F : Z ⇒ R m where Z ⊆ R n and such that ∀z ∈ Z, F (z) = f (z) + W, (1)
where

f : Z → R m is a monotone function and W = [w, w],
W ⊆ R m is a bounded interval of disturbances. The following property of F is a straightforward consequence of (1) and of Definition 1 and is therefore stated without proof. Claim 1: The map F given by ( 1), where f is a monotone function and W ⊆ R m , is a monotone map.

We can now describe the problem under consideration in this section:

• Let us consider a map F : Z ⇒ R m of the form (1)
where the monotone function f : Z → R m is unknown and the disturbances lower and upper bounds w, w ∈ R m are known. • Let us consider a set of data D ⊆ Z × R m generated by the map F :

D = {(z k , y k ) | y k ∈ F (z k ), k ∈ K}, (2)
where K is a finite set of indices. Given the data D, the bounds w, w, and under the sole assumption that f is monotone, we aim at computing an overapproximation of the map F that is as "tight" as possible. A precise notion of tightness will be formally defined in subsection II-B.

Remark 1: In practical applications, we are mostly interested in maps F : Z ⇒ R m with Z ⊆ R n . Note that any such map can be seen as a map from a subset of R n to R m and can thus be cast in our framework. Moreover, working with extended real numbers will facilitate the mathematical analysis of our approach.

A. Consistent and simulating maps

Let us define several notions that will help us formalize our problem. A map will be said to be consistent if it is of the form (1) and is capable of generating the data D. that, for all z ∈ Z, F (z) = f (z) + W ; 2) For all (z, y) ∈ D, y ∈ F (z). We denote the set of maps that are consistent with the data D by C D .

Obviously, there is at least one map that is consistent with D, which is the map F that generated D. In general, there could be more. Then, a map over-approximating all the consistent maps is called a simulating map. Remark 2: It is important to emphasize that while minimality holds only almost everywhere, it is required that the simulation property holds for all z ∈ Z. Hence, a minimal simulating map is always an over-approximation of the unknown map F that generated the data D. The main reason for requiring minimality only almost everywhere is that it will make it possible to give a simple construction of the minimal simulating map as shown in subsection II-B.

In subsection II-B, we will provide a constructive proof of the existence of minimal simulating maps. The following proposition describes the relation between minimal simulating maps.

Proposition 1: Given two maps S m , S ′ m : Z ⇒ R m , the following properties hold: Proof: We prove the first part of the proposition. According to Definition 4, there exists a set Z 0 ⊆ Z of measure zero such that, for all z ∈ Z \ Z 0 , S m (z) ⊆ S ′ m (z), and there exists a set Z ′ 0 ⊆ Z of measure zero such that, for all z ∈ Z \Z ′ 0 , S ′ m (z) ⊆ S m (z). Hence, for all z ∈ Z \(Z 0 ∪Z ′ 0 ), S ′ m (z) = S m (z). The set Z 0 ∪ Z ′ 0 being of measure zero, it follows that S m = S ′ m a.e.. Now, let us show the second part of the proposition. For S ′ m to be a minimal simulating map, S ′ m should meet the two requirements of Definition 4. The first one is satisfied by Fig. 2. Construction of a minimal simulating map. Left: data points in the input space z k ∈ Z, k ∈ K, points larger than (k ∈ K + (z)) / lesser than (k ∈ K -(z)) / incomparable to a given z ∈ Z are represented with squares / circles / triangles, respectively. Right: data points in the output space y k ∈ R m , k ∈ K are surrounded by red intervals representing the bounded disturbance, and the green interval represents the minimal simulating map Sm(z) ⊆ R m given by (3).

• If S m
The light green area represents the approximation of the function f as mentioned in Remark 3.

assumption. The second one can be shown as follows. There exists a set Z 0 ⊆ Z of measure zero such that, for all z ∈ Z \ Z 0 , for all S ∈ S D , S m (z) ⊆ S(z). Also, there exists a set

Z ′ 0 ⊆ Z of measure zero such that, for all z ∈ Z \ Z ′ 0 , S ′ m (z) = S m (z). Hence, for all z ∈ Z \ (Z 0 ∪ Z ′ 0 ), for all S ∈ S D , S ′ m (z) ⊆ S(z). Because Z 0 ∪ Z ′
0 is of measure zero, the second requirement of Definition 4 is satisfied, and S ′ m is a minimal simulating map.

Proposition 1 essentially states that minimal simulating maps are uniquely determined up to a set of measure zero. We can now formally state the problem under consideration in this section:

Problem 1: Given the data D and the disturbance bounds w, w ∈ R m , compute S m : Z ⇒ R m , a minimal simulating map of D.

B. Computation of minimal simulating maps

In the following, we describe a solution to Problem 1. We first establish a characterization of a minimal simulating map and prove some of its properties. Then, a practical approach is provided for computing a minimal simulating map.

1) Characterization and properties:

The following theorem provides an effective characterization of a minimal simulating map:

Theorem 1: Let S m : Z ⇒ R m be the map given for all z ∈ Z by:

S m (z) =   k∈K -(z) y ∈ R m y k + w -w ⪯ y   ∩   k∈K + (z) y ∈ R m y ⪯ y k + w -w   , (3)
where

K -(z) = {k ∈ K | z k ⪯ z}, K + (z) = {k ∈ K | z ⪯ z k }.
(4) Then, S m is a minimal simulating map of D.

Figure 2 shows the construction of the map S m . Let us remark that Theorem 1 also provides a constructive proof of the existence of minimal simulating maps. It follows from (3) that the minimal simulating map S m is interval-valued. Indeed, for all z ∈ Z, S m (z) = [S m (z), S m (z)] with upper and lower bounds S m (z), S m (z) given by

S m (z) = inf y k + w -w | k ∈ K + (z) , (5)
S m (z) = sup y k + w -w | k ∈ K -(z) . (6)
We now prove Theorem 1.

Proof: Let us first show that S m ∈ S D . Let F ∈ C D , then there exists a monotone function f :

Z → R m such that, F (z) = f (z) + W . Let z ∈ Z, y ∈ F (z), then, y ⪯ f (z) + w. From the monotonicity of f , we have f (z) ⪯ f (z k ) for all k ∈ K + (z). Therefore, y ⪯ f (z k ) + w, ∀k ∈ K + (z). Moreover, F ∈ C D implies y k ∈ F (z k ), for all k ∈ K. Therefore, we have f (z k ) + w ⪯ y k , ∀k ∈ K.
Hence, we get from the inequalities above that

y ⪯ y k + w -w, ∀k ∈ K + (z).
Then, from ( 5), y ⪯ S m (z). Similarly, it can be shown that y ⪰ S m (z). Therefore, we have y ∈ S m (z); hence, S m is a simulating map of D. Now, we prove minimality. For any arbitrary point included in S m we will build a consistent map containing this point. Consider the set Z 0 of measure zero, defined as

Z 0 = {z ∈ Z | ∃i ∈ {1, . . . , n}, k ∈ K, z i = z i k }. (7)
This set consists of the union of finitely many hyperplanes. Each one of those hyperplanes is defined by one of the components of a data point. According to Definition 4, minimality has to hold almost everywhere. Then, let us consider z ⋆ ∈ Z \ Z 0 and y ⋆ ∈ S m (z ⋆ ). Let us define a partition of Z in 3 regions defined by

Z 1 = Z + (z ⋆ ) \ Z -(z ⋆ ), Z 2 = Z + (z ⋆ ) ∩ Z -(z ⋆ ), Z 3 = Z \ (Z 1 ∪ Z 2 ),
where

Z -(z ⋆ ) = n i=1 z ∈ Z | z i > sup{z i k | z i k ≤ z ⋆i , k ∈ K} , Z + (z ⋆ ) = n i=1 z ∈ Z | z i < inf{z i k | z i k ≥ z ⋆i , k ∈ K} .
A representation of the regions can be found in Figure 3. Let us remark that by construction

z k ∈ Z 1 ∪ Z 3 for all k ∈ K and that z ⋆ ∈ Z 2 because z ⋆ ∈ Z \ Z 0 .
Then, let us consider the map The area bounded from below by the orange line represents the set Z -(z ⋆ ). The area Z 1 in blue contains all the data points smaller than z ⋆ . The area Z 3 in red contains all the data points larger than z ⋆ or incomparable with z ⋆ . The area Z 2 is colored in green.

F ⋆ : Z ⇒ R m given for all z ∈ Z by F ⋆ (z) = f ⋆ (z) + W where f ⋆ : Z → R m is defined as follows f ⋆ (z) =   
In Figure 4, a representation, in the one-dimensional case, of the maps F, F ⋆ , and S m can be found. Since

z ⋆ ∈ Z 2 and since y ⋆ ∈ S m (z ⋆ ) = [S m (z ⋆ ), S m (z ⋆ )], it follows from (8) that y ⋆ ∈ F ⋆ (z ⋆ ). We will now prove that F ⋆ ∈ C D .
We first prove that function f ⋆ is monotone. Let z, z ′ ∈ Z such that z ⪯ z ′ . There are several possibilities for the position of z and z

′ : i) z, z ′ ∈ Z 1 , ii) z, z ′ ∈ Z 2 , iii) z, z ′ ∈ Z 3 , iv) z ∈ Z 1 , z ′ ∈ Z 2 , v) z ∈ Z 2 , z ′ ∈ Z 3 , vi) z ∈ Z 1 , z ′ ∈ Z 3 . We want to verify that f ⋆ (z) ⪯ f ⋆ (z ′ ) in all those cases. • Case (i): Since z ⪯ z ′ , we have K -(z) ⊆ K -(z ′ ), which implies from (6), S m (z) ⪯ S m (z ′ ). Then, from (8), f ⋆ (z) ⪯ f ⋆ (z ′ ). • Case (ii): From (8), we have f ⋆ (z) = f ⋆ (z ′ ). • Case (iii): Similar to case (i). • Case (iv): It follows from the definition of Z 2 that K -(z ′ ) = K -(z ⋆ ). Then, since z ⪯ z ′ , we have K -(z) ⊆ K -(z ′ ) and thus K -(z) ⊆ K -(z ⋆ ). Hence, from (6), we get S m (z) ⪯ S m (z ⋆ ). Since z ∈ Z 1 and z ′ ∈ Z 2 we get from (8), f ⋆ (z) = S m (z) -w ⪯ S m (z ⋆ ) -w ⪯ f ⋆ (z ′ ). • Case (v): Similar to case (iv). • Case (vi): Since S m ∈ S D and F ∈ C D , we get that S m (z) -w ⪯ f (z) and f (z ′ ) ⪯ S m (z ′ ) -w. Since f is monotone and z ⪯ z ′ , we get f (z) ⪯ f (z ′ ). Finally, since z ∈ Z 1 and z ′ ∈ Z 3 , we get from (8), f ⋆ (z) = S m (z) -w ⪯ S m (z ′ ) -w = f ⋆ (z ′ ).
Hence, we established the monotonicity of function f ⋆ . We now prove that the second requirement of Definition 2 is satisfied. Let us consider k ∈ K, we already know that z k ∈ Z 1 ∪ Z 3 . In the case where z k ∈ Z 1 , we have from ( 6) and [START_REF] Girard | Least-violating symbolic controller synthesis for safety, reachability and attractivity specifications[END_REF],

y k ⪯ S m (z k ) -w + w = f ⋆ (z k ) + w.
Moreover, since S m ∈ S D , we have

y k ⪰ S m (z k ) = f ⋆ (z k ) + w. Therefore, y k ∈ F ⋆ (x k ). A similar statement can be shown when z k ∈ Z 3 . It follows that F ⋆ ∈ C D .
Disturbance Fig. 4. Construction of the map F ⋆ according to [START_REF] Girard | Least-violating symbolic controller synthesis for safety, reachability and attractivity specifications[END_REF].

Hence, we have proved that for all z ⋆ ∈ Z \ Z 0 and for all y ⋆ ∈ S m (z ⋆ ), there exists a map

F ⋆ ∈ C D such that y ⋆ ∈ F ⋆ (z ⋆ ).
Hence, for all map S ∈ S D , we have y ⋆ ∈ S(z ⋆ ). Hence, S m is a minimal simulating map.

We now show some useful properties of minimal simulating maps. The first one states that the minimal simulating map in (3) inherits the monotonicity property:

Proposition 2: The minimal simulating map S m defined in (3) is monotone. Proof: Let z, z ′ ∈ Z such that z ⪯ z ′ . Then, K + (z ′ ) ⊆ K + (z) and K -(z) ⊆ K -(z ′
) which give from ( 5) and ( 6), S m (z) ⪯ S m (z ′ ) and S m (z) ⪯ S m (z ′ ). Then, for all y ∈ S m (z), y ⪯ S m (z) ⪯ y ′ with y ′ = S m (z ′ ) ∈ S m (z ′ ). Similarly, we can show that for all y ′ ∈ S m (z ′ ), y ′ ⪰ S m (z ′ ) ⪰ y with y = S m (z) ∈ S m (z). Therefore, according to Definition 1, S m is monotone.

The second result characterizes the minimal simulating map obtained from the fusion of two data sets: Proposition 3: Given two data sets D ′ , D ′′ , consider the data set D = D ′ ∪ D ′′ , and let S m , S ′ m , S ′′ m be the minimal simulating maps of D, D ′ , D ′′ given by (3). Then, we have

∀z ∈ Z, S m (z) = S ′ m (z) ∩ S ′′ m (z). Proof:
The data set D can be written as

D = {(z k , y k ) | y k ∈ F (z k ), k ∈ K D }, with K D = K D ′ ∪ K D ′′
where K D ′ and K D ′′ are the sets of indices of D ′ and D ′′ . Then, for all z ∈ Z, we have

K - D (z) = K - D ′ (z) ∪ K - D ′′ (z), K + D (z) = K + D ′ (z) ∪ K + D ′′ (z). Substituting in (3), we can easily show that for all z ∈ Z that S m (z) = S ′ m (z) ∩ S ′′ m (z)
. Proposition 3 has several implications. The first one is that the computation of a minimal simulating map for large data sets can be easily parallelized. Indeed, one can split the data into multiple data sets, compute minimal simulating maps for these data sets concurrently, and then compute their intersection. The second implication is that pre-computed minimal simulating maps can be easily refined when new Fig. 5. The figure shows the rectangular partition of the input set Z based on the collected data points and an example of one of the discrete variables q. All the points in yellow, aggregated in q, have the same output.

data becomes available. Indeed, it is sufficient to compute a minimal simulating map for the new data and to intersect with the pre-computed one.

Remark 3: In this section, we provided a solution to Problem 1 for computing a minimal simulating map of D, which provides a tight over-approximation of the map F . Note that similar results can be obtained if we seek an overapproximation of the function f as it can be shown that

∀z ∈ Z, f (z) ∈ [S m (z) -w, S m (z) -w],
which can be seen in Figure 2. If we want to calculate a singlevalued approximation of the function f , we can use a similar approach to the one used in set membership technique [START_REF] Novara | Set membership identification of nonlinear systems[END_REF] and build an estimation of f using the calculated bounds on the minimal simulating map

f (z) = S m (z) -w + S m (z) -w 2 .
It will be apparent from what follows that this estimation is piecewise constant and not continuous. Continuous estimation of f included in the over-approximation map S m is derived in [START_REF] Makdesi | Safe learning-based model predictive control using the compatible models approach[END_REF] using the class of multi-affine functions.

2) Computation: After characterizing minimal simulating maps, we move to address how to compute one in practice. Our approach and the used notations are illustrated in Figure 5, where a six data points case is demonstrated.

For that purpose, we introduce a rectangular partition of the input set Z, induced by the data D. For simplicity, let us

assume that Z is a closed interval of R n , i.e. Z = [α, α].
This is without loss of generality since it is always possible to embed Z in such a set. For each i ∈ {1, . . . , n}, we sort the i th components of all the data points in the input space (i.e., z i k , for k ∈ K), so we have z i

k 1 i ≤ • • • ≤ z i k |K| i
. The sorted values are then used to define the finite partitions

(Z i q i ) q i ∈Q i of [α i , α i ] where Q i = {0, . . . , |K|} and        Z i 0 = [α i , z i k 1 i ), Z i q i = [z i k q i i , z i ki q i +1 ), q i = 1, . . . , |K| -1, Z i |K| = [z i k |K| i , α i ].
Then, let us define

Q = Q 1 × • • • × Q n ,
and let the finite rectangular partition (Z q ) q∈Q of Z be given for q = (q 1 , . . . , q n ) by The function ϕ maps all the points inside a cell (e.g. square in yellow) to a discrete variable. The map σ is the minimal interval-valued simulating map.

Z q = Z 1 q 1 × • • • × Z n q n .
Lemma 1: Let S m be the minimal simulating map of D given by (3). Then, for all q ∈ Q, for all z, z ′ ∈ int Z q , S m (z) = S m (z ′ ).

Proof: Because the partition is defined using the components of data points, we have, for all q ∈ Q, for all

z, z ′ ∈ int Z q , K -(z) = K -(z ′ ) and K + (z) = K + (z ′ ). Subsequently, we obtain by (3), S m (z) = S m (z ′ ).
Building on the property presented in Lemma 1, we can introduce the new map σ m : Q ⇒ R m defined as follows

∀q ∈ Q, σ m (q) = S m (z), for z ∈ int Z q . ( 9
)
We also define a quantization function ϕ m : Z → Q associated to the finite data-induced partition (Z q ) q∈Q as follows

∀z ∈ Z, ∀q ∈ Q, ϕ m (z) = q ⇐⇒ z ∈ Z q . ( 10
)
Proposition 4: Let σ m and ϕ m be given by ( 9) and ( 10), then σ m • ϕ m is a minimal simulating map of D.

Proof: It follows directly from Lemma 1 and from the definition of σ m and ϕ m that σ m • ϕ m = S m a.e.. Then, let us show that σ m • ϕ m is a simulation map. From the construction of the partition (Z q ) q∈Q , it follows that for all q ∈ Q, for all z ∈ Z q , for all z

′ ∈ int Z q , K -(z ′ ) ⊆ K -(z) and K + (z ′ ) ⊆ K + (z). Hence, from (3), S m (z) ⊆ S m (z ′ ). Then, it results that S m (z) ⊆ σ m • ϕ m (z), for all z ∈ Z. Since S m ∈ S D , we have that σ m •ϕ m ∈ S D . Therefore, according to Proposition 1, σ m • ϕ m is a minimal simulating map of D.
Hence, computing the finite partition (Z q ) q∈Q and the map σ m offers an effective way to compute a minimal simulating map and store it. Let us remark that the number of elements in Q is (|K|+1) n . Hence, it grows polynomially with the number of data points |K| and exponentially with the dimension n of the input space Z. Therefore, it is computationally prohibitive to use this approach for large data sets, particularly with highdimensional input spaces. In the following, we tackle this problem by fixing a partition a priori and by finding the best simulating map on this partition.

C. Simulating maps on fixed partitions

Instead of relying on the data set D to partition Z, which makes calculating the simulating map computationally expensive, we will assume that a rectangular partition of Z is given a priori. In the following, we characterize the tightest interval-valued simulating map on this partition. Our approach is illustrated in Figure 6 1) Construction: For each coordinate i ∈ {1, . . . , n}, let be given finite partitions

(D i q i ) q i ∈Q i of [α i , α i ] where Q i = {0, . . . , K i } and    D i 0 = [α i , α i 1 ), D i q i = [α i q i , α i q i +1 ), q i = 1, . . . , K i -1, D i K i = [α i K i , α i ], where α i < α i 1 < • • • < α i K i < α i .
Note that the number of partition elements K i +1 can be different from one coordinate to another. Then, let us define

Q = Q 1 × • • • × Q n ,
and let the finite rectangular partition (D q ) q∈Q of Z be given for q = (q 1 , . . . , q n ) by

D q = D 1 q 1 × • • • × D n q n . Let us define the map σ : Q ⇒ R m given for all q ∈ Q by σ(q) =   k∈K -(z q ) {y ∈ R m | y k + w -w ⪯ y}   ∩   k∈K + (zq) {y ∈ R m | y ⪯ y k + w -w}   (11)
where z q = inf D q , z q = sup D q , and K -, K + are defined as in (4). From ( 5), ( 6) and ( 11), we get that σ is an intervalvalued map: for all q ∈ Q, σ(q) = [σ(q), σ(q)] with σ(q) = S m (z q ) and σ(q) = S m (z q ). We also consider a quantization function ϕ : Z → Q associated to the finite partition (D q ) q∈Q and defined as

∀z ∈ Z, ∀q ∈ Q, ϕ(z) = q ⇐⇒ z ∈ D q . ( 12
)
Theorem 2: Let σ and ϕ be given by ( 11) and ( 12), then the following properties hold:

1) σ • ϕ ∈ S D (Simulation); 2) For any interval-valued map σ ′ : Q ⇒ R m , such that σ ′ • ϕ ∈ S D , it holds for all q ∈ Q, σ(q) ⊆ σ ′ (q) (Minimality).

Proof: Let S m be the minimal simulating map in (3). Let us prove the simulation property. Let z ∈ Z, then z ϕ(z) ⪯ z ⪯ z ϕ(z) , which yields

K -(z ϕ(z) ) ⊆ K -(z) and K + (z ϕ(z) ) ⊆ K + (z).
Hence, by ( 3) and [START_REF] Hewing | Learningbased model predictive control: Toward safe learning in control[END_REF]

, S m (z) ⊆ σ • ϕ(z). Since S m ∈ S D , we get σ • ϕ ∈ S D .
Next, let us prove the minimality property. Let us consider an interval-valued map σ

′ : Q ⇒ R m such that σ ′ • ϕ ∈ S D ,
and let q ∈ Q. From (4), it can be seen that there exists a neighborhood N q of z q such that for all z ∈ N q with z ⪯ z q it holds K + (z) = K + (z q ), which yields by (3) S m (z) = S m (z q ). Then, let z ⋆ ∈ (N q ∩ D q ) \ Z 0 with Z 0 the set of measure zero defined as in [START_REF] Girard | Controller synthesis for safety and reachability via approximate bisimulation[END_REF]. Let us remark that z ⋆ ⪯ z q and hence S m (z ⋆ ) = S m (z q ). From the proof of Theorem 1, we get that there exists

F ⋆ ∈ C D such that S m (z q ) = S m (z ⋆ ) ∈ F ⋆ (z ⋆ ). Then, σ ′ • ϕ ∈ S D gives us that S m (z q ) ∈ σ ′ • ϕ(z ⋆ ). Moreover, z ⋆ ∈ D q , gives us that σ ′ • ϕ(z ⋆ ) = σ ′ (q). Hence, S m (z q ) ∈ σ ′ (q).
Similarly, we can show that S m (z q ) ∈ σ ′ (q). Then since σ ′ is interval-valued, we get that [S m (z q ), S m (z q )] ⊆ σ ′ (q). Then, σ(q) = [S m (z q ), S m (z q )] gives us σ(q) ⊆ σ ′ (q).

It follows from Theorem 2 that it is possible to define a notion of minimal simulating map of D relative to a given partition (D q ) q∈Q of Z. It should be noticed that similar results to Propositions 2 and 3 hold for the simulating map σ • ϕ. Finally, one can check that if the partition (D q ) q∈Q of Z coincides with the data-induced partition defined in the previous section, σ • ϕ is minimal simulating map of D. The difference between the maps σ •ϕ and σ m •ϕ can be examined by comparing Figure 5 and Figure 6.

Remark 4: For a given number of cells, typically smaller than data-induced partition, choosing an optimal partition (D q ) q∈Q is a complicated problem. However, it can be easily shown that the optimal partition (achieving the minimal overapproximation volume) is aligned with some of the data point components. This makes it possible to rely on heuristics to choose good yet sub-optimal partitions. For instance, a possible approach would be to aggregate cells of the datainduced partition where a lot of data is available.

2) Efficient computation: Given the partition (D q ) q∈Q of Z, computing the simulating map defined in the previous section amounts to computing the map σ : Q ⇒ R m given by [START_REF] Hewing | Learningbased model predictive control: Toward safe learning in control[END_REF].

A straightforward algorithm to compute σ is as follows. For each q ∈ Q, we go through all the data points in D and determine K -(z q ) and K + (z q ). Then, one gets σ(q) by ( 11). That makes the overall complexity O(|K| × |Q|), i.e. bilinear with respect to the number of data points and to the number of the partition elements. In this section, we present a more efficient approach to calculating σ.

For simplicity, we assume in the following that for all k ∈ K, z k ∈ int Z, i.e. that no data point lies on the boundary of the input set. Then, since σ(q) = S m (z q ), it follows from (5) that for all q ∈ Q such that q i = K i , for some i ∈ {1, . . . , n}, σ(q) = (+∞, . . . , +∞). Similarly, since σ(q) = S m (z q ), it follows from (6) that for all for all q ∈ Q such that q i = 0, for some i ∈ {1, . . . , n}, σ(q) = (-∞, . . . , -∞).

Then, let us consider the following subsets of Q:

Q = {0, . . . , K 1 -1} × • • • × {0, . . . , K n -1}, Q = {1, . . . , K 1 } × • • • × {1, . . . , K n }.
One needs to compute σ(q) and σ(q) for q ∈ Q and q ∈ Q, respectively. For that purpose, let us first define the functions

σ 0 : Q → R m and σ 0 : Q → R m as follows: σ 0 (q) = inf y k + w -w | z k ∈ cl D up(q) , (13)
σ 0 (q) = sup y k + w -w | z k ∈ cl D lo(q) , (14)
where up(q) = q + 1 n , lo(q) = q -1 n , and

1 n = (1, . . . , 1
). An illustration of the elements up(q) and lo(q) can be seen in Figure 7. To compute σ 0 , σ 0 , we start by initializing σ 0 (q) = (+∞, . . . , +∞) for all q ∈ Q, and σ 0 (q) = (-∞, . . . , -∞) for all q ∈ Q. Then, we go through all the points in the set D; for each entry (z k , y k ) we find all q such that z k ∈ cl D q . Then, we update the value of σ 0 (up(q)) and σ 0 (lo(q)) using ( 13) and ( 14). The partition is stored and sorted component-wise, so to find q, we can do a binary search for each component of z k . Therefore, computing σ 0 , σ 0 is done with a complexity O(|K| × i log(K i )) or equivalently O(|K| × log(|Q|)).

q u ad(q) u ad(q) l ad(q) l ad(q) lo(q) up(q) σ σ Fig. 7. A representation of the sets up(q), u ad(q), lo(q), l ad(q) and the order in which we calculate σ and σ.

We now present a result that will allow us to compute the map σ sequentially.

Proposition 5: Let σ be the interval-valued map given by [START_REF] Hewing | Learningbased model predictive control: Toward safe learning in control[END_REF], its upper and lower bounds σ, σ satisfy:

∀q ∈ Q, σ(q) = min (inf {σ(q ′ )| q ′ ∈ u ad(q)} , σ 0 (q)) , ( 15
) ∀q ∈ Q, σ(q) = max (sup {σ(q ′ )| q ′ ∈ l ad(q)} , σ 0 (q)) , (16)
where u ad(q) = {q ′ ∈ Q| ∃k ∈ {1, . . . , n}, q ′ -q = e k }, l ad(q) = {q ′ ∈ Q| ∃k ∈ {1, . . . , n}, qq ′ = e k }, and e k ∈ R n whose k th component is 1 and all others are 0.

An illustration of the sets u ad(q) and l ad(q) can be seen in Figure 7.

Proof: We prove the property for the upper bound σ; the proof for the lower bound σ follows similarly. Let q ∈ Q, it can be seen that

K + (z q ) = q ′ ∈u ad(q) K + (z q ′ ) ∪ {k ∈ K| z k ∈ D up (q)} .
Then, since σ(q) = S m (z q ), (15) follows directly from (5) and the equality above.

From Proposition 5, we can see that to compute σ we go through all q ∈ Q sequentially in a decreasing order, starting from q = (K 1 , . . . , K n ) as represented by Figure 7. For σ we start from q = (1, . . . , 1) and though all q ∈ Q in an increasing order. Proof: We already showed the complexity of computing σ 0 and σ 0 is O(|K| × log(|Q|)). To compute σ we should go through all the elements q ∈ Q twice, one in decreasing order to compute σ and one in increasing order to compute σ. Therefore, the complexity of computing σ is O(|K| × log(|Q|) + |Q|).

III. DATA-DRIVEN ABSTRACTION

In this section, we show that our approach presented in the previous section can be used for data-driven modeling of discrete-time dynamical systems. We show that simulating maps computed from data can be used to define sound finite state abstractions for dynamical systems with unknown monotone dynamics. Then, we extend our approach to nonmonotone systems with bounded derivatives.

We consider a discrete-time dynamical system:

x(τ + 1) ∈ f (x(τ ), u(τ )) + g(x(τ ), u(τ )) + W (17)
where x(τ ) ∈ X ⊆ R nx and u(τ ) ∈ U ⊆ R nu denote the state and the control input, g :

X × U → R nx is a known function whereas f : X × U → R
nx is an unknown monotone function with respect to z = (x, u), W = [w, w] ⊆ R nx is a bounded interval of disturbances with known bounds. Let be given a set of data D ⊆ X × U × R nx consisting of transitions sampled from the system (17):

D = (x k , u k , x ′ k ) x ′ k ∈ f (x k , u k ) + g(x k , u k ) + W, k ∈ K (18)
where K is a finite set of indices. Because g is known, we can construct an auxiliary data set D ′ as follows:

D ′ = {(x k , u k , y k )| y k = x ′ k -g(x k , u k ), k ∈ K}.
The new auxiliary data set can be seen as a set of data generated by a monotone map of the form [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF].

Therefore, we can use the approach presented in the previous section to compute a simulating map S : X × U ⇒ R nx of data D ′ . Then, a data-driven model of system [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF] can be defined as follows x(τ + 1) ∈ S(x(τ ), u(τ )) + g(x(τ ), u(τ )).

(

) 19
In the following, we formally relate the behaviors of ( 17) and [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF].

A. Alternating simulation relation

Alternating simulation [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF] is a formal relationship between the behaviors of two systems that makes it possible to refine a controller synthesized for one system in order to control the other while preserving guarantees of correctness. Alternating simulation relations are usually defined within the framework of transition systems.

Definition 5: A transition system T is a tuple T = (X, U, ∆, Y, H), where X is a set of states, U is a set of inputs, ∆ : X × U ⇒ X is a transition relation, Y is a set of outputs, and H : X → Y is an output map.

An input u ∈ U is called enabled at x ∈ X if ∆(x, u) ̸ = ∅.
The set of all inputs enabled at x is denoted enab ∆ (x).

We define transition systems T sys = (X, U, ∆ sys , Y, H) and T data = (X, U, ∆ sys , Y, H) associated to [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF] and [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] where the set of states X and inputs U are the same as in [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF] and [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. The transition relation ∆ sys is defined as follows, for all x ∈ X: [START_REF] Makdesi | Data-driven abstraction of monotone systems[END_REF] and ∀u ∈ enab ∆sys (x), ∆ sys (x, u) = f (x, u) + g(x, u) + W.

u ∈ enab ∆sys (x) ⇐⇒ f (x, u) + g(x, u) + W ⊆ X,
(21) Similarly, the transition relation ∆ data is defined as follows, for all x ∈ X:

u ∈ enab ∆ data (x) ⇐⇒ S(x, u) + g(x, u) ⊆ X, (22)
and

∀u ∈ enab ∆ data (x), ∆ data (x, u) = S(x, u) + g(x, u). ( 23
)
Let us remark that ( 20) and [START_REF] Makdesi | Safe learning-based model predictive control using the compatible models approach[END_REF] ensure that an input u is enabled at state x only if it is guaranteed that the next state of ( 17) and ( 19) belongs to the set of states X. The set of outputs Y and the output map H are left unspecified. They can be chosen arbitrarily but are assumed to be the same for T sys and T data . One can, for instance, choose Y = X and H to be the identity map.

In order to relate formally the behaviors of T sys and T data , we recall the notion of alternating simulation relation [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF]:

Definition 6: Let us consider two transition systems

T i = (X i , U i , ∆ i , Y i , H i ) i = 1, 2, sharing the same sets of outputs (Y 1 = Y 2 = Y ). A relation R ⊆ X 1 × X 2 is an alternating simulation relation from T 1 to T 2 if the following conditions are satisfied: 1) for all x 1 ∈ X 1 , there exists x 2 ∈ X 2 such that (x 1 , x 2 ) ∈ R; 2) for all (x 1 , x 2 ) ∈ R, H 1 (x 1 ) = H 2 (x 2 ); 3) for all (x 1 , x 2 ) ∈ R, for all u 1 ∈ enab ∆1 (x 1 ), there ex- ists u 2 ∈ enab ∆2 (x 2 ) such that for all x ′ 2 ∈ ∆ 2 (x 2 , u 2 ), there exists x ′ 1 ∈ ∆ 1 (x 1 , u 1 ) satisfying (x ′ 1 , x ′ 2 ) ∈ R.
Then, we say that:

• T 1 is alternatingly simulated by T 2 , denoted T 1 ⪯ AS T 2 ,
if there exists an alternating simulation relation R from

T 1 to T 2 ; • T 1 is alternatingly bisimilar to T 2 , denoted T 1 ∼ =AS T 2 ,
if there exists a relation R such that R is an alternating simulation relation from T 1 to T 2 and R -1 is an alternating simulation relation from T 2 to T 1 . Proposition 7: Let S ∈ S D ′ , then for any choice of set of outputs Y and of output map H : X → Y , T data ⪯ AS T sys .

Proof: Let us show that

R = {(x data , x sys ) ∈ X × X | x data = x sys }
is an alternating simulation relation from T data to T sys . The first two conditions of alternating simulation follow directly from the form of R and from the fact that T sys and T data have the same sets of states X and of outputs Y , and the same output maps H. Since S ∈ S D ′ , it follows from Definition 3 that

∀x ∈ X, u ∈ U, f (x, u) + W ⊆ S(x, u).
Hence, S(x, u) + g(x, u) ⊆ X implies f (x, u) + g(x, u) + W ⊆ X . Therefore, from ( 20) and ( 22), we have for all x ∈ X, enab ∆ data (x) ⊆ enab ∆sys (x). Moreover, from ( 21) and ( 23), for all u ∈ enab ∆ data (x), it holds ∆ sys (x, u) ⊆ ∆ data (x, u).

Let us now show the third condition of alternating simulation. Let us consider (x sys , x data ) ∈ R , then x sys = x data . Let u data ∈ enab ∆ data (x data ), then for u sys = u data , we have u sys ∈ enab ∆sys (x sys ). Moreover, ∆ sys (x sys , u sys ) ⊆ ∆ data (x data , u data ). Therefore, for all x ′ sys ∈ ∆ sys (x sys , u sys ), there exists x ′ data ∈ ∆ data (x data , u data ) satisfying x ′ sys = x ′ data , and hence (x ′ sys , x ′ data ) ∈ R. From the previous result, it follows that the data-driven model T data can be used to synthesize controllers that can be refined into controllers for the partially unknown system T sys , with formal guarantees of correctness.

B. Symbolic abstraction

We now go one step further by computing symbolic abstractions of T data . This will allow us to use discrete controller synthesis techniques to control the system with formal guarantees on the closed-loop behavior. Let us assume that the sets of states and inputs X and U are closed intervals of R nx and R nu and let (X q ) q∈Q , (U p ) p∈P be given rectangular partitions of X and U as defined for Z in Section II-C. We define a quantization function ϕ : X × U → Q × P associated to these finite partitions as follows,

∀(x, u) ∈ X × U, ∀(q, p) ∈ Q × P, ϕ(x, u) = (q, p) ⇐⇒ x ∈ X q , u ∈ U p . (24)
Then, we can use the approach presented in Section II-C to compute a map σ :

Q×P ⇒ R nx such that σ•ϕ is a simulating map of D ′ .
Let us assume that for all (q, p) ∈ Q × P , we can compute subsets G(q, p) ⊆ R nx such that

∀x ∈ X q , ∀ u ∈ U p , g(x, u) ∈ G(q, p). ( 25
)
Such sets can be computed, for instance, using interval analysis [START_REF] Jaulin | Applied interval analysis[END_REF] or using approaches based on mixed monotonicity or on growth bounds [START_REF] Meyer | Interval reachability analysis: Bounding trajectories of uncertain systems with boxes for control and verification[END_REF].

We define a symbolic transition system T symb = (Q, P, ∆ symb , Y, H symb ) where the set of states and inputs are given by the partitions index sets Q and P , and the transition relation ∆ symb is defined as follows, for all q ∈ Q:

p ∈ enab ∆ symb (q) ⇐⇒ σ(q, p) + G(q, p) ⊆ X, (26)
and

∀p ∈ enab ∆ symb (q), ∆ symb (q, p) = q ′ ∈ Q| (σ(q, p) + G(q, p)) ∩ X ′ q ̸ = ∅ . ( 27
)
We define the set of outputs to be Y = Q and the output map H symb to be the identity map.

Theorem 3: Let S = σ • ϕ ∈ S D ′ , let T sys and T data be defined as Section III-A for the set of outputs Y = Q and the output map H : X → Q, given by ∀x ∈ X, ∀q ∈ Q, H(x) = q ⇐⇒ x ∈ X q .

(28)

Then, the following relation holds:

1) T symb ⪯ AS T data ⪯ AS T sys .
2) If g = 0, then T symb ∼ =AS T data .

Proof: The fact that T data ⪯ AS T sys follows directly from Proposition 7. Then, let us show that

R = {(q, x) ∈ Q × X | x ∈ X q }
is an alternating simulation relation from T symb to T data . Let q ∈ Q and let x ∈ X q , then (q, x) ∈ R and the first condition of alternating simulation holds. Let (q, x) ∈ R, then x ∈ X q , which by [START_REF] Pola | Control of cyber-physical-systems with logic specifications: a formal methods approach[END_REF] gives H(x) = q. Since H symb (q) = q, the second condition of alternating simulation holds.

Let us now show the third condition of alternating simulation. Let (q, x) ∈ R and p ∈ enab ∆ symb (q), then choose u ∈ U p . Since S = σ • ϕ, we have S(x, u) = σ(q, p), and by [START_REF] Milani Fard | Fast and flexible monotonic functions with ensembles of lattices[END_REF] we have g(x, u) ∈ G(q, p). Hence,

S(x, u) + g(x, u) ⊆ σ(q, p) + G(q, p). ( 29
)
Then, it follows from ( 22) and ( 26) that u ∈ enab ∆ data (x). Let x ′ ∈ ∆ data (x, u), from ( 23) and ( 29), we get x ′ ∈ σ(q, p) + G(q, p). Let q ′ ∈ Q such that x ′ ∈ X q ′ , from ( 27), we get that q ′ ∈ ∆ symb (q, p). Moreover, since x ′ ∈ X q ′ , (q ′ , x ′ ) ∈ R.

Hence, we proved that T symb ⪯ AS T data . Now, let us assume that for all x ∈ X, and u ∈ U , g(x, u) = 0. Then, [START_REF] Milani Fard | Fast and flexible monotonic functions with ensembles of lattices[END_REF] holds with G(q, p) = {0}, for all q ∈ Q and p ∈ P . Let us show that in that case, R -1 is an alternating simulation relation from T data to T symb . Let x ∈ X, since (X q ) q∈Q is a partition of X, there exists q ∈ Q such that x ∈ X q . Hence, (x, q) ∈ R -1 and the first condition of alternating simulation holds. The proof that the second condition holds for R -1 is the same as for R.

Then, let us show the third condition of alternating simulation. Let (x, q) ∈ R -1 and u ∈ enab ∆ data (x), then there exists p such that u ∈ U p . Since S = σ•ϕ, we have S(x, u) = σ(q, p). Moreover, since g(x, u) = 0 and G(q, p) = {0}, we get from ( 22) and ( 26) that p ∈ enab ∆ symb (q). Let q ′ ∈ ∆ symb (q, p), from [START_REF] Novara | Set membership identification of nonlinear systems[END_REF] we get that σ(q, p) ∩ X q ′ ̸ = ∅. Hence, S(x, u)∩X q ′ ̸ = ∅. Then, let x ′ ∈ S(x, u)∩X q ′ . From (23), x ′ ∈ ∆ data (x, u). Moreover (x ′ , q ′ ) ∈ R -1 . Hence, we proved that T data ∼ =AS T symb .

Theorem 3 shows that the symbolic abstraction T symb can be used to synthesize controllers for the data-driven system T data and hence also for the partially unknown system T sys . Moreover, when T sys is fully unknown and monotone (i.e. when g = 0), then Theorem 3 shows that working with the symbolic abstraction T symb does not bring any conservatism compared to working with the data-driven model T data .

Remark 5: The fact that T symb has only a finite number of state and input values allows us to use algorithmic techniques to synthesize controllers for various specifications such as safety, reachability, or attractivity [START_REF] Girard | Least-violating symbolic controller synthesis for safety, reachability and attractivity specifications[END_REF] or even more complex specifications expressed e.g. in linear temporal logic [START_REF] Belta | Formal methods for discrete-time dynamical systems[END_REF] or using hybrid automata [START_REF] Sinyakov | Formal controller synthesis from specifications given by discrete-time hybrid automata[END_REF]. Due to the fact that there is an alternating simulation relation between T sys and T symb , the synthesized controller can be refined such that the desired closed-loop behavior of T sys is guaranteed to be correct-bydesign (see e.g. [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF] for more details on discrete controller refinement).

C. Systems with bounded derivatives

We now show how our approach can be adapted to compute data-driven models for systems of the form [START_REF] Khajenejad | Interval observers for simultaneous state and model estimation of partially known nonlinear systems[END_REF] where the unknown dynamics f is not necessarily monotone but has bounded derivatives with known upper and lower bounds. Our construction is inspired by the approach presented in [START_REF] Yang | On sufficient conditions for mixed monotonicity[END_REF] for computing decomposition functions of mixed-monotone functions.

Let us assume that for all x ∈ X, u ∈ U :

∂f i ∂x j (x, u) ∈ [a ij , a ij ], i, j ∈ {1, . . . , n x }, ∂f i ∂u j (x, u) ∈ [b ij , b ij ], i ∈ {1, . . . , n x }, j ∈ {1, . . . , n u }.
where the bounds

a ij , a ij , b ij , b ij ∈ R are known. Let us introduce the auxiliary matrices A -, A + ∈ R nx×nx and B -, B + ∈ R nx×nu
, where for all i, j ∈ {1, . . . , n x }

A - ij = a ij if a ij < 0, 0 otherwise, A + ij = a ij if a ij > 0, 0 otherwise,
and for all i ∈ {1, . . . , n x }, for all j ∈ {1, . . . , n u }

B - ij = b ij if b ij < 0, 0 otherwise, B + ij = b ij if b ij > 0, 0 otherwise.
Then, let the functions f -, f + : X × U → R nx and g -, g + :

X × U → R nu be defined for all x ∈ X, u ∈ U , by:

f -(x, u) = f (x, u) -A -x -B -u, f + (x, u) = A + x -B + u -f (x, u), g -(x, u) = g(x, u) + A -x + B -u, g + (x, u) = g(x, u) + A + x + B + u.
Let us remark that g -, g + are known, while f -, f + are unknown but monotone since it can be readily checked that all their partial derivatives are nonnegative. Given a data set D as in [START_REF] Lavaei | Data-driven synthesis of symbolic abstractions with guaranteed confidence[END_REF], we can define two auxiliary data sets:

D -= {(x k , u k , y - k )| y - k = x ′ k -g -(x k , u k ), k ∈ K}, D + = {(x k , u k , y + k )| y + k = g + (x k , u k ) -x ′ k , k ∈ K}.
We can use the approach presented in Section II to compute simulating maps S -, S + : X × U ⇒ R nx of data D -and D + , respectively. Proposition 8: Let S -∈ S D -and S + ∈ S D + , then let S : X × U ⇒ R nx be given for all x ∈ X, u ∈ U by:

S(x, u) = g -(x, u) + S -(x, u) ∩ g + (x, u) -S + (x, u) .
Then, it holds:

∀x ∈ X, u ∈ U, f (x, u) + g(x, u) + W ⊆ S(x, u). Proof: Let x ∈ X, u ∈ U , from S -∈ S D -, we have f -(x, u) + W ⊆ S -(x, u). Since f (x, u) + g(x, u) = f -(x, u) + g -(x, u), we get f (x, u) + g(x, u) + W ⊆ g -(x, u) + S -(x, u).
Similarly, we can show that

f (x, u) + g(x, u) + W ⊆ g + (x, u) -S + (x, u),
which leads to the result of the proposition.

From Proposition 8, a data-driven model of system ( 17) can be defined as follows

x(τ + 1) ∈ S(x(τ ), u(τ )). (30)
Following similar approaches to that described in Sections III-A and III-B, one can define a data-driven transition system T data and a symbolic abstraction T symb such that T symb ⪯ AS T data ⪯ AS T sys .

IV. DATA-DRIVEN VERIFICATION OF ASSUMPTIONS

In the previous sections, we made several assumptions on the unknown map F : bounds on the disturbances and bounds on the partial derivatives of the function f . It is sometimes possible to derive such bounds a priori. For instance, in some cases, the monotonicity of a function can be inferred from principles of physics. However, there are many situations where these bounds also need to be inferred from data. In this section, we briefly present an approach to compute bounds on the disturbances and on the partial derivatives of the function f directly from data, with probabilistic guarantees. Our approach is based on the scenario approach [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF], a data-driven approach to robust convex optimization.

Let us consider a map F : Z ⇒ R m where Z ⊆ R n and such that for all z ∈ Z, F (z) = f (z)+W where f : Z → R m is a differentiable function with unknown lower and upper bounds a ij , a ij on its partial derivatives ∂f i ∂z j and W , W ⊆ R m is a bounded interval of disturbances with unknown lower and upper bounds w, w. Let us remark that without loss of generality, it is always possible to choose w = -w.

For all i ∈ {1, . . . , m}, let a i = (a i1 , . . . , a in ) and a i = (a i1 , . . . , a in ). Some bounds are consistent with our assumptions if and only if they satisfy:

∀z,z ′ ∈ Z, ∀y ∈ F (z), ∀y ′ ∈ F (z ′ ), y i -y ′i ≤ a i • [z -z ′ ] + + a i • [z -z ′ ] -+ 2w i (31)
where

[z] + = max(z, 0) and [z] -= min(z, 0).
Let us assume that we are given a set of random data generated by the map F :

D = {(z k , y k , z ′ k , y ′ k ) | y k ∈ F (z k ), y ′ k ∈ F (z ′ k ), k ∈ K}
where K is a finite set of indices. We assume that the samples in D are independent and identically distributed. We aim at computing bounds such that (31) holds with probabilistic guarantees. This can be done using the scenario approach [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF], which essentially consists in computing bounds such that the inequality in [START_REF] Saoud | Contract-based design of symbolic controllers for safety in distributed multiperiodic sampled-data systems[END_REF] holds at all points in D. However, to obtain high confidence bounds, one needs to consider large data sets, resulting in large linear programs that can be complicated to solve in practice. For that reason, we present the following two-step approach that allows us to deal with very large data sets in practice. Let us partition the set of D in two subsets D′ and D′′ indexed by K′ and K′′ such that K = K′ ∪ K′′ . In the first step, D′ will be used to estimate the bounds a ij , a ij . In the second step, D′′ will be used to estimate the bound w i . Typically, the number of samples in D′′ will be much larger than that in D′ . We first consider the following linear program:

       min a i ,a i ,w i a i • k∈ K′ [z k -z ′ k ] + + a i • k∈ K′ [z k -z ′ k ] -+ 2| K′ |w i s.t. y i k -y ′i k ≤ a i • [z k -z ′ k ] + +a i • [z k -z ′ k ] -+ 2w i , ∀k ∈ K′ a i ⪯ a i . ( 32
)
The constraints of ( 32) are chosen such that [START_REF] Saoud | Contract-based design of symbolic controllers for safety in distributed multiperiodic sampled-data systems[END_REF] holds at all points in D′ , while the cost is chosen so as to minimize the average value of the right-hand side of [START_REF] Saoud | Contract-based design of symbolic controllers for safety in distributed multiperiodic sampled-data systems[END_REF]. Let a ⋆ i , a ⋆ i denote the optimal values of a i , a i in [START_REF] Sinyakov | Formal controller synthesis from specifications given by discrete-time hybrid automata[END_REF]. The optimal value of w i is disregarded and estimated again in the second step.

Indeed, let us define

w ⋆i = 1 2 max y i k -y ′i k -a ⋆ i • [z k -z ′ k ] + -a ⋆ i • [z k -z ′ k ] -k ∈ K′′ . (33
) Let us remark that while (32) can be solved using a moderately large data set D′ , [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF] can easily be estimated on very large data sets. Moreover, we can provide probabilistic guarantees on the estimated bounds: Proposition 9: Let β ∈ (0, 1) be a confidence parameter and let i ∈ {1, . . . , m}, let a ⋆ i , a ⋆ i and w ⋆i be estimated bounds using [START_REF] Sinyakov | Formal controller synthesis from specifications given by discrete-time hybrid automata[END_REF] and [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF], respectively. Then, with probability at least 1 -β, it holds that

P(y i -y ′i ≤ a ⋆ i •[z-z ′ ] + +a ⋆ i •[z-z ′ ] -+2w ⋆i ) ≥ β 1 | K′′ | . ( 34
)
Proof: Consider the following robust linear program

   min w i w i s.t. y i -y ′i ≤ a ⋆ i • [z -z ′ ] + + a ⋆ i • [z -z ′ ] -+ 2w i , ∀z, z ′ ∈ Z, ∀y ∈ F (z), ∀y ′ ∈ F (z ′ ). (35)
Following the scenario approach, we consider the following associated scenario linear program:

   min w i w i s.t. y i k -y ′i k ≤ a ⋆ i • [z k -z ′ k ] + +a ⋆ i • [z k -z ′ k ] -+ 2w i , ∀k ∈ K′′ . ( 36
)
It is easy to see that [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF] provides the unique solution to [START_REF] Zonetti | Decentralized monotonicity-based voltage control of dc microgrids with zip loads[END_REF]. Moreover, from the results in [START_REF] Campi | The exact feasibility of randomized solutions of uncertain convex programs[END_REF], we get that (34) holds for the solution of [START_REF] Zonetti | Decentralized monotonicity-based voltage control of dc microgrids with zip loads[END_REF].

Hence, we have provided a method to estimate bounds on the disturbances and bounds on the partial derivatives. By using a two-step approach, our method can be used with very large data sets and can thus estimate the bounds with very high confidence. Let us remark that the method can be used to check if the system is monotone by checking if a ⋆ i ⪰ 0, for all i ∈ {1, . . . , m}.

V. NUMERICAL EXAMPLES

In this section, we present three numerical examples. In the first, we test the performance of the introduced overapproximation. In the second, we study the case of a system with an unknown monotone part. Finally, we finish by showing the case of a system with bounded derivatives.

For all of the examples, we use interval domains, Z = [α, α], α, α ∈ R n for the maps we are trying to overapproximate. For i ∈ {1, . . . , n}, we define the finite partition (D i q i ) q i ∈Q i as follows, Q i = {0, . . . , K i } and

         D i 0 = [α i , α i 1 ) D i q i = [ q i -1 K i -1 (α i K i -α i 1 ) + α i 1 , q i K i -1 (α i K i -α i 1 ) + α i 1 ), q i = 1, . . . , K i -1 D i K i = [α i K i , α i ]
For the three examples, we choose

α i 1 = α i + c (α i -α i ), α i K i = α i -c (α i -α i )
, and c is a constant specific to each example.

A. Over-approximating an unknown function

In the first example, we present a set of experiments envisioned to test and visualize the algorithms introduced in this paper for the over-approximation of set-valued maps.

To quantitatively measure the performance of the overapproximation, we can check the execution time and the conservatism in the resulting over-approximation. We calculate the conservatism of the over-approximation using the following performance criterion

µ(D, Q) = q∈Q ′ (vol(Z q ) × σ(q)) q∈Q ′ (vol(Z q ) × vol(W ))
where

Q ′ = q ∈ Q -∞ < σ i (q), σ i (q) < ∞, ∀i ∈ {1, . . . , n} .
The denominator of µ represents the volume of the graph of the unknown map for the part of space where we can find an over-approximation, whereas the numerator represents the volume of our over-approximation. µ can take its value in the interval [1, ∞), and the smaller its value is, the less conservative the over-approximation is.

In this example, we consider a monotone set-valued map

F : [-π, π] × [-π, π] ⇒ R given by F (z) = {2 (sin z 1 + sin z 2 + z 1 + z 2 ) + w | w ∈ [-0.1, 0.1]}
(37) We use F to generate the sets of data D used in the subsequent experiments. First, we visualized the over-approximation. We sampled |K| = 10 6 data points. The parameters of the partition are chosen as follows K 1 = 30, K 2 = 30, c = 0.01. Figure 8 shows the undisturbed function in solid and the over-approximation calculated from data. We see how the undisturbed function is included in the over-approximation. 1) The effect of changing the number of data points: To study the effect of changing the number of data points, we chose and fixed a partition, K 1 = 100, K 2 = 100, c = 0.01. Then, for an increasing number of data points, we calculated the over-approximation of the map F and measured the execution time and the performance criterion. For each number of data points, we redo the experiment a hundred times using different randomly generated data sets. The results of this statistical study of changing the number of data points are shown in Figure 9. We can see from the figure the linear relation between the number of data points and the execution time as predicted in Proposition 6. We also see how the conservatism in the calculated over-approximation decreases with the increase in the number of data points. Note that even if the number of points increases toward infinity, the performance criterion µ will not reach one because we are using a fixed partition. For µ to reach one, both the number of data points and the number of partition elements should go to infinity.

2) The effect of changing the size of the partition: In the second experiment, we fixed the number of data points |K| = 10 6 , and changed the size of the partition. For each size considered, the partition is chosen such that K 1 = K 2 , and c = 0.01. We also redo the calculation for each size a hundred times using different randomly sampled data sets. Figure 10 shows the results of the experiment. Time of execution increases logarithmically when |Q| the number of the elements is small. Then, the relation becomes linear for big values of |Q|. This behavior can be justified by the complexity relation introduced in Proposition 6. Similar to the first experiment, the performance criterion decreases with the increase in the number of partition elements.

3) Comparison with the state-of-the-art robust models: We have already shown both theoretically in Proposition 6 and experimentally the computational complexity of calculating the introduced data-driven model. Let us compare with the set membership approaches in [START_REF] Novara | Set membership identification of nonlinear systems[END_REF] and [START_REF] Canale | Nonlinear model predictive control from data: a set membership approach[END_REF], where an optimization problem is used to find the model. The complexity of this approach grows polynomially with the number of data points. In contrast, our approach scales linearly. The maximum reported number of data points used in [START_REF] Novara | Set membership identification of nonlinear systems[END_REF] is in the tens of thousands. The number of points used is less in [START_REF] Canale | Nonlinear model predictive control from data: a set membership approach[END_REF]. As can be seen, the introduced approach in this paper can handle orders of magnitude more data points. The same can be said about [START_REF] Sadraddini | Formal guarantees in data-driven model identification and control synthesis[END_REF], where a Mixed Integer Linear Programming (MILP) problem is used to find the model, and the reported number of data points used in the case study is 400.

B. Cruise control problem

The second example showcases a system that can be seen as a sum of a known function and an unknown monotone function. We find the abstraction representing the system. Then, we use this abstraction to synthesize a safety controller.

Let us consider two vehicles moving in one lane on an infinite straight road. The leader is uncontrollable (vehicle 2), whereas the follower is controllable (vehicle 1). A discretetime model of this setup is given by equations:

d(k + 1) = d(k) + (v f (k) -v l (l))T 0 v f (k + 1) = v f (k) + γ(u(k), v f (k))T 0 v l (k + 1) = v l (k) + w(k)T 0 (38)
Here u is the control input, and d is the signed distance between the vehicles. v l , v f are the velocities of the leader and follower, respectively. The function γ represents the follower Controllable set of the symbolic controller enforcing the assume-guarantee contract (39) for system (38) vehicle acceleration caused by the control input and the friction forces acting upon it. The term w(k) ∈ [W , W ] accounts for uncertainty in the leader velocity and is considered as a disturbance. Model (38) can be seen as the sum of a known part, which can be inferred from the physics of the system, and an unknown or hard-to-model part, namely the function γ(u(k), v f (k))T 0 . We are also able to study the monotonicity of the function γ(u(k), v f (k))T 0 starting from the physics of the system. The acceleration of the car will increase with the increase of the input, and the friction forces will increase with the increase of the velocity. Therefore, by making the change of variable v ′ f = -v f we can apply the algorithm introduced in this work to find an abstraction of the system.

Function γ is given as follows

γ(u, v) = u -M -1 C (f 0 + f 1 v + f 2 v 2 ). The vector of parameters f = (f 0 , f 1 , f 2 ) ∈ R 3
+ describes road friction and vehicle aerodynamics whose numerical values are taken from [START_REF] Nilsson | Correct-by-construction adaptive cruise control: Two approaches[END_REF]: f 0 = 51 N, f 1 = 1.2567 N s/m, f 2 = 0.4342N s 2 /m 2 . For the car mass, we chose M C = 1370 kg. Other numerical values related to the model are X

= [-50, -5] × [0, 30] × [5, 25], U = [-3, 3], T 0 = 0.7 s, W = -W = 2 m/s 2 .
To build the abstraction we fixed a partition K 1

x = 50, K 2 x = 50, K 3 x = 50, K 1 u = 50, c = 0.005, and sampled a set containing |K| = 10 6 data points. Finding the over-approximation of the unknown monotone part was done in 0.633s.

To show the usefulness of the calculated abstraction, we use it to synthesize a safety controller. In this example, the leader's vehicle is uncontrollable. Therefore, we considered a specification given as an assume-guarantee contract:

∀k ∈ N v l (k) ∈ [v l , v l ] =⇒ (39) ∀k ∈ N, v f (k) ∈ [v f , v f ] ∧ d(k) ∈ [d, d].
If the velocity of the leader remains within the bounds [v l , v l ], then the velocity of the follower and the distance between the two vehicles should remain within the bounds [v f , v f ] and [d, d], respectively. The synthesis of symbolic controllers enforcing assume-guarantee contracts such as (39) has been considered in [START_REF] Saoud | Contract-based design of symbolic controllers for safety in distributed multiperiodic sampled-data systems[END_REF]. For this example, we chose

d = α 1 1 , d = α 1 K 1 , v f = α 2 1 , v f = α 2 K 2 , v l = α 3 1 , v f = α 3 K 3 .
The controllable set of the resulting symbolic controller is shown in Figure 11.

C. Lorenz system

In the last example, we study the Lorenz system, a system with bounded derivatives. We first compute the bounds on the disturbance and on the derivatives of the function representing the system from the data. Then, we find a data-driven abstraction that we use to synthesize a controller for the system. The controlled Lorenz system has three states and one input. We built the abstraction using a large number of points and a partition with a large number of elements, and we showed that we could do it efficiently.

The discrete-time controlled Lorenz system is described by the following model The discretization time was chosen to be T 0 = 0.01s. First, we used the algorithm described in Section IV to calculate the bounds on the disturbance and on the derivatives. In the first step, we used 5 * 10 4 data points sampled randomly. In the second step, we sampled 10 6 data points. The results of the calculated bounds on Jacobian matrices and the disturbance are shown below. Subscript D refers to matrices calculated from data, while their unsubscripted counterparts are the actual matrices we are estimating. The results are rounded to the third significant digit. The execution time of the first step is 12.4s, whereas the execution time for the second step is 0.5s. From the results, we can see that all values but the disturbances added to the second and third states are similar to the actual values. The second and third values of the disturbance are more conservative. If we choose β = 10 -6 , then according to Proposition 9, with a probability at least 1 -β, it holds that the probability of We used these values to build the data-driven abstraction. We sampled 10 8 data points. We choose the parameters of the partition as follows K 1

x = 100, K 2 x = 100, K 3 x = 100, K 1 u = 100, c = 0.05. The execution time to reach the abstraction was 929.2 s. To test the validity of the calculated abstraction, we used it to find the maximal safe controlled invariant of the system. Given a safe set X s ⊆ X, a safe controlled invariant I s is a set included in the safe set I s ⊆ X s , and that can be rendered invariant using a suitable controller. Using the finite-state abstraction, we calculated earlier, we can apply an iterative algorithm to find the maximal safe controlled invariant [START_REF] Tabuada | Verification and control of hybrid systems: A symbolic approach[END_REF]. We chose the safe set to be X s = X. The resulting maximal control invariant is represented in Figure 12.

VI. CONCLUSION

In this work, data-driven over-approximation of monotone maps has been introduced. We have proved that the introduced over-approximating maps are minimal. We have also studied the case of over-approximating a monotone map on a fixed partition. We have developed an efficient algorithm to find those maps and then used this data-driven over-approximation to build models for partially unknown systems, where the unknown part is monotone. Using the data-driven models, symbolic abstractions of the system are calculated, which then can be used to synthesize discrete correct-by-design controllers. We have extended our method to study systems with bounded derivatives and introduced a procedure to calculate those bounds from data. In the numerical examples part, we have shown that the time to calculate the over-approximation grows linearly with the number of data points and sub-linearly with the number of elements in the partition. We have also shown the validity of the calculated abstraction by synthesizing safety controllers. In the future, we want to refine and improve the over-approximation using data collected online. We also want to safely explore unknown parts of the state space from which we do not have data points.

Fig. 1 .

 1 Fig. 1. Example of a set-valued monotone map F : Z ⇒ Y .

Definition 2 :

 2 A map F : Z ⇒ R m is consistent with the data D if the following hold:1) There exists a monotone function f : Z → R m such

  and S ′ m are minimal simulating maps of D, then S m = S ′ m a.e.; • If S m is a minimal simulating map of D, S ′ m ∈ S D and S m = S ′ m a.e., then S ′ m is a minimal simulating map of D.

Fig. 3 .

 3 Fig. 3. Partition of Z used to define the function f ⋆ in (8). The area bounded from above by the yellow line represents the set Z + (z ⋆ ).The area bounded from below by the orange line represents the set Z -(z ⋆ ). The area Z 1 in blue contains all the data points smaller than z ⋆ . The area Z 3 in red contains all the data points larger than z ⋆ or incomparable with z ⋆ . The area Z 2 is colored in green.

Fig. 6 .

 6 Fig. 6. The figure shows the predefined partition of the input set Z.The function ϕ maps all the points inside a cell (e.g. square in yellow) to a discrete variable. The map σ is the minimal interval-valued simulating map.

Proposition 6 :

 6 The map σ can be computed with complexity O(|K| × log(|Q|) + |Q|).

Fig. 8 .Fig. 9 .

 89 Fig. 8. Map F (z) with w = 0 everywhere is represented in solid. The upper and lower bounds of the over-approximation are represented in transparency

Fig. 10 .

 10 Fig. 10. Line in blue represents the average time of execution with respect to the number of cells in the partitions, with bars representing the standard deviation. Line in orange represents the performance criterion µ(D, Q) with respect to the number of cells in the partitions.

Fig. 11 .

 11 Fig. 11.Controllable set of the symbolic controller enforcing the assume-guarantee contract (39) for system(38)

  x(k + 1) = x(k) + σ(x(k) -y(k) + w 1 (k))T 0 y(k + 1) = y(k) + (ρ x(k) -x(k)z(k)y(k) + u(k) + w 2 (k))T 0 z(k + 1) = z(k) + (x(k)y(k) -β z(k) + w 3 (k))T 0 (40)For the particular parameter values of σ = 10, β =8 3 , ρ = 28 the system behaves chaotically[START_REF] Lorenz | Deterministic nonperiodic flow[END_REF].We studied the system on the sets X = [-10, 10] × [-10, 10] × [-10, 10], U = [-200, 200], w = [0.5, 0.5, 0.5].

Fig. 12 .

 12 Fig. 12. The maximal controlled invariant of (40) calculated using the data-driven abstraction

S m (z) -w if z ∈ Z 1 , max(min(y ⋆ , Sm(z ⋆ ) -w), S m (z ⋆ ) -w) if z ∈ Z 2 , Sm(z) -w if z ∈ Z 3 .

This project has received funding from the H2020-EU.1.1. research and innovation programme -ERC-2016-COG under grant agreement No 725144.

A. Makdesi and A.