Motivated by the inherently uncertain nature of real-world systems, significant amounts of recent research has aimed at developing motion planning algorithms that explicitly reason about uncertainty in system dynamics, often with uncertainty modeled probabilistically. When probabilistic uncertainty is introduced, chance-constraints are a natural way to encode a desired level of safety in the motion planning problem. Historically, algorithms for chance-constrained motion planning under dynamics uncertainty have been restricted to the case of linear dynamics and additive Gaussian uncertainty [CCRRT, LQG1, Larss2, aoude2013probabilistically]
. This assumption is popular in part because linear transformations of multivariate Gaussians are still multivariate Gaussian, thus allowing for cumulative density functions to be easily encoded. However, such representations can be limiting and inaccurate as most real world systems are nonlinear and many sources of uncertainty are non-Gaussian. For one, almost all distributions propagated through nonlinear dynamics will produce a non-Gaussian state distribution. To handle nonlinear systems and non-Gaussian uncertainty, some prior works adopt Monte Carlo style approaches and sample the sources of uncertainty to fix the random variables at each of the sample values, generating multiple deterministic dynamical systems. These multiple systems can then be used either as a large number of constraints in an optimization routine or as a sub-routine in a sampling-based motion planner[calafiore2006scenario, Larss1, melchior2007particle]. Such approaches, however, usually become computationally intractable as a relatively large number of samples is needed to certify chance-constraints.
To avoid the computational intractability of sampling, many recent works adopt approaches that can establish upper bounds on state constraint violation using statistical moments of random variables. For example, higher order statistical moments of a random variable can be used in a sums-of-squares (SOS) program to find upper bounds on risk that are optimal w.r.t. the supplied moment information [jasour2018moment, jasour2019risk, bertsimas2005optimal]
. When only low order moments, such as mean and covariance, are used, concentration inequalities can be used to establish simple analytic upper bounds on the probability of state constraints being violated[wang2020non, hewing2019cautious].
To enforce chance-constraints on state, moments of the state distribution need to be known, but it is challenging to determine moments of state when random variables are propagated through nonlinear dynamics. Due to this difficulty, most prior works that utilize moment-based approaches resort to linearizing the dynamics to perform moment propagation [calafiore2006distributionally, dr_rrt, hewing2019cautious, carron2019data]. To address the problem of moment propagation for Gaussian driven stochastic processes, polynomial chaos expansions (PCE) are a well-studied method in the uncertainty quantification community to determine statistics of random variables propagated through some nonlinear function [xiu2010numerical]. A few efforts have been made towards using PCE in control and robotics applications [fisher2011optimal, nakka2019trajectory, kim2012generalized], but its practicality thus far has been limited by: 1) lack of guarantees on approximation error 2) the large number of terms that are often needed 3) the necessity of calling expensive quadrature routines to determine expansion coefficients.
In this paper, we address the problem of moment propagation through nonlinear systems by developing a framework that can search for a new dynamical system in terms of moment state. We present this framework in Section II and show that it can produce a polynomial system or even a linear system at the cost of increasing dimensionality. Given the statistical moments of the disturbance variables, these systems are in closed form and, thus, can be evaluated rapidly. In Section III, we extend the range of applicability of this framework to trigonometric polynomial systems; systems that have sines or cosines in the system state variables but can be transformed into polynomial systems by a change of variables. This is significant as it extends our method to encompass a wide range of systems important in control and robotics such as Dubin’s car. Finally, in Section IV, we provide an example application of this moment propagation framework to develop an improved distributionally robust RRT (DR-RRT) that uses a stochastic Dubin’s car model without relying on linearization as was previously required. In Section V, we provide examples of the accuracy and run-times of our approach compared to linearization and naive Monte Carlo. The source code for TreeRing is available at github.com/wangmengyu96/TreeRing.
Ii Moment Propagation for Polynomial Systems
This section presents the methodology for finding a moment state dynamical system. Throughout this section, we will be concerned with the system:
Where is the
dimensional state random vector,is a dimensional disturbance random vector that is assumed to be independent of , and is a vector-valued function that is polynomial in and . This formalism implicitly models control inputs as deterministic control variables can be represented in as random variables that take on a value with probability one. We assume that the moments of
are known for each time step; we argue this is not a particularly restrictive requirement as moments can be computed via auto-differentiation of the moment generating function of known distributions. While we restrict our focus to discrete time systems in this paper, we would like to note that there are standard discrete time approximations of continuous time stochastic differential equations, such as the Euler-Maruyama and stochastic Runge-Kutta methods, that are analogous to the well known ones for deterministic systems[kloeden2013numerical].
II-A introduces multi-index notation which will be used frequently throughout this paper. II-B then establishes several properties of polynomial systems relevant to moment propagation. In II-C and II-D, we show that a moment state dynamical system can be found if we find a complete moment basis for the set of moments we are interested in. Finally, II-E presents an algorithm for finding complete moment bases, which we call TreeRing.
Ii-a Multi-Index Notation
For simplicity, we adopt multi-index notation. That is, for any dimensional vector and multi-index , we adopt the following short-hand:
All multi-indices are denoted with lower case Greek letters. This provides much simpler expressions for polynomials and for moments of random vectors. For example, if we want to express the moment in multi-index notation, we can simply let and write . One of the advantages of this formalism is the existence of a one-to-one mapping between moments of and multi-indices in , and similarly for , so it is equivalent to talk about one or the other.
Ii-B Properties of Polynomial Systems
We begin by reviewing a few elementary properties of polynomials arising from their ring structure that are relevant to Proposition 3. In Proposition 1 and 2, denotes the ring of polynomials in the vector over the field . In this paper, we will only work with and , but we adopt the more general notation for concision in the propositions. Proposition 1 states a few closure properties; these are of importance because simulating the discrete-time system essentially consists of applying these operations.
(Closure Properties) For any
A useful heuristic for evaluating the complexity of a polynomial expression is its degree. In general, it’s undesirable for the degree of polynomials to increase as the application of operations to high degree polynomials is likely to result in a polynomial with a large number of terms. Proposition2 states the degree of the resulting polynomial after applying the closure operations.
In the case of vector-valued polynomials, the degree operator returns the sum of the vector of degrees of each entry of the vector-valued polynomial which we will refer to as the degree vector. To compute its degree after applying an operation, simply apply the above proposition to each of the scalar-valued entries and sum the resulting values.
Closure under exponentiation is a particularly important property as it implies that raising any polynomial in a random vector to the power of results in a polynomial . By applying the linearity of expectation to , we see that it can be computed in terms of moments of . By simply extending this idea to the case where is a vector-valued function and is raised to a multi-index, we arrive at Proposition 3a which is a simple extension of Proposition 1 in [wang2020non] from the scalar-valued to the vector-valued case. Proposition 3b then states the degree of the resulting polynomial.
a) For any dimensional random vector , any vector-valued polynomial s.t. with degree vector and any , there exists a set of multi-indices, , and coefficients s.t:
While is a vector-valued polynomial in , is a scalar-valued polynomial in Proposition 3.
Ii-C Moment Update Forms
If we treat the state and disturbance vectors as one vector in our polynomial system, we can apply Proposition 3a and see that any moment of
can be computed in terms of moments of the joint distribution ofand . Since we assumed that the disturbance is independent of state, we have that moments of can be expressed in the form of equation 4, the moment update form, by splitting the multi-index s.t. . Note that by Proposition 3b, the polynomial in the moment update form is of relatively low order: to be precise.
Given the polynomial system, equation 1, and some multi-index corresponding to a moment of the state vector, the moment update form (MUF) is:
Throughout this paper, we will denote the multi-index set of the moment update form induced by as and denote the sets of the state and disturbance components respectively by and .
In practice, moment update forms can be easily found with standard symbolic algebra packages. In this paper, we used the the built-in functionality in SymPy to express symbolic expressions as polynomials in a set of variables [meurer2017sympy].
The moment update form essentially states that moments of the state vector at time are linear in the moments of the state vector at time when given moments of external disturbances. Now if there is independence between elements of the state vector, may be factored into the product of lower order moments to arrive at what we will refer to as the reduced moment update form. There are ultimately interesting tradeoffs between using the reduced and non-reduced moment update forms that we will discuss in the following subsection.
Given the polynomial system, equation 1, and some multi-index corresponding to a moment of the state vector, the reduced moment update form is:
Where , the expression involving the product operator equals and for each , can not be factored further.
The reduced moment update form requires the additional step of factorizing , but this can be performed with standard search algorithms to find connected components on the, manually specified, dependence graph of .
The graph is the dependence graph of where is the set of variables in and is the set of undirected edges s.t. i.f.f. and are dependent.
However, pairwise independence of a collection of random variables does not necessarily imply mutual independence. This is a difficult issue to resolve in general, so we make the following assumption on pairwise independence of state variables.
Suppose and are random vectors that are subsets of . Then pairwise independence between elements of and implies and are independent. That is, letting denote probability density functions:
denote probability density functions:
Note that moments of joint distributions can be factored in the same way as probability density functions above. With these facts in mind, the goal now is to find a partition of into multiple random vectors s.t. equation 6 can be applied. It should be easy to see that the set of connected components of the dependence graph forms the correct partition because by definition, pairwise independence corresponds to the existence of an edge. With the example in Fig. 1, we can perform the following factorization:
Thus, for any given multi-index , we can find the sub-graph of the dependence graph containing the variables with non-zero entries in and then find a factorization by finding the connected components of the resulting sub-graph.
Ii-D Moment Bases
Since we assume the moments of the external disturbances are known, the (possibly reduced) moment update form seems to suggest there may be a way to recursively determine the desired moments of the state vector across a finite time horizon when given some initial condition. In fact, this is possible if we find a complete moment basis. Suppose we are interested in determining a set of moments of the state vector ; we will denote the set of multi-indices corresponding to these moments by . Such sets of multi-indices will be referred to as moment bases. Given a moment basis, we denote the moment state of the system by:
This notation also lends itself to expressing the set of all moments of that are needed. Recall that for any given , denotes the set of disturbance multi-indices in the (possibly reduced) MUF, so for each we can denote the set of disturbances needed for that particular by:
By taking the union across , we arrive at the set of disturbance moments needed. For brevity, we will denote:
We are now ready for the definitions of complete moment bases w.r.t. the (possibly reduced) MUF. When a moment basis is not complete, we will want to find a completion.
is a complete moment basis w.r.t. the moment update form if , the multi-index set in the moment update form corresponding to the state vector, , is a subset of .
is a complete moment basis w.r.t. the reduced moment update form if , , .
A completion of a moment basis denoted is a super-set of that is a complete moment basis w.r.t. the (possibly reduced) moment update form.
In both cases, completeness of the moment bases is a sufficient condition that the (possibly reduced) MUFs of the elements in the moment basis form a moment state dynamical system. This should be clear as completeness of implies that every moment in can be expressed explicitly in terms of and moments of .
A moment state dynamical system given a complete moment basis is:
Where is a vector-valued function formed by the (possibly reduced) moment update forms of .
Ultimately, both Definition 4 and 5 simply try to capture the idea that every moment at time step can be computed in terms of moments in the moment state at the previous time step. However, the distinction is important because a complete moment basis w.r.t. the un-reduced MUF has linear time-varying dynamics. This can be easily seen by inspection of equation 4 because moments of are assumed to be known and thus can be treated as constants. However, empirically, we find that using the reduced MUF has the advantage of producing smaller completions. It also makes intuitive sense that using reduced MUFs will result in smaller completions as the reduced MUF tries to express moments in terms of the lowest order moments possible.
Ii-E TreeRing: An Algorithm for Finding Completions of Moment Bases
If a moment basis is incomplete, a simple step to take would be to simply expand the moment basis by adding the multi-indices that are not in to and determining their MUFs. Recursively doing so produces an algorithm analogous to tree search, where the nodes are moments of the state vector that show up in MUFs, and the nodes are expanded if the moment is not in the current moment basis. We call this algorithm TreeRing as it is essentially tree search, but over the ring of polynomials. TreeRing in Algorithm 1 uses the reduced MUF, but an algorithm using the un-reduced MUF can be arrived at by simply removing line number 5, letting , and replacing with everywhere. If the algorithm terminates, then we have that the resulting is a complete moment basis. By finding all of the (possibly reduced) MUFs associated with moments in the complete moment basis, we arrive at our moment state dynamical system.
Iii Extension to Trigonometric Polynomial Systems
Many dynamical systems in control and robotics are not polynomial, but can be made polynomial by a change of variables replacing sines and cosines with new indeterminants. We will refer to such systems as trigonometric polynomial systems. For example, consider the following Dubin’s car model with uncertainty in actuation (note that the control inputs are absorbed as constant shifts in and for brevity):
We can transform the above system into a polynomial system by making the substitutions and .
Update relations for and were arrived at by using the trigonometric sums formulas to expand out and into the expressions shown in equation 13 where and . To apply the methods developed for polynomial systems, however, we need to be able to compute the moments of and . In III-A, we show that these trigonometric moments
can be computed directly in terms of the characteristic functions of the original random variables. This retains the generality of our results on polynomial systems as extensive catalogues of characteristic functions for distributions are available.
Iii-a Computing Trigonometric Moments
In this subsection, we show moments of the form:
for any random variable , can be computed in terms of the characteristic function of , defined as:
The characteristic function is simply the Fourier transform of the PDF of, so it makes intuitive sense that there exists deep connections between it and trigonometric functions. We provide explicit closed form solutions for the pure moment case, i.e when one of or equals zero, and provide an algorithmic approach for finding an expression using symbolic algebra in the case when both . The first major insight is that moments of and may be related to the characteristic function via Euler’s identity:
And so we immediately have a relation that can be used to compute the first trigonometric moments:
The trigonometric power formulas can then be used to relate quantities of the form to for any and similarly for the sin function. Applying the expectation operator to the trigonometric power formulas, applying the linearity of expectation and substituting in equation 17, we have for some coefficients that equations 18 hold. The exact values of the coefficients can be found in [trig_power]; they are omitted here for concision.
In the general case when , we can substitute in the exponential forms of and to get:
Observe that the expression:
can be written as a polynomial in and by the closure property of . An alternative view is that such an expression would be a Laurent polynomial in (Laurent polynomials are generalizations of polynomials that allow for negative powers). In either case, since the following equality holds :
By expanding equation 19b into a (possibly Laurent) polynomial and applying the linearity of expectation, we arrive at the weighted sum of evaluated at multiple values. Like before, the (possibly Laurent) polynomial can be found using symbolic algebra.
Iv Example Application: DR-RRT with Stochastic Dubin’s Dynamics
As an example application of this moment propagation framework, we present a version of distributionally robust RRT (DR-RRT) with uncertainty propagated through nonlinear Dubin’s car dynamics. We make two major improvements to the original DR-RRT formulation: 1) we propagate moments through nonlinear dynamics as opposed to linear systems in the original formulation and 2) we use a different method to assess risk that is strictly less conservative [dr_rrt]. Throughout this section, the robot is modelled with the stochastic Dubin’s car described in equation 13 so denotes the state of the robot and denotes its position. In this section, for any random vector , will denote its mean vector and covariance matrix. Algorithm 2 describes this nonlinear DR-RRT algorithm; the key points of interest are the StochasticSteer and RiskBound functions which are described in the following subsections.
Iv-a Stochastic Steering
It is common in RRTs, especially for ground vehicles, to perform steering between nodes with Dubins paths [karaman2011anytime, dubins1957curves]. In our formulation, a deterministic Dubin’s path is calculated, then control inputs for steering, , to track the Dubin’s path at a constant speed are computed with the discretized deterministic Dubin’s dynamics. The steering disturbance in our stochastic system given by Eq 13, is then offset by to account for both the control and uncertainty in a single term. In an offline step, we run TreeRing on the system in equation 13 with an initial moment basis, , containing the first and second moments of and :
And we arrive at the following completion w.r.t. the reduced MUF :
We use this completion as the completion w.r.t the un-reduced MUF also has the moments:
requiring more equations in the system. The reduced MUF for each is then found using symbolic algebra to arrive at our moment state dynamical system. In this way, we can compute a sequence of mean vectors and covariance matrices to steer the vehicle from an existing node to a new one. To see the full set of equations associated with these moment bases, see the example provided in the Github repository.
Iv-B Bounding Risk
In the environment, polytopic obstacles are represented by the intersection of half-spaces defined by a collection of linear inequalities:
We are interested in ensuring that across a time step horizon, the probability of collision with an obstacle is less than some . This risk can be upper bounded as such by applying Boole’s Inequaltiy twice:
So it is sufficient to establish upper bounds on for all times and for all obstacles. Since the probability of the intersection of a collection of events is less than the probability of any individual event occurring, for the obstacle, we have that:
That is, the probability of the robot being in the obstacle is upper bounded by the probability of violating any one of the linear inequalities that define the obstacle; so we can upper-bound the probability of being in the obstacle with an upper-bound on the probability of being in the least risky half-space that defines it. Departing from previous works, we directly find an upper bound on risk as opposed to checking a sufficient condition (more discussion in IV-C) [calafiore2006distributionally, dr_rrt]. By Cantelli’s inequality, for any measurable function , whenever :
Note that the restriction that is not particularly restrictive as it is only not satisfied when the average case corresponds to collision. Thus, we simply let the risk bound be whenever . In our case, we have that for each linear inequality. The mean and covariance can be computed in terms of and with:
Note that since the Cantelli bound holds for any measurable function
, in principle, we can extend this formulation to the case of non-convex obstacles described as the sub-level sets of polynomials. However, since we need the mean and variance of, doing so would require higher moments of , so more moments need to be propagated.
Iv-C Comparison with Prior DR Constraints
In prior works [calafiore2006distributionally, dr_rrt], given a chance-constraint of across the whole horizon, a maximum allowable risk for each time step and each obstacle is assigned s.t. . The following sufficient condition for the chance-constraint to be satisfied is then checked for each linear inequality:
And collision is defined as the above condition not being satisfied. We would like to note that this sufficient condition is ultimately equivalent to the bound we employ as simple manipulation of the equation results in:
Thus, our approach uses the same fundamental inequality. However, instead of allocating a priori, we determine upper bounds on risk and ensure that the summation is less than the overall chance-constraint thus allowing for less conservative results.
V Comparison with Monte Carlo and Linearization
This section uses the Dubin’s car model to compare the proposed method for exact nonlinear moment propagation against naive Monte Carlo and moment propagation with a linearized system. We linearize the Dubin’s car model about the initial state and, since we apply Euler integration, the matrices and vector are scaled by the time step to arrive at the following discrete time system:
denotes the identity matrix. The mean vector and covariance matrix thus have the following dynamics:
The distributions of the disturbances are chosen to be and . Fig. 3 compares the results of our nonlinear moment propagation method against naive Monte Carlo and the standard method of propagating a linearized system. We note that the Monte Carlo results closely match our moment propagation method, while the second moments of the linearized system diverge quite significantly from those obtained from the other two methods. The nonlinear moment propagation method is also very fast; in our C++ implementation, propagation over time steps takes only 35.2 ms for a computation time of microseconds per time step. This is on par with speeds for well-implemented moment propagation methods for linear systems and many orders of magnitude faster than naive Monte Carlo.
In this paper, we present a new framework for nonlinear moment propagation with wide applicability to problems in robotics and control, namely to trigonometric polynomial systems. The key result is an algorithm that can search for a moment state dynamical system that can be used to compute moments of the state vector given the moments of external disturbances. We demonstrate its application to a Dubin’s car model with actuation uncertainty which is then used in an improved distributionally robust RRT. Finally, numerical experiments show that our method is as fast as moment propagation techniques that employ linearization while being exact. As TreeRing
can produce relatively large systems of equations, making manual conversion to code both tedious and prone to error, future work should consider developing code generation capabilities. There are many other potential applications of this exact moment propagation technique including but not limited to Gaussian process motion planning, trajectory optimization, and nonlinear estimation which should be investigated further in future works.