Voronoi Partition-based Scenario Reduction for Fast Sampling-based Stochastic Reachability Computation of LTI Systems

Voronoi Partition-based Scenario Reduction for Fast Sampling-based Stochastic Reachability Computation of LTI Systems

Hossein Sartipizadeh, Abraham P. Vinod, Behçet Açıkmeşe, and Meeko Oishi This material is based upon work supported by the National Science Foundation, the Air Force Office of Scientific Research, and the Office of Naval Research. Hossein Sartipizadeh and Behçet Açıkmeşe were supported by Air Force Research Laboratory grant FA8650-15-C-2546 and the Office of Naval Research (ONR) Grant No. N00014-15-IP-00052. Vinod and Oishi were supported under NSF Grant Number CMMI-1254990, NSF Grant No. IIS-1528047, and AFRL Grant No. FA9453-17-C-0087. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.
H. Sartipizadeh (corresponding author) is with University of Texas at Austin, TX, US. Email: hsartipi@utexas.edu.
A. Vinod and M. Oishi are with the Electrical & Computer Engineering, University of New Mexico, Albuquerque, NM, US. Email: aby.vinod@gmail.com; oishi@unm.edu.
B. Açıkmeşe is with the Department of Aeronautics & Astronautics in the University of Washington, Seattle, WA. Email: behcet@uw.edu.

In this paper, we address the stochastic reach-avoid 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 sampling-based 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 sampling-based 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 partition-based to check the reach-avoid 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 sampling-based method. We propose a systematic approach for selecting these representative cells and provide the flexibility to trade-off the number of cells needed for accuracy with the computational cost.

1 Introduction

Reach-avoid analysis is an established verification tool for discrete-time 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 reach-avoid 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 chance-constrained optimization [4], Fourier transform-based verification [10, 11], Lagrangian approaches [5], and semi-definite programming [12]. Currently, the largest system verified is a -dimensional chain of double integrators [10, 11] using Fourier transform-based techniques. Existing methods impose a high computational complexity which makes them unrealistic for real-time applications.

In this paper, we reconsider the sampling-based approach, proposed in [4]. Similar sampling-based approach has been used successfully in robotics [13] and in stochastic optimal control [14, 15, 16, 17]. In the sampling-based stochastic reach-avoid problem, we sample the stochastic disturbance to produce a finite set of scenarios, and then formulate a mixed-integer linear program (MILP) to maximize the number of scenarios that satisfy the reach-avoid 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 two-fold. We first provide a lower bound on the number of scenarios needed to probabilistically guarantee a user-specified upper bound on the approximation error with a user-specified 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 Voronoi-based undersampling technique that underapproximates the MILP-based 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 partition-based 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 discrete-time stochastic LTI system,


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,


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 reach-avoid problem

We are interested in the terminal time problem [1]. As in [1], we seek open-loop 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 open-loop control , as the probability that the state trajectory remains inside the safety set and reaches the target set at time ,


with . The stochastic reach-avoid problem is formulated as:

Problem 1.

Open-loop terminal time problem: {maxi} —s— U∈U^Nr_x_0^U(S, T) p^∗( x_0)= Problem 1 is equivalent to (see [1, Sec. 4]), {maxi} —s— U∈U^N1_S( x_0)E_z^x_0,U[ z ],p^∗( x_0)= where is a Bernoulli random variable with a discrete probability measure induced from .

Remark 1.

In Problem 1, is trivially zero when , irrespective of the choice of the controller.

In [4, 13], a mixed-integer 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 reach-avoid constraint set as


with , , and is constructed using and . Using the “big-M” approach [18, 13, 4], and a sampling-based empirical mean for that replaces by a finite set of random samples,


we obtain a MILP approximation to Problem 1.

Problem 2.

