Voronoi Partitionbased Scenario Reduction for Fast Samplingbased Stochastic Reachability Computation of LTI Systems
Abstract
In this paper, we address the stochastic reachavoid problem for linear systems with additive stochastic uncertainty. We seek to compute the maximum probability that the states remain in a safe set over a finite time horizon and reach a target set at the final time. We employ samplingbased methods and provide a lower bound on the number of scenarios required to guarantee that our estimate provides an underapproximation. Due to the probabilistic nature of the samplingbased methods, our underapproximation guarantee is probabilistic, and the proposed lower bound can be used to satisfy a prescribed probabilistic confidence level. To decrease the computational complexity, we propose a Voronoi partitionbased to check the reachavoid constraints at representative partitions (cells), instead of the original scenarios. The state constraints arising from the safe and target sets are tightened appropriately so that the solution provides an underapproximation for the original samplingbased method. We propose a systematic approach for selecting these representative cells and provide the flexibility to tradeoff the number of cells needed for accuracy with the computational cost.
1 Introduction
Reachavoid analysis is an established verification tool for discretetime stochastic dynamical systems, which provides probabilistic guarantees on the safety and performance [1, 2, 3, 4, 5]. This paper focuses on the finite time horizon terminal hitting time stochastic reachavoid problem [1] (referred to here as the terminal time problem), that is, computation of the maximum probability of hitting a target set at the terminal time, while avoiding an unsafe set during all the preceding time steps using a controller that satisfies the specified control bounds.
The solution to the terminal time problem relies on dynamic programming [1, 6, 7], hence a variety of approximation methods have been suggested in literature. Researchers have looked for scalable approaches to solve this problem using approximate dynamic programming [8, 9], Gaussian mixtures [8], particle filters [9, 4], convex chanceconstrained optimization [4], Fourier transformbased verification [10, 11], Lagrangian approaches [5], and semidefinite programming [12]. Currently, the largest system verified is a dimensional chain of double integrators [10, 11] using Fourier transformbased techniques. Existing methods impose a high computational complexity which makes them unrealistic for realtime applications.
In this paper, we reconsider the samplingbased approach, proposed in [4]. Similar samplingbased approach has been used successfully in robotics [13] and in stochastic optimal control [14, 15, 16, 17]. In the samplingbased stochastic reachavoid problem, we sample the stochastic disturbance to produce a finite set of scenarios, and then formulate a mixedinteger linear program (MILP) to maximize the number of scenarios that satisfy the reachavoid constraints [13, 4]. As expected, the approximated probability will converge to the true terminal time probability as the number of scenarios increases. However, the computational complexity of MILP increases exponentially with the number of binary decision variables (the number of scenarios) [18, Rem. 1] making the MILP formulation practically intractable.
The main contributions of this paper are twofold. We first provide a lower bound on the number of scenarios needed to probabilistically guarantee a userspecified upper bound on the approximation error with a userspecified confidence level using concentration techniques. Using Hoeffding’s inequality, we demonstrate that the number of scenarios that need to be considered is inversely proportional to the square of the desired upper bound on the estimate error. Next, we propose a Voronoibased undersampling technique that underapproximates the MILPbased solution in a computationally efficient manner. This approach allows us to partially mitigate the exponential computational complexity, and provides flexibility to select the number of partitions based on the allowable online computational complexity. We demonstrate the application of the proposed method in a problem of spacecraft rendezvous and docking.
The organization of the paper is as follows: Problem formulation and preliminary definitions are stated in Section 2. Lower bound on the required number of scenarios for the prescribed confidence level is given in Section 3. Section 4 presents the proposed partitionbased method and the approximate solution reconstruction. The performance of the proposed method is investigated on a spacecraft rendezvous maneuvering and docking in Section 5.
2 Problem formulation
We presume and are sets of real and natural numbers, with a length vector of real numbers, and the set of natural numbers between and . For , denotes the transpose of . Vector with all elements 1 is denoted .
2.1 System description
Consider a discretetime stochastic LTI system,
(1) 
with state , input , disturbance at time instant , and matrices assumed to be of appropriate dimensions. We assume that is an independent and identically distributed (i.i.d) random variable with a PDF . Note that we require only to be a probability density function from which we can draw samples, and do not require it to be Gaussian. The system (1) over a time horizon with length can be alternatively written in a “stacked” form,
(2) 
with , , and the concatenated state, input, and disturbance vectors over a length horizon [15]. The matrices , , and may be obtained from the system matrices in (1) (see [15]). Due to the stochastic nature of , the state and the concatenated state vector are random. We define as the probability measure associated with the random vector , which is induced from the probability measure of the concatenated disturbance vector and (2). By the i.i.d. assumption on , is characterized by .
2.2 Stochastic reachavoid problem
We are interested in the terminal time problem [1]. As in [1], we seek openloop control laws, to assure tractability (at the cost of conservativeness [10, 4, 11]). We define the terminal time probability, , for a given initial state and an openloop control , as the probability that the state trajectory remains inside the safety set and reaches the target set at time ,
(3) 
with . The stochastic reachavoid problem is formulated as:
Problem 1.
Remark 1.
In Problem 1, is trivially zero when , irrespective of the choice of the controller.
In [4, 13], a mixedinteger linear program (MILP) was formulated as an approximation of Problem 1 when the safe and the target sets are polytopic. Note that restriction of the safe and the target sets to polytopes is not severe since convex and compact sets admit tight polytopic underapproximations [19, Ex. 2.25]. We will denote the safe set , the target set , and the reachavoid constraint set as
(4a)  
(4b)  
(4c) 
with , , and is constructed using and . Using the “bigM” approach [18, 13, 4], and a samplingbased empirical mean for that replaces by a finite set of random samples,
(5) 
we obtain a MILP approximation to Problem 1.
Problem 2.
A MILP approximation to Problem 1 is given by {maxi*}—s— U∈UN1K∑_i=1^Kz^(i)\addConstraintX^(i)= G_xx_0 + G_u U + G_wW^(i), i∈N_[1,K] \addConstraintFX^(i)≤h + M(1z^(i))1, i∈N_[1,K] \addConstraintz^(i)∈{0,1}, i∈N_[1,K] with the optimal value denoted by , optimal control input , some large positive number, and that are concatenated disturbance realizations sampled from , based on the probability law .
As observed in [4, 18, 13], takes the value if and only if for all . For any sampled trajectory that violates the reachavoid constraint , we have which is encoded by . This concept is illustrated in Figure 1; red crosses indicate the sampled trajectories which fail to remain in and, therefore, their corresponding binary variables are zero. From [4], we have,
(6) 
However, Problem 2 becomes computationally intractable for large values of since the worstcase time complexity of MILP problems exponentially increases in the number of binary decision variables [18, Rem. 1].
2.3 Problem statements
Based on Problem 2, we define the random vector , the concatenation of an i.i.d process consisting of Bernoulli random variables . By definition, the probability measure associated with is .
We will address the following questions:
Question 1.
Given a violation parameter and a risk of failure , characterize the sufficient number of scenarios to guarantee or equivalently , for all .
Question 2.
Given scenarios as characterized in Question 1 to meet desired specifications, construct an underapproximate MILP with scenarios (hence binary decision variables) that results in the terminal time probability estimate with .
We will address Question 1 using Hoeffding’s inequality. By solving Question 1, we seek sufficient number of scenarios that leads to a desired upper bound on the likelihood that the approximate solution exceeds the true solution by some threshold (). Note that a smaller implies less conservatism as well as higher accuracy in estimation, but requires more scenarios for a fixed , as claimed in next section. Then we address Question 2 using Voronoi partitions to reduce the number of scenarios while preserving the original specifications and .
2.4 Voronoi partition and data clustering
Here we introduce some preliminaries on Voronoi partition that we will use in the rest of this paper. Given a set of seeds (centres) , , a Voronoi partition partitions the space to cells such that any point in , , is closer to than the other seeds. Given a set of points in , we use to show the partition of through a Voronoi partition with seeds . We define each cell (may also be referred to as partition in this paper) of as
(7) 
where is the distance of from in any valid metric (Euclidean norm is used in this paper). We denote the number of elements in by .
A given set with points in can be clustered in clusters by finding a set of seeds that minimizes the withincluster sum of squares, as proposed by means method/Lloyd’s algorithm [20]. The withincluster sum of squares, denoted by , simply represents the total sum of squared deviations of points in cells from their seeds, i.e., given a set of points , a set of seeds , and a partition set ,
(8) 
where is an indicator function which is one if belongs to cell , and zero otherwise. Let
(9) 
be the set of optimal seeds. Given , cells of represent the optimal clusters of . Although solving (9) is an NPhard problem [21], efficient algorithms exist to compute a local minima. Starting from an initial guess of , a successive algorithm with time complexity for iterations can be used to update each seed by replacing it with the centroid of its cluster elements [22]. This process continues until convergence or reaches its maximum value. Note that the set of optimal seeds is not necessarily a subset of . In addition, is nonincreasing in , and the number of required clusters can be decided by making a tradeoff between (representing accuracy) and (representing computational complexity).
Lemma 1.
Seed selection and Voronoi configuration are preserved under translation.

