Data-driven Reachable Set Estimation with Tunable Adversarial and Wasserstein Distributional Guarantees*
arXiv:2604.12654v1 [math.OC] 14 Apr 2026
Georgios Pantazis and Michelle S. Chong Abstract— We study finite horizon reachable set estimation for unknown discrete-time dynamical systems using only sampled state trajectories. Rather than treating scenario optimization as a black-box tool, we show how it can be tailored to reachable set estimation, where one must learn a family of sets based on whole trajectories, while preserving probabilistic guarantees on future trajectory inclusion for the entire horizon. To this end, we formulate a relaxed scenario program with slack variables that yields a tunable trade-off between reachable set size and out-of-sample trajectory inclusion over the horizon, thereby reducing sensitivity to outliers. Leveraging the recent results in adversarially robust scenario optimization, we then extend this formulation to account for bounded adversarial perturbations of the observed trajectories and derive a posteriori probabilistic guarantees on future trajectory inclusion. When probability distribution shifts in the Wasserstein distance occur, we obtain an explicit bound on how gracefully the theoretical probabilistic guarantees degrade. For different geometries, i.e., p-norm balls, ellipsoids, and zonotopes, we derive tractable convex reformulations and corroborate our theoretical results in simulation.
I. I NTRODUCTION A. Motivation Reachable sets describe the potential states a dynamical system can attain over a finite horizon under admissible uncertainty, initial conditions and inputs, and form the foundation for verification, safety analysis, and constrained control [1]. While reachable set estimation is mature for well-modeled systems [2], many modern applications rely on partially known dynamics or are affected by uncertainty sources difficult to model without additional information. At the same time, trajectory data are increasingly abundant, collected from simulation or experiments. These trends have motivated data-driven reachability methods that infer reachable sets directly from sampled trajectories while retaining guarantees that extrapolate beyond the observed trajectories [3], [4]. Recent work on reachable state estimation spans several complementary approaches. One approach estimates reachable sets directly from data using functional or geometric set estimators, including methods based on Christoffel functions with finite-sample guarantees [5] and support vector data description for Lipschitz nonlinear systems [4]. Another approach formulates reachable set estimation as a scenario *This work is funded through the CETPartnership’s ProRES project under the Joint Call 2024, co-funded by the European Commission (Grant Agreement No. 101069750) and participating national funding organizations. G. Pantazis and M. S. Chong are with the Dynamics and Control (D&C) section at the Department of Mechanical Engineering, Eindhoven University of Technology, the Netherlands. Email:{G.Pantazis33, M.S.T.Chong}@tue.nl.
program [6], thereby leveraging randomized optimization to learn reachable sets with probabilistic guarantees for unseen trajectories [7], [8]. [9] studies direct data-driven reachableset approximation using Christoffel functions and conformal prediction, yielding finite-sample coverage guarantees. The work in [10] studies verification of candidate neural reachable tubes via scenario optimization and conformal prediction. More recently, the work in [11] studies datadriven reachability through the holdout method, deriving a posteriori guarantees for reachable sets and reachable tubes and comparing them with wait-and-judge bounds [12]. In parallel, probabilistic reachability has been studied through stochastic models, including Gaussian processes [13] and learning-based models of human-in-the-loop systems [14], [15]. More recently, robustness to training–deployment mismatch has been addressed through surrogate modeling and conformal prediction under distribution shift based on KLdivergence [16]. Despite this progress, three practical challenges are still rarely addressed simultaneously in reachable set estimation: 1) A tunable trade-off between set tightness and reliability, 2) robustness to perturbed trajectory data, and 3) guarantees that degrade gracefully under distribution shifts based on the Wasserstein distance. This paper addresses these challenges in a unified framework with the scenario approach as its backbone. B. Main Contributions Specifically, our contributions with respect to the related literature are as follows: 1) Prior works that leverage scenario optimization for reachable set estimation [7], [8] rely on scenario programs with hard constraints for the inclusion of observed trajectories. [10] focuses on verification of candidate neural reachable tubes based on the samplingand-discarding approach [17], while [11] develops probabilistic guarantees for reachable sets based on the holdout method, comparing the results with the bounds in wait-and-judge scenario optimization [12]. This paper instead leverages the scenario approach with relaxation [18] to develop a reachable set estimation framework with auxiliary decision variables that relax selected trajectories. In this way, it bridges [18] with datadriven reachable set estimation and yields a tunable trade-off between reachable set size and the theoretical guarantees for inclusion of future trajectories over the considered horizon.
2) Existing reachable set estimation methods, based on scenario optimization, treat the sampled trajectories as absolute truth and provide guarantees only for the nominal distribution [7], [8], [10], [11]. Recent robust reachability methods based on robust conformal prediction rely on surrogate modeling and calibration [16] to address this. We instead leverage results from adversarially robust scenario optimization [19] to propose a method for reachable set estimation under perturbed sampled trajectories. Specifically, we model each sampled trajectory as subject to bounded adversarial perturbations, known or estimated, and derive a posteriori bounds on the probability of future trajectory inclusion by solving a robustified version of our problem. If the probability distribution of the state trajectories shifts upon deployment, but lies in a 1-Wasserstein ball around the nominal distribution, we obtain an explicit out-ofdistribution bound. 3) Different geometric shapes of candidate reachable sets yield different reformulations. We thus perform a parametric analysis of the size vs. future trajectory inclusion trade-off for p-normed balls, ellipsoids and zonotopes. The obtained reformulations are convex and amenable to efficient numerical computation. The rest of the paper is organized as follows: In Section II we formulate the problem of learning a reachable set with tunable trade-off between size and future trajectory inclusion. In Section III we develop the framework on adversarially robust reachable set estimation with probabilistic guarantees. In Section IV we derive tractable convex formulations for different set geometries such as p-normed balls, ellipsoids, and zonotopes. In Section V, we focus on the distributional robustness of the estimated reachable sets in the Wasserstein distance. Section VI corroborates our results via simulation studies, Section VII concludes the paper, while Section VIII is the Appendix. Notation: We consider the index set T = {0, . . . , T }, where T denotes the horizon length, and N = {1, . . . , N } is the index set of sample trajectories with N being the number of samples. For a collection of states, we use the equivalent ⊤ ⊤ notation (x⊤ = (xk )k∈T , where xk ∈ Rnx is a 0 , . . . , xT ) column We denote Qvector. Q the cartesian product among L sets L−1 Yi as i=0 Yi . PN = i∈N P denotes the cartesian product of N copies of the same probability distribution P. Vol(B) denotes the volume of a set B. A⊤ denotes the transpose of a matrix A, while det(A) denotes its determinant. Sm ≻0 denotes the set of positive definite symmetric matrices of dimension m × m. Let Bnp x be the p-normed ball of dimension nx , and ∥y∥p the p-norm of a vector y. R+ is the set of nonnegative reals, and 1n the vector of ones of dimension n. The indicator function 1[C] takes the value 1 when condition C is true and zero otherwise. II. P ROBLEM F ORMULATION We consider a system with (possibly) unknown dynamics: xk+1 = f (xk ),
(1)
where xk ∈ Xk ⊆ Rnx denotes the state of the system and Xk denotes the set in which the state xk takes values. The function f is a continuous function that maps a state xk to a unique next state xk+1 . The concatenation of these states defines a trajectory x = (xk )k∈T . In (1), the function f is considered to be unknown and incorporates uncertainty in the state dynamics such as model mismatch or disturbances. If the dynamics f are explicitly affected by an uncertain disturbance, a special case of the random operator f is obtained, i.e., f (xk ) = f˜(xk , wk ), where wk ∈ Rnw is an uncertain disturbance parameter. Note that no assumption on the boundedness of the support set of wk for each timestep k ∈ T is considered, thus allowing for possibly unbounded disturbances. Since we have no knowledge of the dynamics, we assume (i) access to N different state trajectories x(i) = (xk )k∈T , i ∈ N obtained from simulations, synthetic data and/or past experiments. For the case f (xk ) = f˜(xk , wk ), these samples represent the trajectories due to different realizations of the initial condition x0 and the disturbance wk , k ∈ T . The state is thus modeled as a random vector in X = Q trajectory (T +1)nx X ⊂ R that satisfies (1) and is drawn from k k∈T a (possibly) unknown probability distribution P. We impose the following standing assumption: Assumption 1: We have access to N independent sample trajectories x(i) from the probability distribution P. □ Note that, in Assumption 1, the probability distribution P is used to describe the entire trajectory vector. As such correlations among the elements of the state trajectories as a result of the dynamics in (1) are allowed. A. Tunable Reachable Set Estimation: Size vs Robustness We consider a parametric family of sets to be used to estimate the reachable region for the states xk at time k ∈ T defined by Rk (θk ) = {zk ∈ Rnx : gk (zk , θk ) ≤ 0}, where θk ∈ Θk ⊆ Rnθ , and gk : Rnx × Rnθ → Rq is a set of constraints that enforce consistency with the specific choice of reachable set geometry. Depending on this choice, gk takes different forms and the decisions θk have a different geometric meaning. Remark 1: For example, considering the family of p-norm balls for reachable sets, we have that gk (zk , θk ) = ∥zk − ck ∥p − rk , where θk = (ck , rk ) ∈ Rnx +1 with ck being the center of the reachable set and rk its radius. □ We propose the following optimization program over a horizon of length T : " # X X PN,ρ : min S Rk (θk ) + ρ ξi θ∈Θ,ξ∈RN ≥0
k∈T (i)
i∈N
s.t. g x , θ ≤ ξi ,
∀i ∈ N ,
where θ = (θk )k∈T is a vector of concatenated decision variables for the candidate reachable set, ξ = (ξi )i∈N the vector of relaxation variables and g is obtained by taking the maximum over the horizon, i.e., g(x, θ) = maxk∈T gk (xk , θk ). S(Rk (θk )) denotes an appropriate size proxy of the set Rk (θk ), e.g., its volume Vol(Rk (θk )). An alternative size
proxy can be the radius for the p-normed ball or the scaling factors of an ellipsoid or zonotope (see Section IV). To ensure that our approach is computationally feasible, we will focus on convex reachable sets and their corresponding reformulations. The optimal solution of program PN,ρ is ∗ denoted by θN,ρ and feeding this solution to the vector-to-set mapping Rk (θk ) produces the estimated reachable sets for ∗ each time step k ∈ T , denoted by R∗k,N = Rk (θk,N,ρ ). We ∗ ∗ denote the collection of these sets by RN = {Rk,N }k∈T . Remark 2: Taking the limit ρ → ∞, problem PN,ρ recovers the reachable set formulation in [7]. Its main advantage is the introduction of relaxation variables ξi for each sampled trajectory, which provide additional flexibility. In particular, varying ρ allows one to trade off the size of the estimated reachable set against consistency with the collected state trajectories. Larger values of ρ place greater emphasis on trajectory consistency and less on the minimization of reachable set size. As shown in the subsequent developments, this trade-off also affects the probabilistic guarantees for inclusion of future state trajectories over the horizon. □ Having obtained our optimal collection of reachable sets R∗k for each k ∈ T by solving PN,ρ , we now wish to evaluate how well they generalize to new yet unseen state trajectories. Thus, we define the probability of trajectory exclusion. ∗ Definition 1: Let R∗k,N = Rk (θk,N,ρ ) be the estimated reachable set at time k ∈ T obtained from PN,ρ . Then, the probability of trajectory exclusion is defined as: n o V(R∗N ) := P x ∈ X : ∃ k ∈ T such that xk ∈ / R∗k,N , where x = (xk )k∈T is a new yet unseen state trajectory. □ The terms trajectory exclusion and trajectory violation are used interchangeably. In words, V(R∗N ) is the probability that a new trajectory is drawn from P according to the dynamics in (1), such that for some timestep k ∈ T , the state xk is not contained in the estimated set R∗k,N . Though general, this violation metric does not account for state trajectories, perturbed due to noise or poisoning from adversaries, which can occur upon deployment.
III. A DVERSARIALLY ROBUST R EACHABLE S ETS We consider that each sampled state trajectory can be perturbed, due to noise or the presence of an adversary, from its nominal value by a radius R according to a well-defined distance metric. The radius R is considered to be either a designer’s choice or an approximation based on adversarial observations, and depicts the trust the designer puts in the sampled trajectory data. As such, the set of adversarial perturbations around a nominal sample x is given by Ax := {x̃ : d(x, x̃) ≤ R}, where d : Rn × Rn → R≥0 is a distance metric such as the p-normed difference d(x, x̃) = ∥x − x̃∥p . We thus propose an extension of our reachable set estimation problem PN,ρ that accounts for trajectory perturbations upon
deployment: " PN,ρ,A :
X
min
θ∈Θ,ξ∈RN ≥0
# X S Rk (θk ) + ρ ξi
k∈T
s.t. g x̃(i) , θ ≤ ξi ,
i∈N
∀x̃(i) ∈ Ax(i) , ∀i ∈ N .
In our setting, we assume that either Ax is known or an approximation Âx can be obtained. More information with regard to methods used to construct such adversarial set approximations can be found in [19]. Hence, program P̂N,ρ,A can be approximated by the program: " # X X P̂N,ρ,A : min S Rk (θk ) + ρ ξi θ∈Θ,ξ∈RN ≥0
k∈T
(i,j)
s.t. g x
, θ ≤ ξi ,
i∈N
∀j ∈ Mi , ∀i ∈ N .
where Mi is the set of vertices of a polytopic approximation of the adversarial set Âx(i) . In case the adversarial set is a known polytope, the approximation is exact. Note that due to the auxiliary variables ξi , an optimizer of P̂N,ρ,A always exists. Depending on the geometric shapes selected for Rk (θk ), k ∈ T , which results in different constraints g(x, θ) (see Remark 1), solving P̂N,ρ,A can be computationally challenging. To this end, we impose the following assumption. Assumption 2: For each k ∈ T , the size proxy S(Rk (θk )) is convex in θk . Moreover, for each fixed x, the constraint function g(x, θ) is convex in θ. □ Assumption 2 ensures that the scenario program P̂N,ρ,A is convex. However, convexity of the cost and constraints alone does not guarantee uniqueness of the optimizer (see [20, Table 1]). Hence, we assume the following: Assumption 3: Problem P̂N,ρ,A admits a unique optimal ∗ solution θN,ρ,A . □ Assumption 3 implies that a unique reachable set R∗k,N,ρ,A is returned for each time step k ∈ T . In case multiple solutions exist, the use of a convex tie-break rule, such as the minimum norm or lexicographic order, can be used to extract a unique solution. Based on this setting, we introduce the concept of adversarial probability of trajectory exclusion: Definition 2: Let R∗k,N,ρ,A be the unique estimated reachable set at time k ∈ T obtained from P̂N,ρ,A and R∗N = {R∗k,N,ρ,A }k∈T . Then, the adversarial probability of trajectory exclusion is defined as: n o VA (R∗N ) := P x : ∃ x̃ ∈ Ax , ∃k ∈ T s.t. x̃k ∈ / R∗k,N,ρ,A . Note that the collection of reachable sets obtained by solving P̂N,ρ,A can be written as R∗k,N,ρ,A = {xk ∈ Rnx : ∗ gk (xk , θk,N,ρ,A ) ≤ 0}. We now introduce the notion of adversarial complexity appropriately adapted from [19] into our reachable set estimation framework: Definition 3: Given the collection of state trajectories {x(i) }i∈N , define the adversarial complexity s∗A as the number of indices i ∈ N such that at least one of the following holds: (i) There exists x̃ ∈ Âx(i) such that for some k ∈ T it ∗ holds that gk (x̃k , θk,N,ρ,A ) > 0;
(ii) There exists x̃ ∈ Âx(i) such that for some k ∈ T it ∗ holds that gk (x̃k , θk,N,ρ,A ) = 0; (iii) There exists x̃ ∈ Ax(i) such that for some k ∈ T it ∗ holds that gk (x̃k , θk,N,ρ,A ) > 0. Definition 3 follows by the observation that {x : ∃k ∈ T such that xk ∈ / R∗k,N,ρ,A } = {x : ∃k ∈ T such that gk (xk , θk ) > 0} and then invoking Definition 7 in [19]. The adversarial complexity s∗A denotes the amount of state trajectories which are important for the construction of the reachable sets across the horizon. In general, a relatively low value of s∗A for a fixed N implies a better ability of the reachable sets to include yet unseen trajectories (see Theorem 1). Furthermore, we impose the following non-accumulation assumption. Assumption 4: For every decision θ ∈ Θ, it holds that: n o P x ∈ X : ∃ x̃ ∈ Âx such that g(x̃, θ) = 0 = 0. □ Assumption 4 prevents degeneracies in the scenario problem and is a standard assumption for uncertainties following a continuous probability distribution, where exact equalities occur with probability zero. We now define the equations for the violation levels to be used for our probabilistic result in Theorem 1, based on [19] and [18]. Definition 4: For each ν ∈ {0, 1, . . . , N −1}, consider the polynomial equation in the variable t ≥ 0: N −1 4N i i−ν N N −ν β X i i−ν β X t − t = 0. t − 2N i=ν ν 6N ν ν i=N +1 (2) Let its two nonnegative roots be t(ν) ≤ t(ν). For ν = N , consider 4N β X i 1− t i−N = 0, (3) 6N N i=N +1
which has a unique nonnegative root t(N ), and set t(N ) := 0. Define, for all ν ∈ N , ϵ(ν) := 1 − t(ν) and ϵ(ν) := max{0, 1 − t(ν)}. Equations (2) and (3) analyzed in [18], [19] establish a connection among the number of samples N , the complexity s∗A (written as ν), and the confidence parameter β, and determine the theoretical violation levels ϵ and ϵ. To solve (2) and (3) efficiently, a bisection algorithm can be used [18, Appendix]. The following result then holds: Theorem 1: Consider system (1), its estimated reachable ∗ set R∗k,N,ρ,A := Rk (θk,N,ρ,A ) with R∗N := {R∗k,N,ρ,A }k∈T obtained via P̂N,ρ,A under Assumptions 1-4. Let s∗A be the adversarial complexity from Definition 3. Then, for any confidence parameter β ∈ (0, 1), it holds that: o n PN {x(i) }i∈N : ϵ(s∗A ) ≤ VA (R∗N ) ≤ ϵ(s∗A ) ≥ 1 − β, where the violation levels ϵ(·) and ϵ(·) are obtained from (2) and (3). □ Proof : See Appendix. ■ Theorem 1 states that for a number of samples N and a parameter β ∈ (0, 1), with confidence at least 1 − β, the adversarial probability of trajectory exclusion VA (R∗N )
lies between ϵ(s∗A ) and ϵ(s∗A ), bounds that exhibit a monotonically increasing relation with respect to the adversarial complexity s∗A . IV. T RACTABLE R EFORMULATIONS In order, for Theorem 1 to be useful, problem P̂N,ρ,A should be tractable. This section derives tractable reformulations of P̂N,ρ,A for different choices of geometric shapes, commonly used in reachable set estimation. We impose one of the following assumptions with regard to the structure of the reachable set: Assumption 5: For each k ∈ T , the reachable sets are pnormed balls, i.e., Rk (ck , rk ) := {x ∈ Rnx : ∥x − ck ∥p ≤ rk } with radius rk ∈ R≥0 , and centre ck ∈ Rnx . □ Assumption 6: For each k ∈ T , the reachable sets are ellipsoids, i.e., Rk (Ck , bk ) := {x ∈ Rnx : ∥Ck x + bk ∥2 ≤ x and bk ∈ Rnx . □ 1}, where Ck ∈ Sn≻0 Assumption 7: For each k ∈n T , the reachable sets are zonotopes, i.e., Rk (ck , ak ) := x ∈ Rnx : ∃ζ ∈ Rm , x = o ck + Gk ζ, |ζ| ≤ ak , where Gk ∈ Rnx ×m is a generator matrix, ζ ∈ Rm a scaling factor, ck ∈ Rnx the centre and ak ∈ Rm + the half-widths. The inequality |ζ| ≤ ak is understood component-wise. □ A. Reformulations based on volume Denote the set of indices K := {(j, k, i) : i ∈ N , k ∈ T , j ∈ Mi }. Then, the following exact reformulations are obtained. Lemma 1: i) Under Assumption 5, consider S(Rk (ck , rk )) = Vol(Rk (ck , rk )). Furthermore, denote c = (ck )k∈T , r = (rk )k∈T , ξ = (ξi )i∈N . Then problem P̂N,ρ,A can be written as: ! X X n n x x min Vol(B )r + ρ ξi (c,r,ξ)∈Rnx T ×RT ×RN ≥0 ≥0
s.t. x(i,j) − ck k
p
p
k
k∈T
i∈N
≤ rk + ξi , ∀(j, k, i) ∈ K.
ii) Under Assumption 6, consider S(Rk (Ck , bk )) = Vol(Rk (Ck , bk )) and denote C = (C0 ; . . . ; CT ), b = (bk )k∈T and ξ = (ξi )i∈N . Then, P̂N,ρ,A can be written as: ! X X Vol(Bnx ) 2 +ρ ξi minnx n Q det(Ck ) S≻0 ×R x , (C,b)∈ i∈N k∈T k∈T
ξ∈RN ≥0
(i,j)
s.t. Ck xk + bk ≤ 1 + ξi , ∀(j, k, i) ∈ K. 2 Proof: See Appendix. The exact reformulations in Lemma 1 can still be challenging to solve in general. To this end, alternative proxies of size are often used to obtain a computationally tractable reachable region. Note that such regions also enjoy the theoretical guarantees of Theorem 1, though with possibly different adversarial complexity. For example, under Assumption 6, consider the monotone proxy S(Rk ) = log(Vol(Rk )) in
problem P̂N,ρ,A for each k ∈ T . Then P̂N,ρ,A can be written as: ! X X min − log det(Ck ) + ρ ξi (C,b)∈ Q Sn≻0x ×Rnx k∈T
ξ∈RN ≥0 (i,j) s.t. Ck xk + bk
i∈N
k∈T
2
≤ 1 + ξi , ∀(j, k, i) ∈ K.
The reformulation is based on the fact that Ck is positive definite, which implies that det(Ck ) > 0. The logarithm of volume of an ellipsoid then takes the form log(Vol(Rk (Ck , bk ))) = log(Vol(Bn2 x )/ det(Ck )) = log(Vol(Bn2 x ))−log(det(Ck )). Since Vol(Bn2 x ) is a constant, the optimizer of the problem at hand is not affected by the translation this term imposes. B. Reformulations based on size proxies The corresponding proxy reformulations of some common shapes, used in reachable set estimation, are established in the following lemma: Lemma 2: i) Under Assumption 5, consider the proxy S(Rk ) = rk and denote c = (ck )k∈T , r = (rk )k∈T , ξ = (ξi )i∈N . Then, problem P̂N,ρ,A can be written as: ! X X min rk + ρ ξi N (c,r,ξ)∈Rnx T ×RT ≥0 ×R≥0
s.t. x(i,j) − ck k
p
k∈T
i∈N
≤ rk + ξi , ∀(j, k, i) ∈ K.
ii) For ellipsoidal sets of the form Rk (ck , sk ) = {x : x ∥Hk (x − ck )∥2 ≤ sk }, where Hk ∈ Sn≻0 is fixed for each k ∈ T , we consider the proxy S(Rk ) = sk and denote c = (ck )k∈T , s = (sk )k∈T and ξ = (ξi )i∈N . Then, P̂N,ρ,A takes the form: ! X X min sk + ρ ξi N {(c,s,ξ)∈Rnx T ×RT ≥0 ×R≥0
s.t. H x(i,j) − c k k k
k∈T 2
i∈N
≤ sk + ξi , ∀(j, k, i) ∈ K.
iii) Under Assumption 7, consider S(Rk ) = 1⊤ m ak with Gk ∈ Rnx ×m fixed, and denote c = (ck )k∈T , a = (ak )k∈T . For each (j, k,mi) ∈ K, let ζi,k,j = ζi,k,j,1 . . . ζ ∈ R and consider ζ := i,k,j,m ζi,k,j (j,k,i)∈K . Then, P̂N,ρ,A takes the form: X X min 1⊤ ak + ρ ξi m mT +N (c,a,ξ,ζ)∈Rnx T ×R≥0 ×Rm|Mi |T N k∈T i∈N (i,j) s.t. ck + Gk ζi,k,j = xk , ∀(j, k, i) ∈ K, − (ak + ξi 1m ) ≤ ζi,k,j ≤ ak + ξi 1m , ∀(j, k, i) ∈ K. Proof : See Appendix. Remark 3: Note that the penalty parameter ρ should be tuned in accordance with the proxy and geometry used for reachable set estimation, as they directly affect the associated scenario problem. Section VI.A supports this observation by showing that, for zonotopes, larger values of ρ are often required to enforce consistency with the sampled trajectories
than for balls or ellipsoids. However, this comes with the benefit of yielding a tighter estimated set. □ V. ROBUSTNESS TO DISTRIBUTIONAL SHIFTS We first define a distance metric between distributions. Due to its optimal transport properties [21], we use the socalled 1-Wasserstein distance defined as follows: Definition 5: Let (∆, G) be a measurable space endowed with a metric d : ∆ × ∆ → R+ , and let P, Q be probability measures on (∆, G). The 1-Wasserstein distance between P and Q is defined as: Z W (P, Q) := inf d(δ, δ ′ ) π(dδ, dδ ′ ), π∈Π(P,Q)
∆×∆
where Π(P, Q) is the set of all joint probability distributions with marginals P and Q, respectively. □ Definition 6: The 1-Wasserstein ball centered at P with radius µ̃ is defined as: Bµ̃ (P) := { Q probability on (∆, G) : W (P, Q) ≤ µ̃ } . □ We now define the following out-of-distribution probability metric over the Wasserstein ambiguity set. Definition 7: Consider a Wasserstein ball Bµ̃ (P). Then, the out-of-distribution probability of trajectory exclusion for the collection of sets R∗N = {R∗k,N,ρ,A }k∈T is defined as: Vµ̃,P (R∗N ) := sup P̃ x ∈ X : ∃k ∈ T s.t. xk ∈ / R∗k,N,ρ,A . P̃∈Bµ̃ (P)
Based on Definition 7, the following result then holds: Theorem 2: Consider system (1) and the estimated collection of reachable sets R∗N obtained by solving P̂N,ρ,A under Assumptions 1-4. Consider the ball of distributions Bµ̃ (P) of radius µ̃ > 0 around P and the radius of sampled trajectory deviations R > 0. Let s∗A be the adversarial complexity from Definition 3 and ϵ(·) the violation level obtained from (2) and (3). Then, for any β ∈ (0, 1), with confidence at least 1 − β, it holds that: µ̃ Vµ̃,P (R∗N ) ≤ ϵ(s∗A ) + . R Proof : We follow similar steps as in the proof of Theorem 1, and then apply Theorem 5 in [19]. ■ Theorem 2 implies that adversarial robustness against perturbations of state trajectories inherits distributional robustness properties. Note that the guarantees are similar to those of Theorem 1 with an additional term that depends on the ratio between the Wasserstein radius and the adversarial radius. VI. N UMERICAL S IMULATIONS A. Adversarially robust reachable set estimation We consider the following stochastic dynamical system: xk+1 = Axk + Bϕ(Cxk ) + wk ,
(4)
where ϕ(z) = a tanh(z) with a = −0.9, and the matrices: 0.95 0.10 0.18 A= ,B = ,C = 1 0 . −0.20 0.85 0.06 We consider a horizon of length T = 25. We draw initial conditions x0 from the box [−0.6, 0.6] × [−0.45, 0.45]
(a) Adversarially robust reachable set based on zonotopes.
(b) Adversarially robust reachable set based on ellipsoids.
(c) Adversarially robust reachable set based on Euclidean balls.
Fig. 1: Adversarially robust data-driven estimation of reachable sets for different geometries, i.e., zonotopes, ellipsoids and Euclidean balls. We consider T = 25 and illustrate three timestamps at k ∈ {5, 10, 20}. We learn each set using N = 1000 different state trajectories, perturbed by the set ∆, thus giving rise to the red cross points which denote adversarial samples. We consider a penalty parameter ρ ∈ {0.5, 1, 2, 5} and illustrate the geometric shapes for these values. Note that as ρ increases, the set increases as well. This is because the relaxation variables of P̂N,ρ,A are penalized more to satisfy the observed state trajectories. As such, ρ acts as a tuning parameter trading set size for inclusion of future state trajectories.
Fig. 2: Adversarial theoretical trajectory violation level (dashed lines) vs its empirical counterpart (solid lines) for different geometric shapes. As the number of samples N increases, the theoretical and empirical probabilities improve significantly. Note that zonotopes require a higher number of sampled trajectories to provide the same guarantees with ellipsoids and balls, though at the benefit of tighter zonotopic regions.
and the disturbance wk from a uniform probability distri(i) bution with support [−0.05, 0.05]2 . By feeding wk and (i) x0 through the dynamics, we obtain the sample trajectory (i) x(i) = {xk }k∈T . Using N such trajectories, we obtain the multi-sample {x(i) }i∈N and we model the presence of possible data perturbations by considering a perturbation (i) (i) to each sampled trajectory, i.e., x̃k = xk + δ, where ∆ = {δ ∈ R2 : ∥δ∥∞ ≤ γ}, with γ = 0.03. Since ∆ is a rectangle, we can use the vertices to define the set M in P̂N,ρ,A . We now define the empirical adversarial probability of trajectory exclusion as follows:
V̂adv =
1
N test X
Ntest m=1
(m,j)
1[max max gk (xk k∈T j∈M
∗ , θk,N,ρ,A ) > 0],
where Ntest is the number of test samples and M is the index set denoting the vertices of the perturbation set ∆. Leveraging the results from Lemma 2 we solve problem P̂N,ρ,A for system (4) for three different geometric shapes, namely, Euclidean balls, ellipsoids and zonotopes. In Figure 1 we set a horizon of T = 25 and illustrate the learned sets at k ∈ {5, 10, 20}, trained from N = 1000 state trajectories, including adversarially perturbed samples from ∆ (red crosses). Note that perturbations are considered for each sample. We chose to depict only 20 of them in Figure 1 for illustrative purposes. Shapes are shown for ρ ∈ {0.5, 1, 2, 5}. Larger ρ penalizes the relaxation variables more, yielding larger sets that better cover the data. As such, ρ trades set size against out-of-sample performance. Figure 2 shows the adversarial theoretical and the adversarial empirical probability of violation for different number of samples N ∈ {500, 1000} and penalty parameters ρ ∈ {0.5, 1, 2, 5, 7}. Note that the empirical probability of trajectory violation is always below the corresponding theoretical bound.
Fig. 3: Out-of-distribution theoretical trajectory violation level (dashed lines) and out-of-distribution empirical trajectory violation level for reachable sets of different geometries.
B. Distributional robustness through adversarial training In this section we generate 5 different experiments using N = 1000 different trajectory samples per experiment. We consider the same system (4) as in Section VI-A with the initial condition x0 now drawn from a Gaussian distribution N (µx0 , diag(σx0 )) where µx0 = (0, 0)⊤ , and 2 σx0 = (0.3, 0.225)⊤ and wk ∼ N (µw , σw I2 ), with µw = ⊤ (0, 0) and σw = 0.0167. We now consider Ntest = 3000 state trajectories drawn from the perturbed distributions N (µ̂x0 , diag(σ̂x0 )) where µ̂x0 = (0.01, −0.01)⊤ , 2 I2 ), and σ̂x0 = (0.315, 0.23625)⊤ and ŵk ∼ N (µ̂w , σ̂w ⊤ with µ̂w = (0.002, −0.002) and σ̂w = 0.0175. With these samples, we calculate the out-of-distribution empirical probability of violation and compare it with the theoretical guarantees in Theorem 2, where µ̃ = 0.0243 is an upper bound that guarantees that the test distributions lie within this Wasserstein ball. The results are depicted in Figure 3. Note that the empirical out-of-distribution probability of trajectory violation is always below the corresponding theoretical trajectory violation bound for all considered geometries. It is important to note that the distributional robustness is a property of adversarial training. Thus, for significantly different distributions, the second term in the bound of Theorem 2 can increase significantly unless R is tuned differently. Future work will focus on establishing a distributionally robust methodology incorporated in the scenario methodology such that the distributional robustness is explicitly guaranteed. Figure 4 shows the relative cumulative reachable set size for each P geometry, where P Sizeball (ρ) := Pk∈T rk (ρ), Sizeell (ρ) := k∈T sk (ρ), and Sizezono (ρ) := k∈T 1⊤ m ak (ρ). Each quantity is normalized by its corresponding value at ρ = ρ0 = 0.5. Note that as ρ increases, all geometries interestingly exhibit a similar relative increase in cumulative size. VII. C ONCLUSION We propose a methodology for reachable set estimation, based on recent results in the scenario approach, that learns
Fig. 4: Relative cumulative reachable set size over time for each geometry. Interestingly, as ρ increases, all shapes exhibit similar relative increase in cumulative size.
convex reachable sets directly from trajectory data, explicitly balancing set size and generalization properties of the set through a penalty hyperparameter. Furthermore, we obtain a posteriori robustness certificates to bounded adversarial perturbations and Wasserstein distribution shifts. A key next step is to make distributional robustness explicit by embedding a Wasserstein ambiguity set into the training procedure, rather than relying on robustness inherited from adversarial training. VIII. A PPENDIX Proof of Theorem 1: The following equalities hold: n o VA (R∗N ) := P x : ∃ x̃ ∈ Ax , ∃ k ∈ T s.t. x̃k ∈ / R∗k,N,ρ,A n ∗ = P x : ∃ x̃ ∈ Ax , such that max gk (x̃k , θk,N,ρ,A ) > 0} k∈T n ∗ = P x : ∃ x̃ ∈ Ax , such that g(x̃, θN,ρ,A ) > 0}. Application of Theorem 1 in [19] concludes the proof. ■ Proof of Lemma 1: i) The proof follows by setting θ = (i,j) (ck )k∈T ; (rk )k∈T , g(x(i,j) , θ) = maxk∈T ∥xk −ck ∥p − rk in P̂N,ρ,A and taking the volume of the p-normed ball Vol(Rk ) = αnx ,p rknx . ii) The proof follows by setting (i,j) g(x(i,j) , θ) = maxk∈T ∥Ck xk + bk ∥2 − 1 in P̂N,ρ,A and αnx ,2 . taking Vol(Rk ) = det(Ck ) Proof of Lemma 2: i) The proof is similar to the corresponding proof of Lemma 1. ii) If Hk is fixed for all k, (i,j) setting g(x, θ) = maxk∈T ∥Hk (xk − ck )∥2 − sk yields the desired reformulation. iii) For each time k, fix a generator matrix Gk ∈ Rnx ×m . We model the reachable set as the zonotope Rk (ck , ak ) = {x ∈ Rnx : ∃ζ ∈ Rm , x = ck + Gk ζ, |ζ| ≤ ak }, with center ck ∈ Rnx and half-widths m ak ∈ Rm + . Introducing the auxiliary variables ζi,k,j ∈ R , then the above set is equivalent to satisfying the constraints (i,j) ck + Gk ζi,k,j = xk and − (ak + ξi 1m ) ≤ ζi,k,j ≤ ak + ξi 1m componentwise. ■ R EFERENCES [1] M. Althoff, G. Frehse, and A. Girard, “Set propagation techniques for reachability analysis,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 4, pp. 369–395, 2021.
[2] M. S. Chong, “On the ellipsoidal bounds of the reachable set for a class of time-delayed nonlinear systems with bounded input,” 2022 European Control Conference (ECC), pp. 882–887, 2022. [3] A. Alanwar, A. Koch, F. Allgöwer, and K. H. Johansson, “Datadriven reachability analysis from noisy data,” IEEE Transactions on Automatic Control, vol. 68, no. 5, pp. 3054–3069, 2023. [4] Z. Wang, B. Chen, R. M. Jungers, and L. Yu, “Data-driven reachability analysis of Lipschitz nonlinear systems via support vector data description,” Proceedings of the 62nd IEEE Conference on Decision and Control (CDC), pp. 7043–7048, 2023. [5] A. Devonport, F. Yang, L. E. Ghaoui, and M. Arcak, “Data-driven reachability and support estimation with Christoffel functions,” IEEE Transactions on Automatic Control, vol. 68, no. 9, pp. 5216–5229, 2023. [6] G. Calafiore and M. Campi, “The scenario approach to robust control design,” IEEE Transactions on Automatic Control, vol. 51, no. 5, pp. 742–753, 2006. [7] A. Devonport and M. Arcak, “Estimating reachable sets with scenario optimization,” Proceedings of the 2nd Conference on Learning for Dynamics and Control, vol. 120, pp. 75–84, 10–11 Jun 2020. [8] E. Dietrich, A. Devonport, and M. Arcak, “Nonconvex scenario optimization for data-driven reachability,” 6th Annual Learning for Dynamics & Control Conference, pp. 514–527, 2024. [9] A. Tebjou, G. Frehse, and F. Chamroukhi, “Data-driven reachability using Christoffel functions and conformal prediction,” Proceedings of the Twelfth Symposium on Conformal and Probabilistic Prediction with Applications, vol. 204, pp. 194–213, 2023. [10] A. Lin and S. Bansal, “Verification of neural reachable tubes via scenario optimization and conformal prediction,” Proceedings of the 6th Annual Learning for Dynamics & Control Conference, vol. 242, pp. 719–731, 2024. [11] E. Dietrich, R. Devonport, S. Tu, and M. Arcak, “Data-driven reachability with scenario optimization and the holdout method,” 2025 IEEE 64th Conference on Decision and Control (CDC), pp. 3925–3931, 2025. [12] M. C. Campi and S. Garatti, “Wait-and-judge scenario optimization,” Mathematical Programming, vol. 167, no. 1, pp. 155–189, 2018. [13] P. Griffioen and M. Arcak, “Data-driven reachability analysis for gaussian process state space models,” 62nd IEEE Conference on Decision and Control (CDC), pp. 4100–4105, 2023. [14] V. Govindarajan, K. Driggs-Campbell, and R. Bajcsy, “Data-driven reachability analysis for human-in-the-loop systems,” 56th IEEE Conference on Decision and Control (CDC), pp. 2617–2622, 2017. [15] J. Choi, S. Byeon, and I. Hwang, “Data-driven forward stochastic reachability analysis for human-in-the-loop systems,” 62nd IEEE Conference on Decision and Control (CDC), pp. 1730–1735, 2023. [16] N. Hashemi, L. Lindemann, and J. V. Deshmukh, “Statistical reachability analysis of stochastic cyber-physical systems under distribution shift,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 43, no. 11, pp. 4250–4261, 2024. [17] M. C. Campi and S. Garatti, “A sampling-and-discarding approach to chance-constrained optimization: Feasibility and optimality,” Journal of Optimization Theory and Applications, vol. 148, no. 2, pp. 257–280, 2011. [18] ——, “A theory of the risk for optimization with relaxation and its application to support vector machines,” Journal of Machine Learning Research, vol. 22, no. 288, pp. 1–38, 2021. [19] M. C. Campi, A. Caré, L. G. Crespo, S. Garatti, and F. A. Ramponi, “Risk analysis and design against adversarial actions,” arXiv preprint arXiv:2505.01130, May 2025. [20] P. Mestres, A. Allibhoy, and J. Cortés, “Regularity properties of optimization-based controllers,” European Journal of Control, vol. 81, p. 101098, 2025. [21] P. Mohajerin Esfahani and D. Kuhn, “Data-driven distributionally robust optimization using the wasserstein metric: Performance guarantees and tractable reformulations,” Mathematical Programming, vol. 171, no. 1–2, pp. 115–166, 2018.