A MILP approximation to Problem 1 is given by {maxi*}—s— UUN1K∑_i=1^Kz^(i)\addConstraintX^(i)= G_xx_0 + G_u U + G_wW^(i), i∈N_[1,K] \addConstraintFX^(i)≤h + M(1-z^(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 reach-avoid 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,


However, Problem 2 becomes computationally intractable for large values of since the worst-case time complexity of MILP problems exponentially increases in the number of binary decision variables [18, Rem. 1].

Figure 1: Sampling-based approach illustration for reach-avoid problem in D. Reach-avoid set is distinguished with lines. Black dots and red crosses represent the sampled state trajectories corresponding to sampled disturbance for a given and that succeed and fail to remain in , respectively. Empirical mean of remaining in reach-avoid set is obtained by dividing the number of samples inside to the total number of samples.

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 under-approximate 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


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 within-cluster sum of squares, as proposed by -means method/Lloyd’s algorithm [20]. The within-cluster 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 ,


where is an indicator function which is one if belongs to cell , and zero otherwise. Let


be the set of optimal seeds. Given , cells of represent the optimal clusters of . Although solving (9) is an NP-hard 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 non-increasing in , and the number of required clusters can be decided by making a trade-off between (representing accuracy) and (representing computational complexity).

Lemma 1.

Seed selection and Voronoi configuration are preserved under translation.

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

  2. For any , let . Then .


The proof is straight forward since both and Voronoi partitions (as given in (8) and (7)) are functions of the relative distance of the points in each cell to their seed. ∎

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 ,

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


Let the optimal solution to Problem 1 be (which may not be equal to ). For ,


Using (12a) and (12b),


Adding and subtracting , we have


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 MILP-based 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 Partition-based 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 real-time applications. In this section, we address Question 2 by proposing a partition-based method which provides an underapproximation to Problem 2 with flexible computational complexity, as opposed to the sampling-based 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 partition-based terminal time problem is {maxi*}—s— UUN1K∑_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 reach-avoid 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 reach-avoid 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.

Figure 2: Partitioning the state uncertainty region of Figure 1 to cells using a Voronoi partition. Larger dots indicate the selected Voronoi seeds. Samples inside each cell are closer to their own seed than other 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 .


Results directly from Lemma 1. Since there is no uncertainty in and , in (2) can be interpreted as a translation term. ∎

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




From (16) we conclude that for all , ,


By adding to the right and left sides of (17), we have


Consequently, for all , the following holds for any initial state and input trajectory,


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.

Figure 3: Buffering process. Cell is shown in green. If , the state trajectory corresponding to the seed of , remains in the buffered constraint , the state trajectory of every sample in will satisfy the original constraint .
Remark 2.

Computing for and is a sorting problem in D and can be executed by worst time complexity of .

Theorem 2.

Let be a set of disturbance samples mapped through the prediction mapping and let be a set of selected seeds. Let , , denote the number of elements of and define as in Lemma 4. Problem 3 provides a lower bound for Problem 2.


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 Voronoi-based 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 reach-avoid 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.

Let and be the optimal values of Problem 2 and Problem 3, respectively, with corresponding optimal solutions and . Define as




i) Since is the optimal terminal time probability with samples and is the evaluation of an open-loop 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

  Input: LTI system (1), safe set , target set , initial state .
  Offline (independent of ):
  1. Generate by taking i.i.d. samples from .

  2. Construct with .

  3. Select based on the required time complexity or from the vs. curve.

  4. Compute , the optimal seeds of , by a clustering method.

  5. Determine with cells from (7).

  6. Compute importance rate vector with the number of elements of for .

  7. Compute for using Lemma 4.

  Online (depends on ):
  1. Solve Problem 3 for .

  2. Compute from (20).

  Output: .
Algorithm 1 Proposed Voronoi-based reach-avoid solution

Algorithm 1 describes the proposed Voronoi-based method to solve the open-loop 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 trade-off. 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 real-time process or based on the vs. , one can compute the Voronoi-based 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 line-of-sight cone, in which accurate sensing of the other vehicle is possible. The relative dynamics are described by the Clohessy-Wiltshire-Hill (CWH) equations as given in [25],


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 discrete-time LTI system,


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],


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 open-loop controller-based reach-avoid 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 mixed-integer linear program, the time complexity exponentially increases exponentially with the number of cells [18, Rem. 1].

(a) (b) (c)
Figure 4: Mean and standard deviation of (a) within-cluster sum of squares (used to select , offline) (b) terminal time probability (c) online run time with increasing number of cells , obtained from experiments with original scenarios.

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 trade-off 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 partition-based method (green stars) with cells. Green regions show the uncertainty regions of Voronoi method at different time instants obtained by 2000 original scenarios.

Method Terminal reach-avoid probability Online run time (s)
Algorithm 1
0.83 0.2
0.8492 0.6
0.8604 2.7
Particle filter [4, 13] (Problem 2), - -
Fourier transform [10] 0.862 66
Table 1: Terminal reach-avoid probability estimate and computation time of existing methods and Algorithm 1
Figure 5: Position trajectory for Fourier algorithm given in [10] and the proposed Voronoi partition-based method with 40 cells.

6 Conclusion

In this paper we presented a novel partition-based method for under-approximating 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 real-time systems.


  • [1] S. Summers and J. Lygeros, “Verification of discrete time stochastic hybrid systems: A stochastic reach-avoid 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 reach-avoid sets for discrete-time 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 multi-stage and stage-wise 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 particle-based iterative scheme,” IEEE Trans. Cybern., pp. 1–13, 2015.
  • [10] A. Vinod and M. Oishi, “Scalable Underapproximation for the Stochastic Reach-Avoid Problem for High-Dimensional 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 reach-avoid 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 particle-control approximation of chance-constrained 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, “NP-hardness of euclidean sum-of-squares clustering,” Machine Learning, vol. 75, no. 2, pp. 245–248, May 2009. [Online]. Available: https://doi.org/10.1007/s10994-009-5103-0
  • [22] J. A. Hartigan and M. A. Wong, “Algorithm as 136: A k-means 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, McGraw-Hill 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.
Comments 0