Let be the set of optimal seeds of , calculated as in (9). For any fixed vector , represents the optimal seeds of .

For any , let . Then .
3 Scenarios required to meet given failure tolerance
Given i.i.d. for and with probability measure , we define the concatenation of these random variables . The probability measure associated with is . We have the following inequalities from Hoeffding [23, Thm. 1].
Lemma 2.
(Hoeffding’s inequality) Define , and . For any ,
(10) 
Theorem 1.
Given a violation parameter , risk of failure , initial state , and the optimal solution to Problem 2, we have the risk of failure , if
(11) 
Proof.
Let the optimal solution to Problem 1 be (which may not be equal to ). For ,
(12a)  
(12b)  
(12c) 
(13) 
Adding and subtracting , we have
(14) 
where (14) follows from (12c).
Thus, is a subset of
by (12b) and (14).
Using the fact that for any two sets , Hoeffding’s inequality (10), , and (13), we have
To obtain the desired probabilistic guarantee , we require . Solving for , we obtain . ∎
Theorem 1 addresses Question 1. Specifically, choosing at least scenarios, Theorem 1 guarantees that the probability of the event that the MILPbased estimated terminal time probability (the optimal solution to Problem 2) exceeds the true terminal time probability (the optimal solution to Problem 1) by more than is less than (a small value). Here, both and are provided by the user. Note that although and are functions of the time horizon , is independent of the choice of . A similar bound is used for the application of aircraft conflict detection in [24].
4 Partitionbased sample reduction
As implied from the concentration probability bounds given in Theorem 1, Problem 2 typically needs a large number of samples to provide a precise approximation for Problem 1 with a small deviation and small risk of failure . Therefore, solving Problem 2 can be computationally expensive or even intractable for realtime applications. In this section, we address Question 2 by proposing a partitionbased method which provides an underapproximation to Problem 2 with flexible computational complexity, as opposed to the samplingbased approach, presented in Problem 2. To this end, we propose the following MILP problem with binary variables, where can be significantly smaller than and is selected by the user.
Problem 3.
The partitionbased terminal time problem is {maxi*}—s— U∈UN1K∑_j=1^^Kα^(j)^z^(j)\addConstraint^X^(j)= G_xx_0 + G_u U + ψ^(j), j∈N_[1,^K] \addConstraintF^X^(j)≤h ε^(j)+ M(1^z^(j))1, j∈N_[1,^K] \addConstraint^z^(j)∈{0,1}, j∈N_[1,^K] with the optimal value denoted by . Here, is some large positive number, and for are selected representatives (seeds) of uncertainty, computed in a prediction mapping with
is the importance rate of the seed with and . For , is an appropriately designed buffer that guarantees the solution of Problem 3 is a lower bound on the solution of Problem 2.
In Problem 3, the state uncertainty is characterized by seeds where the seed represents scenarios of . Then the reachavoid constraints are only checked at the selected seeds instead of being checked at every scenario, which reduces the number of binary variables and constraints.
4.1 Seed Selection and Buffer Computation
Given a sample set , , and , we define as the set of sampled state trajectories. We desire that the elements remain in reachavoid set . The set can be partitioned into cells, where each cell consists of some of the random state trajectories and is represented by a seed. Figure 2 shows a D partition with 11 seeds.
Lemma 3.
Let be the set of optimal seeds of with that minimizes . Then with
represents the set of optimal seeds for sampled state trajectory set .
Proof.
According to Lemma 3, although the state trajectory is an optimization variable, it can be clustered through the prediction mapping offline, independent of the choice of and .
Lemma 4.
Let a set of points , a set of selected seeds , as defined in Lemma 3, and their Voronoi partition with cells be given. Then, for , every point remains in the original constraint set if , the seed, remains in the buffered constraint set with
(15) 
and
(16) 
Proof.
From (16) we conclude that for all , ,
(17) 
By adding to the right and left sides of (17), we have
(18) 
Consequently, for all , the following holds for any initial state and input trajectory,
(19) 
Denoting and with and , when , it is concluded from (19) that , . Since (19) is valid for all and , the proof is completed. ∎
The buffering concept is illustrated in Figure 3.
Remark 2.
Computing for and is a sorting problem in D and can be executed by worst time complexity of .
Theorem 2.
Proof.
According to the definition of the buffers, if seed , , remains in the buffered constraint set, all points of cell remain in the original constraint set. Otherwise, at most samples belonging to cell may violate the original constraints. Since in Problem 3 the worst case is considered by weighting the seed with , provides a lower bound on . ∎
Remark 3.
If initial state is uncertain, the proposed method can be applied by defining .
Obviously, having more cells results in a higher accuracy and when number of cells tends to the number of samples, tends to . However, this improved accuracy comes at a higher computational cost. We thus have to select by trading off accuracy and computational cost.
4.2 Tightening the Voronoibased terminal time probability estimate
Given obtained from Problem 3, a tighter underapproximation on can be recalculated by simply checking the percentage of the original sampled trajectories that remain in reachavoid set , after applying to the stochastic system. In contrast to solving a large MILP (as done in Problem 2) whose computational complexity grows exponentially with , this improved estimate (see (20)) is obtained by a policy evaluation that has a computational complexity of . The following theorem presents the probability underapproximation proposed in this paper.
Theorem 3.
Proof.
i) Since is the optimal terminal time probability with samples and is the evaluation of an openloop controller over these samples, we conclude that . Equality holds if .
ii) Let , the subset of which were deemed safe by Problem 3. By definition of , for every . In other words, since is the set of original scenarios that fall in the cell, whenever the solution of Problem 3 deems the representative seed safe, all the scenarios within it are safe. Thus, we have at least as big as since there might be other cells that were not deemed safe by Problem 3 but contains scenarios that might be safe . Hence, . ∎
4.3 Implementation
Algorithm 1 describes the proposed Voronoibased method to solve the openloop terminal time problem. Given a sample set with random samples directly drawn from , one can construct and find its optimal seeds. The means method can be used to find the seeds of a Voronoi partition as explained in Section 2.4. In addition, in order to determine the number of required seeds, can be used as a measure of variability of points in a cluster. A smaller implies more compact clusters which reduces the size of defined buffers and the average number of samples in each cell. Note that by increasing , clusters become smaller and the precision of Problem 3 grows. However, eventually, the improvement precision is insignificant compared to the imposed computational complexity. Therefore, we compute WSS as a function of , to explore this tradeoff. We propose that the “knee” of the curve provides an efficient compromise between precision and computational complexity. It is shown experimentally in the next section that the knee of vs. curve can be a good representative of the knee on the vs. . As a result, can be computed and selected in advance.
After selecting based on the required running time for the realtime process or based on the vs. , one can compute the Voronoibased partition and the number of elements in each cell, and then compute the buffers from Lemma 4. All these steps are executed offline (independent of ), while solving Problem 3 and probability reconstruction using (20) is done online (dependent on ).
5 Illustrative Example: Spacecraft Rendezvous
We consider the spacecraft rendezvous example discussed in [4]. In this example, two spacecraft are in the same elliptical orbit. One spacecraft, referred to as the deputy, must approach and dock with another spacecraft, referred to as the chief, while remaining in a lineofsight cone, in which accurate sensing of the other vehicle is possible. The relative dynamics are described by the ClohessyWiltshireHill (CWH) equations as given in [25],
(21) 
The position of the deputy is denoted by when the chief, with the mass kg, is located at the origin. For the gravitational constant and the orbital radius of the spacecraft , represents the orbital frequency. In this example, the spacecraft is in a circular orbit at an altitude of km above the earth.
We define as the system state and as the system input, then discretize the dynamics (21) with a sampling time of s to obtain the discretetime LTI system,
(22) 
The additive stochastic noise, modeled by the Gaussian i.i.d. disturbance , with , and , accounts for disturbances and model uncertainty.
We define the target set and the safe set as in [4],
(23)  
(24) 
with a horizon of . We consider the initial position km, the initial velocity km/s and . The terminal time probability for this problem using existing approaches [10, 4] is known to be , which we assume to be the best openloop controllerbased reachavoid probability estimate.
We set as the number of original samples to estimate the terminal time probability, with guarantees afforded by Theorem 1, and run 100 random experiments in which in each experiment is generated randomly. Simulations are carried out using CVX [26] on a GHz processor Intel Core i with GB RAM. Figure 4a shows the curve (mean value and standard deviation of the results of the 100 experiments) with up to cells. Figure 4b shows the terminal time probability approximation provided by Algorithm 1. As proposed, the “knee” of Figure 4a coincides the “knee” of Figure 4b; improvements in the accuracy of Algorithm 1 are insignificant beyond . In practice, can be selected from one single experiment in which is calculated for a random disturbance set for up to 100 (maximum allowable) cells. The computation of curve shown in Figure 4a, for one experiment using means method, took only about s, hence is reasonable for offline computation to select . The reported time includes the computation time for solving means with to cells. This time, which is associated with offline step, can be further reduced by changing the step size of variation (horizontal axis) or calculating for arbitrary s (e.g. finer steps at the beginning and coarser steps at the end). The run time for the online component of Algorithm 1 is shown in Figure 4c. Since Problem 3 is a mixedinteger linear program, the time complexity exponentially increases exponentially with the number of cells [18, Rem. 1].
We see in Figure 4b, that partitions with to cells provide a reasonable estimate of the terminal time probability, without significant loss of precision. The computed terminal time probability and the mean value of the online run time for , and are reported in Table 1. As desired, Algorithm 1 provides a flexible tradeoff between the accuracy and computation time by selecting a suitable partition, and can be significantly faster than the existing Fourier transform approach [10] and particle filter [4, 13] (which fails to deal with large s due to the exponential complexity of MILP problem).
Figure 5 shows the position trajectory, associated with and , obtained by the Fourier method [10] (blue dots) and the proposed Voronoi partitionbased method (green stars) with cells. Green regions show the uncertainty regions of Voronoi method at different time instants obtained by 2000 original scenarios.
6 Conclusion
In this paper we presented a novel partitionbased method for underapproximating the terminal time probability through sample reduction. By using Hoeffding’s inequality, we provided a bound on the required number of scenarios to achieve a desired probabilistic bound on the approximation error. Furthermore, we proposed a method which clusters the taken scenarios in few cells, each cell represented by a seed, where the number of cells is selected by the user in a systematic manner using the trend of a given curve or based on the desired running time. The proposed method scales easily with dimension since the clustering computational complexity increases linearly with the dimension of data. In addition, the simulation results confirm that the proposed method significantly decrease the running time, and therefore, it can be easily applied to realtime systems.
References
 [1] S. Summers and J. Lygeros, “Verification of discrete time stochastic hybrid systems: A stochastic reachavoid decision problem,” Automatica, vol. 46, no. 12, pp. 1951–1961, 2010.
 [2] B. HomChaudhuri, A. P. Vinod, and M. Oishi, “Computation of forward stochastic reach sets: Application to stochastic, dynamic obstacle avoidance,” in American Control Conf., Seattle, WA, 2017.
 [3] N. Malone, K. Lesser, M. Oishi, and L. Tapia, “Stochastic reachability based motion planning for multiple moving obstacle avoidance,” in Proc. Hybrid Syst.: Comput. and Ctrl., 2014, pp. 51–60.
 [4] K. Lesser, M. Oishi, and R. S. Erwin, “Stochastic reachability for control of spacecraft relative motion,” in Proc. IEEE Conf. Dec. & Ctrl. IEEE, 2013, pp. 4705–4712.
 [5] J. Gleason, A. Vinod, and M. Oishi, “Underapproximation of reachavoid sets for discretetime stochastic systems via Lagrangian methods,” in IEEE Conf. Dec. Ctrl., 2017. [Online]. Available: https://arxiv.org/abs/1704.03555.
 [6] A. Abate, M. Prandini, J. Lygeros, and S. Sastry, “Probabilistic reachability and safety for controlled discrete time stochastic hybrid systems,” Automatica, vol. 44, no. 11, pp. 2724–2734, 2008.
 [7] A. Abate, S. Amin, M. Prandini, J. Lygeros, and S. Sastry, “Computational approaches to reachability analysis of stochastic hybrid systems,” in Proc. Hybrid Syst.: Comput. and Ctrl., 2007, pp. 4–17.
 [8] N. Kariotoglou, K. Margellos, and J. Lygeros, “On the computational complexity and generalization properties of multistage and stagewise coupled scenario programs,” Syst. and Ctrl. Lett., vol. 94, pp. 63–69, 2016.
 [9] G. Manganini, M. Pirotta, M. Restelli, L. Piroddi, and M. Prandini, “Policy search for the optimal control of Markov Decision Processes: A novel particlebased iterative scheme,” IEEE Trans. Cybern., pp. 1–13, 2015.
 [10] A. Vinod and M. Oishi, “Scalable Underapproximation for the Stochastic ReachAvoid Problem for HighDimensional LTI Systems Using Fourier Transforms,” IEEE Ctrl. Syst. Letters., vol. 1, no. 2, pp. 316–321, 2017.
 [11] ——, “Scalable underapproximative verification of stochastic LTI systems using convexity and compactness,” in Proc. Hybrid Syst.: Comput. and Ctrl., 2018, pp. 1–10.
 [12] D. Drzajic, N. Kariotoglou, M. Kamgarpour, and J. Lygeros, “A semidefinite programming approach to control synthesis for stochastic reachavoid problems,” in Int’l Workshop on Applied Verification for Continuous and Hybrid Syst., 2016, pp. 134–143.
 [13] L. Blackmore, M. Ono, A. Bektassov, and B. C. Williams, “A probabilistic particlecontrol approximation of chanceconstrained stochastic predictive control,” IEEE Trans. Robot., vol. 26, no. 3, pp. 502–517, 2010.
 [14] H. Sartipizadeh and T. L. Vincent, “A new robust mpc using an approximate convex hull,” Automatica, 2018.
 [15] H. Sartipizadeh and B. Açıkmeşe, “Approximate convex hull based sample truncation for scenario approach to chance constrained trajectory optimization,” in Proc. American Ctrl. Conf., 2018, pp. 4700–4705.
 [16] G. C. Calafiore and L. Fagiano, “Stochastic model predictive control of LPV systems via scenario optimization,” Automatica, vol. 49, no. 6, pp. 1861–1866, 2013.
 [17] G. C. Calafiore and M. C. Campi, “The scenario approach to robust control design,” IEEE Trans. Autom. Ctrl., vol. 51, no. 5, pp. 742–753, May 2006.
 [18] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” Automatica, vol. 35, no. 3, pp. 407–427, 1999.
 [19] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge Univ. Press, 2004.
 [20] J. A. Hartigan, Clustering algorithms. Wiley, 1975.
 [21] D. Aloise, A. Deshpande, P. Hansen, and P. Popat, “NPhardness of euclidean sumofsquares clustering,” Machine Learning, vol. 75, no. 2, pp. 245–248, May 2009. [Online]. Available: https://doi.org/10.1007/s1099400951030
 [22] J. A. Hartigan and M. A. Wong, “Algorithm as 136: A kmeans clustering algorithm,” J. Royal Statistical Society. Series C (Applied Statistics), vol. 28, no. 1, pp. 100–108, 1979.
 [23] W. Hoeffding, “Probability Inequalities for Sums of Bounded Random Variables,” J. Amer. Statistical Asso., vol. 58, no. 301, pp. 13–30, 1963.
 [24] M. Prandini, J. Hu, J. Lygeros, and S. Sastry, “A probabilistic approach to aircraft conflict detection,” IEEE Trans. Intelligent Transportation Syst., vol. 1, no. 4, pp. 199–220, 2000.
 [25] W. E. Weisel, Spaceflight dynamics. New York, McGrawHill Book Co, 1989, vol. 2.
 [26] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.