{'fraction_non_alphanumeric': 0.048434366805491284, 'fraction_numerical': 0.013625379196873875, 'mean_word_length': 4.534073125622421, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SOCKS: A Stochastic Optimal Control and Reachability Toolbox Using Kernel Methods

We present SOCKS, a data-driven stochastic optimal control toolbox based in kernel methods. SOCKS is a collection of data-driven algorithms that compute approximate solutions to stochastic optimal control problems with arbitrary cost and constraint functions, including stochastic reachability, which seeks to determine the likelihood that a system will reach a desired target set while respecting a set of pre-defined safety constraints. Our approach relies upon a class of machine learning algorithms based in kernel methods, a nonparametric technique which can be used to represent probability distributions in a high-dimensional space of functions known as a reproducing kernel Hilbert space. As a nonparametric technique, kernel methods are inherently data-driven, meaning that they do not place prior assumptions on the system dynamics or the structure of the uncertainty. This makes the toolbox amenable to a wide variety of systems, including those with nonlinear dynamics, black-box elements, and poorly characterized stochastic disturbances. We present the main features of SOCKS and demonstrate its capabilities on several benchmarks.


INTRODUCTION
As modern dynamical systems increasingly incorporate learning enabled components, human-in-the-loop elements, and realistic stochastic disturbances, they become increasingly resistant to traditional controls techniques, and the need for algorithms and tools which can handle such uncertain elements has also grown. Because of the inherent complexity of these systems, control algorithms based in machine learning are becoming ever more prevalent, and frameworks such as reinforcement learning (RL) and deep neural network controllers have seen widespread popularity in this areain part because they allow for approximately optimal controller synthesis using a data-driven exploration of the state space and do not rely upon model-based assumptions. Data-driven control techniques present an attractive approach to stochastic optimal control due to their ability to handle dynamical systems which are resistant to traditional modeling techniques, as well as systems with learning-enabled components and black-box elements.
We present SOCKS, a toolbox for data-driven optimal control based in kernel methods. The algorithms in SOCKS use a technique known as kernel embeddings of distributions, a nonparametric technique which is rooted in functional analysis and a class of machine learning techniques known collectively as kernel methods [41,43,46]. Kernel distribution embeddings have been applied to modeling of Markov processes [23,44], robust optimization [57], and statistical inference [45]. In addition, these techniques have also been applied to solve stochastic reachability problems [50,54], forward reachability analysis [52], and to solving stochastic optimal control problems [29,51]. Because these techniques are inherently data-driven, SOCKS can accommodate systems with nonlinear dynamics, black-box elements, and arbitrary stochastic disturbances.
Data-driven stochastic optimal control is an active area of research [17,18], and provides a promising avenue for controls problems which suffer from high model complexity or system uncertainty, such as robotic motion planning [25,30] and model predictive control [39,40]. Recently, approaches using Gaussian processes [16,36] and kernel methods [23,29,51] have also been explored. In SOCKS, we implement the algorithms in [51], which uses data consisting of observations of the system evolution to compute an implicit approximation of the dynamics in a reproducing kernel Hilbert space (RKHS). The novelty of the approach in [51] is that it exploits the structure of the RKHS to approximate the stochastic optimal control problem as a linear program that converges in probability to the original problem, and computes an approximately optimal controller without invoking a model-based approach.
The application areas of stochastic optimal control are often strongly motivated by a need for assurances of safety, which presents a need for optimal control techniques which can account for predefined safety constraints. In the reinforcement learning community, this need has led to the development of learning frameworks which enable guided state space exploration strategies [e.g. 20,38], as well as toolsets which implement safety constraint satisfaction as part of the learning loop, such as Safety Gym [37]. SOCKS can also be used to provide assurances of safety using an established framework known as stochastic reachability [5,49], which seeks to determine the likelihood of satisfying a set of pre-specified safety constraints (also called the safety probability). Numerous toolsets for stochastic reachability have been developed, including [11,15,26,27,42,47,56] (see [2][3][4] for a detailed comparison). In SOCKS, we use the algorithms developed in [50,54], which compute an approximation of the stochastic reachability safety probability using kernel embeddings of distributions. A recent addition to SReachTools [56], presented in [53], implements one of the existing stochastic reachability algorithms in SOCKS, but does not consider the stochastic optimal control problem.
Lastly, SOCKS implements an algorithm for forward reachability analysis presented in [52]. This technique is useful for analyzing systems with black-box elements, such as deep neural network controllers. Because it employs a data-driven approach, it is agnostic to the structure of the network, and can be used for neural network verification. Several toolboxes for reachability analysis and verification of deep neural networks have been presented in [19,24,55]. However, many existing toolsets rely upon prior knowledge of the network structure (such as knowledge of the activation functions), which may not be available without prior knowledge of the system. Because our approach is data-driven, we do not exploit the structure of the system or the network.
The rest of the paper is outlined as follows. In Section 2, we describe the class of systems that SOCKS is designed to handle as well as the problems we consider. In Section 3, we give an overview of the kernel-based techniques used by SOCKS. Section 4 describes the main features of the toolbox. In Section 5, we demonstrate the algorithms in SOCKS on several examples, including a a nonholonomic target-tracking scenario, a realistic satellite rendezvous and docking scenario, a double integrator system to demonstrate stochastic reachability, and on a forward reachable set estimation problem for a neural-network controlled system. Concluding remarks are presented in Section 6.

PRELIMINARIES
We use the following notation throughout: Let ( , E) be an arbitrary measurable space where E is the -algebra on . If is topological and E is the -algebra generated by all open subsets of , then E is called the Borel -algebra and is denoted B X . Let (Ω, F , P) be a probability space, where F is the -algebra on Ω and P : F → [0, 1] is a probability measure on (Ω, F ). A measurable function : Ω → is called an -valued random variable. The image of P under , P( −1 ), ∈ E, is called the distribution of . A sequence of -valued random variables = { | = 0, 1, . . .} is called a stochastic process with state space ( , E).
We define a stochastic kernel according to [13].

System Model & Data
Consider a discrete-time stochastic system, where (X, B X ) is a Borel space, (U, B U ) is a compact Borel space, and are independent and identically distributed (i.i.d.) random variables defined on the measurable space (W, B W ). The system evolves over a time horizon = 0, 1, . . . , , ∈ N, from an initial condition 0 ∈ X, which may be drawn from an initial distribution P 0 on X, with inputs chosen from a Markov policy .

Definition 2 (Markov Policy). A Markov policy is a sequence
We denote the set of all Markov policies as Π, and for simplicity, we assume the policy is stationary, meaning 0 = 1 = · · · = −1 . We can represent the system in (1) as a Markov control process [7]. Definition 3 (Markov Control Process). A Markov control process is a 3-tuple H = (X, U, ), consisting of a Borel space (X, B X ), a compact Borel space (U, B U ), and a stochastic kernel We consider the case where the stochastic kernel is unknown, meaning we have no prior knowledge of the statistical features of (· | , ) or the dynamics in (1), but assume that a sample S collected i.i.d. from is available. We make this scenario more explicit via the following assumptions. Assumption 1. The stochastic kernel is unknown.
where and are randomly sampled from a probability distribution on X × U and ∼ (· | , ).

Problem Definitions
2.2.1 Stochastic Optimal Control. Consider the following stochastic optimal control problem, which seeks to minimize an arbitrary, bounded cost function subject to a set of constraints.
Problem 1 (Stochastic Optimal Control). Let H be a Markov control process as in Definition 3, and define the functions 0 : X × U → R, called the objective or cost function and : X × U → R, = 1, . . . , , called the constraints. We seek a policy ∈ Π that minimizes the following optimization problem: We impose the following mild simplifying assumption which allows us to separate the cost with respect to and . Assumption 3. The cost and constraint functions : X ×U → R, = 0, 1, . . . , , can be decomposed as: Several commonly-known cost functions obey this assumption, such as the quadratic LQR cost function.
The primary difficulty in solving Problem 1 is due to Assumption 1, and also because we seek a distribution which minimizes the objective. Because is unknown, the integral with respect to in (3) is intractable. Thus, we form an approximation of Problem 1 by computing an empirical approximation of the integral operator with respect to using the sample S. Following [51], we can view this as a learning problem by embedding the integral operator as an element in a high-dimensional space of functions known as a reproducing kernel Hilbert space. Details regarding the kernelbased stochastic optimal control method are provided in [51] and in Appendix A.

Backward Stochastic Reachability.
We also consider a special case of the stochastic optimal control problem in (3), known as the terminal-hitting time stochastic reachability problem. As defined in [49], the goal is to compute a policy ∈ Π that maximizes the likelihood that a system H will remain within a pre-defined safe set K ⊆ B X for all time < , and reach some target set T ⊆ B X at time = . We define the safety probability as: The solution to the stochastic reachability problem is typically formulated as a dynamic program using indicator functions. Define the value functions * : X → [0, 1], = 0, 1, . . . , by the backward recursion, where ∈ X. Then * 0 ( 0 ) = sup ∈Π 0 (K, T ). Problem 2 (Terminal-Hitting Time Problem). We seek to compute an approximation of the policy ∈ Π that maximizes the safety probabilities in (5), and converges in probability to the true solution, * = arg sup ∈Π 0 (K, T ).
Similar to Problem 1, the backward recursion in (6) is intractable due to Assumption 1. We can use the same technique as Problem 1 in order to approximate the value functions in (6), and thereby obtain an approximation of the safety probabilities in (5). Remark 1. We note that our toolbox can be used to solve other stochastic reachability problems, including the first-hitting time problem as defined in [49] and the max and multiplicative problems defined in [5]. We focus on the terminal-hitting time problem in the current work for simplicity.

Forward Stochastic
Reachability. The forward reachable set F is defined as the set of all states that the system in (1) can reach after time steps from an initial condition 0 ∈ X. As shown in [52], we can view the problem of estimating the forward reachable set F as a support estimation problem, where the support is the smallest closed set F ⊂ X such that P ( ∈ F ) = 1, where is a random variable representing the state of the system at time and P is the distribution of .
Problem 3 (Forward Reachability). We seek to determine the support of P , the state distribution over (X, B X ) after time steps.
We formulate Problem 3 as learning a classifier , where The difficulty in computing the reachable set classifier is due to the fact that the dynamics in (1) and the stochastic kernel is unknown (by Assumption 1). Thus, we approximate the classifier function as an element in an RKHS, and use a sample collected from P in order to estimate .

EMBEDDING DISTRIBUTIONS IN A HILBERT SPACE OF FUNCTIONS
In this section, we provide an overview of the machine learning techniques used by SOCKS to solve Problems 1 and 2. Details are provided in [50,51,53]. The method used to solve Problem 3 is based in a similar framework, with additional details in [52].
Let : X × X → R be a positive definite kernel function.
Let H denote a Hilbert space of functions from X to R, equipped with the inner product ⟨·, ·⟩ H . Definition 5 (RKHS, [6]). A Hilbert space H of functions X → R is called a reproducing kernel Hilbert space if there exists a positive definite kernel called the reproducing kernel, such that the following properties hold: (1) ( , ·) ∈ H for all ∈ X, and (2) ( ) = ⟨ , ( , ·)⟩ H for all ∈ H and ∈ X.
Alternatively, by the Moore-Aronszajn theorem [6], for any positive definite kernel function , there exists a unique RKHS with reproducing kernel . For instance, a commonly-used kernel function is the Gaussian RBF kernel ( , ′ ) = exp(−∥ − ′ ∥ 2 2 /(2 2 )), > 0. We define the reproducing kernel : X × X → R on X with the associated RKHS H and the kernel : The second property in Definition 5 is called the reproducing property and is key to our approach. In short, it allows us to evaluate any function in H as a Hilbert space inner product. We use this property to evaluate the integral terms in Problems 1 and 2 by embedding the integral operator with respect to the stochastic kernel as an element in an RKHS.
We now define the following: In order to embed a distribution in H , we make the following mild assumption: Assumption 4. The kernel is bounded and measurable with respect to X.
For every ( , ) ∈ X × U, (· | , ) is a probability measure on X. We denote by P the set of probability measures on X conditioned on X × U, of which the probability measures (· | , ) generated by are a part. If the following necessary and sufficient condition is satisfied, which holds due to Assumption 4, then there exists an element ( , ) ∈ H called the kernel distribution embedding, which is a mapping, : Then by the reproducing property, we can evaluate the integral of any function ∈ H as an RKHS inner product, However, in practice, we do not have access to the true embedding ( , ) due to Assumption 1. Thus, we seek an empirical estimateˆ( , ) of ( , ) computed using the sample S.

Empirical Distribution Embeddings
Following [22], we can compute an estimateˆ( , ) using S as the solution to the following regularized least-squares problem: where Q is a vector-valued RKHS [31], and > 0 is the regularization parameter. According to [31], by the representer theorem, the solution to (14) has the following form: where ∈ R is a vector of real-valued coefficients. By substituting (15) into (14) and taking the derivative with respect to , we obtain the following closed-form solution: where Φ and Ψ are called feature vectors with elements Φ = ( , ·) and Ψ = ( , ·) ( , ·). For simplicity, let (ΨΨ ⊤ + ) −1 and define such thatˆ( , ) = Φ ⊤ ( , ). Then by the reproducing property, for any function ∈ H , we can approximate the integral of with respect to (· | , ) as an RKHS inner product: where is a vector with elements = ( ).
In addition, the empirical estimateˆ( , ) converges in probability to the true embedding ( , ) as the sample size increases and the regularization parameter is decreased at an appropriate rate [see 22,46]. Additional details regarding the convergence properties of the embedding are provided in Appendix B.

FEATURES
We implemented SOCKS in Python, which has several available libraries for machine learning and reinforcement learning, such as Tensorflow [1], Keras [12], Scikit-Learn [34], PyTorch [33], and OpenAI Gym [10]. We utilize the Open AI Gym framework to be compatible with several existing libraries. This makes SOCKS comparable to several existing machine learning frameworks and promotes a more direct comparison with state-of-the-art machine learning and reinforcement learning algorithms.

Generating Samples
The algorithms in SOCKS are data-driven, which means they rely upon a sample of system observations S as in Assumption 2. Thus, we have implemented several sampling functions in SOCKS in order to generate samples from a system via simulation when a priori data is unavailable.
The process for generating samples consists of defining a sample generator, a function which generates a tuple contained within the sample S. For example, to generate a sample S = {( , , )} =1 as in (2) Here, the sample_generator function generates a single observation of the system, and the sample function computes a collection of = 100 observations taken from the system. SOCKS implements several commonly-used sample generators, including a one-step sample generator (shown above) and a trajectory generator, which generates samples of trajectories over multiple time steps of the form S = {( , , )} =1 , where ∈ X are the initial conditions, = { 1 , . . . , } is the sequence of states at each time step over the time horizon ∈ N, and = { 0 , . . . , −1 } is a sequence of control actions taken from a policy .

Stochastic Optimal Control
SOCKS can be used to solve the stochastic optimal control problem in (3). Given a sample S as in (2), we can approximate the integrals in (3) using an estimateˆ( , ) of the kernel distribution embedding ( , ), which can then be computed as Hilbert space inner products.
# Setup code omitted. policy = KernelControlFwd( cost_fn=cost_fn, constraint_fn=constraint_fn, ) policy.train(S, A) Here, KernelControlFwd class defines the algorithm, where we compute the optimal policy by minimizing the cost forward in time at each time step. The variables S and A define a sample S taken from the system as in (2) and a collection of admissible control actions in U, respectively. The cost_fn and constraint_fn are user-defined functions which return a real value. We can also solve the stochastic optimal control problem via dynamic programming (backward in time) by using KernelControlBwd in place of KernelControlFwd.

Stochastic Reachability
We can solve the terminal-hitting time stochastic reachability problem using SOCKS. Given a sample S as in (2), we can compute an empirical estimateˆ( , ) of ( , ), and (assuming the stochastic reachability value functions , = 1, . . . , , are in H ) we can approximate the stochastic reachability backward recursion by approximating the value function expectations in (6) via Hilbert space inner products with the estimateˆ( , ). In other words, we define the approximate value functions¯ * : X → [0, 1], = 0, . . . , , and form an approximation of the stochastic reachability backward recursion, given by, where * ( ) = 1 T ( ), and K, T ⊆ B X are the safe set and target set, respectively. Then as shown in [50,54], the solution to the approximate backward recursion,¯ * 0 ( 0 ), is an approximation of the maximal stochastic reachability safety probabilities. See [50] for more details.
# Setup code omitted. alg = KernelMaximalSR( time_horizon=time_horizon, constraint_tube=constraint_tube, target_tube=target_tube, problem="THT", ) alg.fit(S, A) Pr = alg.predict(T) Here, KernelMaximalSR is the stochastic reachability algorithm class, time_horizon is the number of time steps, S is a sample taken i.i.d. from the system, A is a collection of admissible control actions, T is a collection of test (or evaluation) points, i.e. the points where we seek to evaluate the safety probabilities, target_tube and constraint_tube are sets defining T and K, indexed by time, and "THT" specifies that we wish to solve the terminal-hitting time problem. In order to solve the first-hitting time problem, we simply replace "THT" with "FHT".
This means we can evaluate the safety probabilities for a system under a given policy, and enables an analysis of the likelihood of respecting a set of pre-defined safety constraints given by K, based in the same data-driven framework as the stochastic optimal controller synthesis in (33). The primary difference between our approach and existing tools such as [37], is that our approach is not based in reinforcement learning, and does not guard against unsafe exploration of the state space (while collecting the sample S), a well-known problem in safe RL [cf. 20].

Forward Reachability
We also implemented a forward reachable set estimator in SOCKS from [52]. Let P be some distribution on the state space X, and let S = { } =1 be a sample taken i.i.d. from P . The approximate forward reachable set classifierF is an estimate of the support of P and is computed as the solution to the following regularized least-squares problem: where > 0 is the regularization parameter, and is a separating kernel [see 14,52]. An RKHS H with kernel separates all subsets ⊂ X if there exists a function ∈ H such that for all ∉ , ( ) ≠ 0, and ( ′ ) = 0 for all ′ ∈ . The Abel kernel ( , ′ ) = exp(−∥ − ′ ∥ 2 / ), > 0, is a separating kernel, and is implemented in SOCKS. Note that a Gaussian RBF kernel is not a separating kernel, since constant functions are not included in a Gaussian RKHS [48].
The approximate forward reachable set is then given bỹ where is a threshold parameter, typically computed as = 1 − min 1≤ ≤˜( ), where ∈ S. The approximate forward reachable set classifier can accommodate non-convex regions, and the approximation converges almost surely to the true classifier. However, the approximation obtained via the algorithm is not a guaranteed under-or over-approximation, though it does admit finite sample bounds [14]. See [52] for more details.

Batch Processing
The primary computational hurdle of the kernel-based approach in SOCKS is the matrix inverse term in (17), which is O ( 3 ) in general, where is the sample size. Thus, the computation time scales polynomially as a function of the sample size. In addition, the optimal control algorithms frequently involve storing very large, dense matrices that scale as a function of , (the number of evaluation points) and (the number of admissible control actions). The large matrix sizes can lead to memory storage issues on systems with low available memory. In order to account for this, SOCKS implements a batch processing variant for algorithms with large sample sizes, which computes the solution in smaller "chunks". This does not affect the result, but leads to longer computation times, since we must compute multiple matrix multiplications rather than a single multiplication with a large matrix.

Dynamical System Modeling
OpenAI Gym currently implements several classical controls problems, including an inverted pendulum, a cart-pole system, and a "mountain car". These systems are contained within environments, which encapsulate the dynamics, constraints, and cost for the problem. Building on OpenAI Gym's standard framework, we have implemented a new type of learning environment, DynamicalSystem, which makes defining systems with dynamics easier. In addition, we have implemented several benchmark systems in SOCKS that involve classical controls problems which are not included in Ope-nAI Gym, including: (i) a satellite rendezvous and docking problem based on Clohessy-Wiltshire-Hill (CWH) dynamics, (ii) an -D stochastic chain of integrators, (iii) a nonholonomic vehicle, (iv) a point-mass system, (v) a benchmark quadrotor example [21], (vi) a planar quadrotor system, (vii) a translational oscillation with rotational actuation (TORA) system. We plan to add additional benchmarks, since OpenAI Gym can also be used to simulate hybrid dynamics, partially observable systems, and more.
Simulating a DynamicalSystem can be done easily. For example, we can evaluate the policy computed via the solution to the stochastic optimal control problem. Here, policy is the result of the stochastic optimal control algorithm. Behind the scenes, the simulation solves an initial value problem at each time step using Scipy's ODE solver. The input model is a zero-order hold. This means we can simulate continuous dynamics in discrete time with different sampling times without redefining or parameterizing the dynamics.

NUMERICAL EXPERIMENTS
All experiments were performed on an AWS cloud computing instance. The toolbox and code to reproduce all results and analysis is available at https://github.com/ajthor/socks.

Nonholonomic Vehicle
We consider a target tracking problem using nonholonomic vehicle dynamics, given by: where X ⊆ R 3 , U ⊆ R 2 , and we constrain the input such that 1 ∈ [0.1, 1], 2 ∈ [−10, 10]. We then discretize the dynamics in time with sampling time and apply an affine stochastic disturbance ∼ N (0, Σ), Σ = 0.01 . The goal is to minimize the distance to an object moving along the v-shaped trajectory shown in blue in Figure 1. We then collect a sample S = {( , , )} =1 , of size = 2500, and compute the optimal control actions using both the optimal control and dynamic programming algorithms in SOCKS with = 2. The results are shown in Figure 1, and computation time took approximately 0.204 seconds for the optimal control algorithm and 12.003 seconds for the dynamic programming algorithm. We can see that the system more The goal is to reach a small region close to the origin (green square) while remaining within a line of sight cone (blue triangle). The trajectory of the system using the approximate stochastic optimal control algorithm is shown in blue.
closely meets the terminal constraint using the dynamic programming algorithm, but the computation time increases dramatically, since we must compute a sequence of value functions.

Satellite Rendezvous and Docking
We consider an example of spacecraft rendezvous and docking, in which one spacecraft must dock with another while remaining within a line of sight cone. The Clohessy-Wiltshire-Hill dynamics are given by, , and parameters , . From [28], the dynamics in (24) can be written as a discrete- We apply the stochastic optimal control algorithm in SOCKS to the CWH system using the kernel bandwidth parameter = 0.

Double Integrator System
We consider a stochastic double integrator system in order to showcase the stochastic reachability analysis. The dynamics of the system with sampling time are given by, where X ⊆ R 2 , U ⊂ R, and ∼ N (0, Σ), Σ = 0.01 , is a random variable with a Gaussian distribution. We collect a sample S = {( , , )} =1 , = 2,500, taken i.i.d. from , a representation of the dynamics in (25) as a stochastic kernel, such that and are sampled uniformly from X and U, respectively, with in the range [−1.1, 1.1] 2 and in the range [−1, 1], and draw from (· | , ). The safe set is defined as K = [−1, 1] 2 and the target set is defined as T = [−0.5, 0.5] 2 . We then computed the stochastic reachability safety probabilities for both the terminal-hitting time problem and the first-hitting time problem at = 10,000 evaluation points using SOCKS and validated the result using Monte-Carlo. The computation time was ≈ 3 seconds for both problems. The result is shown in Figure 3.

Forward Reachability
We then demonstrate the forward reachable set algorithm in SOCKS for a translational oscillations by a rotational actuator (TORA) system [52]. The dynamics of the system are given by, We presume an initial distribution that is uniform over the region

Scalability & Computation Time
We now present a brief discussion of the scalability and computational complexity of the algorithms. As shown in [50,51,53], the sample size used to compute the empirical distribution embed-dingˆ( , ) presents the most significant computational burden, and is generally O ( 3 ) due to the presence of the matrix inversion in (16). We demonstrate this empirically for the algorithms presented in SOCKS in Figure 5. We calculated the computation times for the algorithms over 16 runs for different values of , and computed the statistical average and the 95% confidence interval. The black dots indicate the empirically measured times, and the blue bars indicate the confidence interval for our algorithm. We can see that the computation times scale polynomially with the sample size , as expected. This can be prohibitive, since the quality of the kernel-based approximation improves as the sample size tends to infinity. Nevertheless, several approximative speedup techniques (e.g. [22,35]) have been developed to alleviate the computational burden, and have been shown to reduce the computational complexity to O ( log ). These techniques are not currently implemented in SOCKS, but we plan to include them as part of a future release.
As mentioned in [50], the complexity of the kernel-based algorithms scales roughly linearly with the dimensionality of the system. This is primarily due to the fact that the system dimensionality only plays a role in the kernel evaluations, and does not significantly affect the computation of the empirical embeddingˆ( , ). We demonstrate this empirically for an -dimensional stochastic chain of integrators system as in (25) (see [50]), where we choose a fixed sample size = 1000 and vary the system dimensionality from = 50 to = 1000. The computation times are shown in Figure 6. We can see in Figure 6 that as the system is dimensionality is increased, the computation time increases roughly linearly. However, as mentioned in [50], for high-dimensional systems, the sample size needed to fully characterize the dynamics and the uncertainty increases as the system dimensionality increases, which can be prohibitive for the reasons mentioned above.

CONCLUSION & FUTURE WORK
In this paper, we introduced SOCKS, a toolbox for approximate stochastic optimal control and approximate stochastic reachability based on a data-driven statistical learning technique known as kernel embeddings of distributions. We have demonstrated the capabilities of the toolbox on a simple stochastic chain of integrators example, a more realistic satellite rendezvous and docking example, a target tracking scenario using nonlinear nonholonomic vehicle dynamics, and on a forward reachable set estimation problem. The approaches used in SOCKS are scalable, computationally efficient, and model-free.
We plan to introduce additional kernel-based algorithms and features to SOCKS, such as the capability to handle neural network reachability analysis, trajectory optimization, and chanceconstrained optimization.

ACKNOWLEDGMENTS
This material is based upon work supported by the National Science Foundation under NSF Grant Numbers CNS-1836900 and CMMI-2105631. 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. The NASA University Leadership initiative (Grant #80NSSC20M0163) provided funds to assist the authors with their research, but this article solely reflects the opinions and conclusions of its authors and not any NASA entity.

A STOCHASTIC OPTIMAL CONTROL
In this section, we give an overview of the stochastic optimal control algorithm. Additional details are provided in [51]. Given a sample S as in (2), taken i.i.d. from , we can compute an empirical esti-mateˆ( , ) of the conditional distribution embedding ( , ). By assumption 3, we assume that the objective and constraints can be decomposed as ( , ) = ( ) + ( ). Assuming the objective function 0 and constraints , = 1, . . . , , are elements of the RKHS H , we can approximate the expectation with respect to (· | , ) for all , = 0, 1, . . . , . Then, following [51], we propose the following form for the approximation of the policy. Given a set A ⊂ U of admissible control actions, A = {˜} =1 , ∈ N, we havê We can write (29) using the policy approximation in (30), and obtain where Υ is a feature vector with elements Υ = (˜, ·). We use this representation in the optimal control problem in order to approximate the objective and constraints. The following problem is approximately equivalent to the optimal control problem in (3), Furthermore, as the number of observations in the sample S and the number of admissible control actions in A tends to infinity, the solution to the approximate problem converges in probability to the true solution [51]. See Appendix B for a more detailed discussion of convergence. However, the optimization problem is unbounded below. In order to ensure feasibility, following [51], we add an additional constraint such that ( ) lies in the probability simplex ( ) ∈ { | 1 ⊤ = 1, 0 ⪯ }, where 1 is a vector of all ones. Thus, the approximate optimal control problem becomes This is a linear program, and can be solved efficiently, e.g. via several commonly-used interior point or simplex algorithms [9].
Additionally, this representation can be used to solve a backward-intime stochastic optimal control problem (dynamic programming). See [51] for more details.

B STABILITY & CONVERGENCE
We now seek to characterize the quality of the approximation and the conditions for its convergence. The convergence properties of kernel distribution embeddings are well-studied in literature. See e.g. [22,32,45,46] for more information. However, the nuances of the different convergence results for kernel distribution embeddings means that the results do not always generalize well to all problems under different kernel choices. For instance, the result in [46,Theorem 6] shows that the empirical estimateˆ( , ) converges in probability to the true embedding ( , ) at a rate of O (( ) −1/2 + 1/2 ), where is the regularization parameter in (14) and is the sample size, but [46] assumes that the RKHS is finite-dimensional, which does not hold for common kernel choices such as the Gaussian RBF kernel. Thus, we present convergence results for the approximate stochastic optimal control problem based in the theory of algorithmic stability [8]. Our result is close to the result presented in [32].
We first seek to characterize the convergence of the estimatê ( , ) in (16) to its actual counterpart ( , ) in (10). For simplicity of notation, we define the operator : X → R for all ∈ X via ( ′ ) = ( , ′ ).
Recall that the estimateˆis in a vector-valued RKHS Q and is the solution to the regularized least-squares problem in (14). Let : H × H → R be a real-valued cost function, defined by Let Z = X × U × X. We define the loss function : Q × Z → R, given by (ˆ, ( , , )) = ( ,ˆ( , )).
The risk, denoted by (ˆ), measures the expected loss (error) of the solutionˆto the regularized least-squares learning problem, and is defined as However, we cannot compute the risk directly since is unknown by Assumption 1. Thus, we seek to bound the risk by its empirical counterpart. Given a sample S ∈ Z as in (2), the empirical risk, denoted by S (ˆ), also known as the empirical error, measures the actual loss of the learning problem, and is defined as We useˆto denote the solution to the regularized least squares problem in (14), and letˆ\ denote the solution when a single observation is removed from S and letˆdenote the solution when the th observation is changed. We useˆ\ andˆin the following to assess the stability of the learning algorithm under minor changes to the sample S used to construct the estimateˆ.
We present the following definition, modified from [8], which allows us to characterize the stability of the learning algorithm with respect to the regularized least-squares problem in (14).
Definition 7 ( -admissible, [8,Definition 19]). A loss function on Q × Z is -admissible with respect to Q if the associated cost function is convex with respect to its first argument and the following condition holds, for all ′ ∈ H and 1 , 2 ∈ D, where is the domain of the first argument of .
We now seek to verify that the loss function pertaining to the regularized least-squares problem is -admissible with respect to Q. To this aim, we present the following proposition. is -admissible with respect to Q.
The proof follows [8, Lemma 20], which shows that the loss function using a Hilbert space norm is -admissible, where depends on the choice of kernel and the corresponding Hilbert space of functions H .
We now present the following definition from [8], modified to our particular formulation, which bounds the maximum difference in the loss function under minor variations to the sample S.
As a consequence of the above definitions, [8] shows that the regularized least-squares problem in the scalar RKHS case has uniform stability. We modify [8,Theorem 22] to a vector-valued RKHS in the following theorem. Theorem 1. Let H be an RKHS with kernel and Q be a vectorvalued RKHS of functions on X × U mapping to H . Let be bounded by < ∞, and let be a -admissible loss function with respect to Q. Then the learning algorithm given bŷ = arg min has uniform stability with respect to with We can ensure boundedness of the kernel of Q using the principle of uniform boundedness (also known as the Banach-Steinhaus theorem), since the kernel is bounded by . Then the proof follows directly from [8,Theorem 22].
We use this result to show that the regularized least-squares problem in (14) has uniform stability with respect to .
Theorem 2 ([8, Theorem 12]). Let be an algorithm with uniform stability with respect to a loss function such that 0 ≤ (ˆ, ( , , )) ≤ , for all ( , , ) ∈ Z and all sets S. Then for any ≥ 1 and any ∈ (0, 1) the following bounds hold with probability 1 − of the random draw of the sample S: (ˆ) ≤ S (ˆ) + 2 + (4 + ) Thus, using Theorem 2 with given by (45) from Theorem 1, we have that for any ≥ 1 and any ∈ (0, 1), with probability 1 − , the risk is bounded by: which shows that as the sample size increases, the empirical embeddingˆ( , ) in (16) converges in probability to the true embedding ( , ) in (10). Thus, the approximation of the expectations in (33) converge in probability to the true expectations, and the approximate optimization problems computed using the estimatê ( , ) converge to the true optimization problems as increases. Similarly, the approximate policyˆ( ) in (30) has the form of an empirical conditional distribution embedding as in (16), which suggests that the approximate policy (and consequently the approximately optimal control action) obtained via (33) also converges in probability as the number of admissible control actions in (30) is increased.