Nonlinear Uncertainty Control with Iterative Covariance Steering

TL;DR

This paper introduces an iterative covariance steering (iCS) method for controlling nonlinear stochastic systems to reach a desired state distribution while satisfying probabilistic constraints, demonstrated on a noisy double integrator.

Contribution

The paper presents a novel iterative approach for nonlinear uncertainty control that linearizes the problem and solves it as a convex program, enabling effective distribution steering.

Findings

Successfully controls a nonlinear system with probabilistic constraints.

Demonstrates effectiveness on a double integrator with quadratic drag.

Provides a practical iterative algorithm for nonlinear stochastic control.

Abstract

This paper considers the problem of steering the state distribution of a nonlinear stochastic system from an initial Gaussian to a terminal distribution with a specified mean and covariance, subject to probabilistic path constraints. An algorithm is developed to solve this problem by iteratively solving an approximate linearized problem as a convex program. This method, which we call iterative covariance steering (iCS), is numerically demonstrated by controlling a double integrator with quadratic drag force subject to additive Brownian noise while satisfying probabilistic path constraints.

Figures

Equations

$$\mathrm{d}{x}_t=f(x_t,u_t,t)\mathrm{d}t+G_t\mathrm{d}{w}_t,t\in [t_0,t_f],$$

$$\mathbb{E}(x_{t_0})={\bar{x}}_0,\mathbb{E}({\tilde{x}}_{t_0}{\tilde{x}}_{t_0}^{\text{mbox}\mathsf{T}})=P_{x_0}.$$

$$\mathbb{P}(x_t\in {\mathcal{X}}_t)\ge 1-p_{x,t},\mathbb{P}(u_t\in {\mathcal{U}}_t)\ge 1-p_{u,t},$$

$$\mathbb{E}(x_{t_f})={\bar{x}}_f,\mathbb{E}({\tilde{x}}_{t_f}{\tilde{x}}_{t_f}^{\text{mbox}\mathsf{T}})=P_{x_f},$$

$$J(u)=\int_{t_{0}}^{t_{f}}\big{[}\ell(\bar{u}_{t},\bar{x}_{t})+\mathbb{E}(\tilde{x}_{t}^{\mbox{\tiny\sf T}}Q_{x,t}\tilde{x}_{t}+\tilde{u}_{t}^{\mbox{\tiny\sf T}}Q_{u,t}\tilde{u}_{t})\big{]}\mathrm{d}t.$$

$$\sigma :=t_f-t_0.$$

$$\mathrm{d}{x}_{\tau }=\sigma f(x_{\tau },u_{\tau },\tau )\mathrm{d}\tau +\sqrt{\sigma }{G}_{\tau }\mathrm{d}{w}_{\tau },\tau \in [0,1],$$

$$J(u)=\sigma\int_{0}^{1}\big{[}\ell(\bar{u}_{\tau},\bar{x}_{\tau})+\mathbb{E}(\tilde{x}_{\tau}^{\mbox{\tiny\sf T}}Q_{x,\tau}\tilde{x}_{\tau}+\tilde{u}_{\tau}^{\mbox{\tiny\sf T}}Q_{u,\tau}\tilde{u}_{\tau})\big{]}\mathrm{d}\tau.$$

$$\mathrm{d}{x}_{\tau }\approx (A_{\tau }^i{x}_{\tau }+B_{\tau }^i{u}_{\tau }+r_{\tau }^i)\mathrm{d}\tau +\sqrt{\sigma }{G}_{\tau }\mathrm{d}{w}_{\tau },$$

$$A_{\tau }^i:=\sigma {\frac{\partial f}{\partial x}|}_{({\hat{x}}_{\tau }^i,{\hat{u}}_{\tau }^i)},B_{\tau }^i:=\sigma {\frac{\partial f}{\partial u}|}_{({\hat{x}}_{\tau }^i,{\hat{u}}_{\tau }^i)},$$

$$r_{\tau }:=\sigma f({\hat{x}}_{\tau }^i,{\hat{u}}_{\tau }^i,\tau )-A_{\tau }^i{\hat{x}}_{\tau }^i-B_{\tau }^i{\hat{u}}_{\tau }^i.$$

$${\tau }_k:=\frac{k}{N},k\in {\mathbb{N}}_0^N.$$

$$u_{\tau }=u_k,\tau \in [{\tau }_k,{\tau }_{k+1}),k\in {\mathbb{N}}_0^{N-1}.$$

$$\begin{gathered} x_{k+1}={\Phi }^i({\tau }_{k+1},{\tau }_k)x_k \\ +{\int }_{{\tau }_k}^{{\tau }_{k+1}}{\Phi }^i({\tau }_{k+1},\tau )(B_{\tau }^i{u}_k+r_{\tau }^i)\mathrm{d}\tau \\ +\sqrt{\sigma }{\int }_{{\tau }_k}^{{\tau }_{k+1}}{\Phi }^i({\tau }_{k+1},\tau )G_{\tau }\mathrm{d}{w}_{\tau },k\in {\mathbb{N}}_0^{N-1}, \end{gathered}$$

$$\frac{\partial }{\partial \tau }{\Phi }^i(\tau ,s)=A_{\tau }^i{\Phi }^i(\tau ,s),{\Phi }^i(\tau ,\tau )=I.$$

$$x_{k+1}=A_k^i{x}_k+B_k^i{u}_k+r_k^i+\sqrt{\sigma }{G}_k^i{w}_k,$$

$$A_k^i$$

$$B_k^i$$

$$r_k^i$$

$$\Sigma ={\int }_{{\tau }_k}^{{\tau }_{k+1}}{\Phi }^i({\tau }_{k+1},\tau )G_{\tau }{G}_{\tau }^{\text{mbox}\mathsf{T}}{\Phi }^{i\text{mbox}\mathsf{T}}({\tau }_{k+1},\tau )\mathrm{d}\tau ,$$

$$\begin{array}{rlrl}& A_k^i:=A_{{\tau }_k}^i\mathrm{d}\tau +I,& B_k^i& :=B_{{\tau }_k}^i\mathrm{d}\tau ,\\ & G_k^i:=\sqrt{\mathrm{d}\tau }{G}_{{\tau }_k},& r_k^i& :=r_{{\tau }_k}^i\mathrm{d}\tau .\end{array}$$

$$A_{k_1,k_0}^i:=\{\begin{array}{ll}{A}_{k_1}^i{A}_{k_1-1}^i\cdots A_{k_0}^i,& k_1\ge k_0,\\ I,& k_1<k_0,\end{array}$$

$$B_{k_1,k_0}^i:=A_{k_1,k_0+1}^i{B}_{k_0},G_{k_1,k_0}^i:=A_{k_1,k_0+1}^i{G}_{k_0}^i,$$

$$r_{k_1,k_0}^i:=A_{k_1,k_0+1}^i{r}_{k_0}.$$

$$U_k$$

$$W_k$$

$$x_k={\bar{A}}_k^i{x}_0+{\bar{B}}_k^i{U}_{k-1}+{\bar{r}}_k^i{1}_k+\sqrt{\sigma }{\bar{G}}_k^i{W}_{k-1},$$

$${\bar{B}}_k^i$$

$${\bar{G}}_k^i$$

$${\bar{r}}_k^i$$

$$X={\mathcal{A}}^i{x}_0+{\mathcal{B}}^{i}U+R^i+\sqrt{\sigma }{\mathcal{G}}^{i}W.$$

$${\mathcal{Q}}_x:=\mathrm{blkdiag}(Q_{x,{\tau }_0},\dots ,Q_{x,{\tau }_{N-1}},0_{n_x}),$$

$${\mathcal{Q}}_u:=\mathrm{blkdiag}(Q_{u,{\tau }_0},\dots ,Q_{u,{\tau }_{N-1}}).$$

Full text

I INTRODUCTION

Guidance and control design has generally followed the standard approach where an open-loop reference optimal control is solved with respect to the nonlinear dynamics and then a feedback controller is subsequently designed with respect to the dynamics linearized about the reference trajectory. Hence, there is an implicit unidirectional dependence of the feedback controller on this reference trajectory, but there is no direct mechanism from which the reference trajectory is affected by the closed-loop system behavior. Intuitively, if we could explicitly couple the design of the reference trajectory with the design of the feedback controller, then since we are optimizing over a larger set we could improve closed-loop system performance. The situation becomes more complex with the introduction of state constraints and uncertainty. If the closed-loop statistics of a system are not considered, then the reference trajectory design must be conservative to satisfy the constraints. For systems that are significantly influenced by uncertain external forces, the conservatism of this approach may lead to greatly increased control cost or even infeasibility.

In this paper we consider the system state to be a random vector which evolves according to a nonlinear stochastic differential equation with additive Brownian noise. By letting the state to be a random vector, the control problem can be formulated as one of simultaneously steering each sample trajectory, and as a consequence, we can analytically study the difference between open and closed-loop control [1]. This machinery will serve as our mechanism to understand the coupling between the reference trajectory and the feedback controller. We assume that the state is normally distributed at the initial time, and we will design a control that steers the mean and the covariance of the initial state distribution to some given terminal values at the final time. This problem is referred to as the nonlinear covariance steering (CS) problem. Since the state is assumed to be normally distributed, and thus unbounded, we must treat state constraints probabilistically. That is, the probability that the constraints are satisfied must be greater than some prespecified value. Since, by construction, these constraints may not be met for every sample path, they are often referred to as chance constraints [2, 3].

The special case of linear time-varying stochastic systems with additive Brownian noise has been extensively studied in the literature. It has been shown that if the system is controllable, then the state covariance is also controllable [1]. That is, for an initial covariance $P_{x_{0}}>0$ at time $t_{0}$, there exist a state feedback gain defined on the interval $[t_{0},t_{f}]$ that steers the covariance to any final value $P_{x_{f}}>0$ for any time $t_{f}>t_{0}$. The solution to the optimal linear CS problem with expected quadratic cost was given by Chen et al. [4, 5, 6], and the solution was found to be closely related to the classical linear quadratic feedback control. The discrete linear CS problem with quadratic cost has also been studied and a similar close connection to linear quadratic control has been shown [7].

It is well known that, for linear systems, and in the absence of any constraints, the mean and the covariance have independent dynamics, and that the mean state is controlled by the mean open-loop control and the covariance is controlled by the state feedback gain. It follows that, without constraints, we can consider the mean steering and the covariance steering as separate problems but that, when there are constraints, the mean and covariance are coupled. In other words, for linear systems, the reference trajectory explicitly depends on the closed-loop behavior of the system when the state or control is constrained. For the discrete-time linear case with convex chance constraints, the chance-constrained CS problem can be cast as a deterministic convex optimization problem [8]. This work was later extended to include non-convex chance constraints using mixed-integer programming [9].

In this paper, we propose an algorithmic approach to solve the nonlinear CS problem by iteratively solving the linear CS problem with respect to the reference trajectory of the previous step. This algorithm, which we will refer to as iterative CS (iCS), is a natural extension of linear CS in the spirit of other well known successive approximation methods such as differential dynamic programming (DDP) [10] and iterative LQG (iLQG) [11], which both compute a feedback control by backwards propagating an approximation of the value function. For iCS, we similarly approximate the nonlinear dynamics about a reference trajectory, but, in contrast, the control updates are found by solving a convex program, which has the benefit of allowing direct consideration of probabilistic constraints at the cost of computation time and restricts the approximation of the dynamics to first order.

To the best of our knowledge, there are currently no known methods to solve the nonlinear CS problem. Furthermore, in contrast to the existing literature on chance constrained linear CS, we begin our analysis with a continuous stochastic system and describe an exact discretization procedure.

I-A Notation and Preliminaries

For a random vector $z$ on a probability space $(\Omega,\mathscr{F},\mathbb{P})$, we denote the expectation of a function $f$ of $z$ as $\mathbb{E}[f(z)]$. The mean value of $z$ is denoted by $\bar{z}:=\mathbb{E}(z)$, and the difference from the mean as $\tilde{z}:=z-\mathbb{E}(z)$. The covariance of $z$ is denoted by $P_{z}:=\mathbb{E}(\tilde{z}\tilde{z}^{\mbox{\tiny\sf T}})$. The complement of an event $A\subseteq\Omega$ is denoted $A^{c}=\Omega\setminus A$, and we use the shorthand $\{\omega\in\Omega:z(\omega)\in B\}=\{z\in B\}$ to denote events. Dependence of a quantity $y$ on time $t$ is denoted by $y_{t}$. For a square matrix $A$, we write $A>0$ $(\geq 0)$ if $A$ is positive (semi-)definite, i.e., $x^{\mbox{\tiny\sf T}}Ax>0$ $(\geq 0)$ for all nonzero real vectors $x$. The set of natural numbers, including zero, is written as $\mathbb{N}_{0}$, and $\mathbb{N}_{+}=\mathbb{N}_{0}\setminus\{0\}$. We will denote by $\mathbb{N}_{0}^{m}$ the set of natural numbers up to, and including, a positive integer $m$ (similarly for $\mathbb{N}_{+}^{m}$).

II PROBLEM FORMULATION

Consider the nonlinear stochastic differential equation

[TABLE]

where $x_{t}\in\mathbb{R}^{n_{x}}$ is the state, $u_{t}\in\mathbb{R}^{n_{u}}$ is the control input, and $w_{t}$ is an $n_{w}$-dimensional standard Brownian motion. At time $t_{0}$, the state $x_{t_{0}}$ is assumed to be normally distributed with fixed mean and covariance

[TABLE]

At each time, the state and control are constrained in probability to given convex sets

[TABLE]

where $0<p_{x,t}<1/2$ and $0<p_{u,t}<1/2$ are prescribed maximum probabilities of failure. The constraints (3) are referred to as chance constraints. We wish to find a control that brings the state $x_{t}$ to a final distribution at time $t_{f}$ with given mean and covariance

[TABLE]

where $P_{x_{f}}$ is a given positive-definite matrix, while minimizing the cost functional

[TABLE]

Here $Q_{x,t}>0$ and $Q_{u,t}\geq 0$ are weight matrices, and $\ell$ is an integrable function that is convex in $\bar{u}_{t}$ and $\bar{x}_{t}$, and the optimization is performed over the control $u$. In summary, we are interested in solving the following problem.

Problem 1

Nonlinear Covariance Steering. Find a control $u_{t}^{*}$ to minimize the cost (5) subject to the dynamics (1), terminal state constraints (4), and chance constraints (3).

In the remainder of this section, we will develop a linear approximation of (1) in the neighborhood of a given reference. Then, after discretizing the linearized system, we will focus our analysis on the discrete linear system.

II-A Time Normalization and Linearization

We begin by normalizing the time domain $[t_{0},t_{f}]$ to the unit interval using the dilation coefficient [12]

[TABLE]

Let $\tau\vcentcolon=(t-t_{0})/\sigma\in[0,1]$ be the normalized time, from which it follows that $\sigma=\mathrm{d}t/\mathrm{d}\tau$. Since the Brownian motion increment $\mathrm{d}w_{t}$ has variance $\mathrm{d}t$, we scale the diffusion term in (1) by $\sqrt{\sigma}$ to obtain $\mathrm{d}w_{\tau}$ with variance $\mathrm{d}\tau$ (i.e., $\mathrm{d}w_{t}$ is identically distributed with $\sqrt{\sigma}\mathrm{d}w_{\tau}$). The time normalized system is then given by

[TABLE]

and the time normalized cost is given by

[TABLE]

Next, we linearize (7) about a given reference trajectory $(\hat{x}^{i}_{\tau},\hat{u}^{i}_{\tau})$, where $i\geq 1$ is an index to count successive linearizations. This procedure results in the linear stochastic system

[TABLE]

where

[TABLE]

[TABLE]

II-B Discrete Approximation

Let $0=\tau_{0}<\tau_{1}<\dots<\tau_{N}=1$ be a partition of the interval $[0,1]$, where

[TABLE]

Henceforth, we will write $x_{k}\vcentcolon=x_{\tau_{k}}$ and $u_{k}\vcentcolon=u_{\tau_{k}}$. We use a zero-order-hold (ZOH) discretization of the control given by

[TABLE]

Substituting $u_{k}$ in (9), we obtain the solution [13]

[TABLE]

where $\Phi^{i}(\tau,s)$ is the state transition matrix for system (9), which satisfies

[TABLE]

For $k\in\mathbb{N}_{0}^{N-1}$, we rewrite (14) as

[TABLE]

where $w_{k}\in\mathbb{R}^{n_{x}}$ are independent and identically distributed $\mathcal{N}(0,I)$ and where

[TABLE]

The stochastic integral in (14) is a zero-mean Gaussian random vector with covariance

[TABLE]

and therefore the matrix $G^{i}_{k}$ is selected such that $G^{i}_{k}w_{k}\in\mathbb{R}^{n_{x}}$ is distributed $\mathcal{N}(0,\Sigma)$. Since $\Sigma$ may have rank greater than $n_{w}$ (e.g., $\Sigma$ is full rank if $(A^{i}_{\tau},G_{\tau})$ is controllable), the matrix $G^{i}_{k}$ must be $n_{x}\times n_{x}$. Furthermore, since the solution to $G^{i}_{k}G^{i\mbox{\tiny\sf T}}_{k}=\Sigma$ is not unique, there exist a continuum of coefficient matrices $G^{i}_{k}$ such that the resulting state processes are equidistributed.

The discrete formulation (14) is an exact representation of the linear system (9) with ZOH control. However, since the previous integrals may be difficult to compute, the following first-order approximation is commonly used:

[TABLE]

We remark that the method of discretization chosen does not affect any of the following discussion.

Following [7, 8, 9], we will rewrite the discrete system (16) as a single linear equation. Define the state transition matrix from step $k_{0}$ to step $k_{1}$ as

[TABLE]

and the corresponding transitions for the control, noise, and affine terms as

[TABLE]

[TABLE]

Define concatenated control and disturbance vectors at step $k$ as

[TABLE]

Then the state at step $k$ can be written as

[TABLE]

where $1_{k}\in\mathbb{R}^{k}$ is a column vector of ones, $\bar{A}^{i}_{k}\vcentcolon=A^{i}_{k-1,0}$, and

[TABLE]

In terms of the concatenated state vector $X\vcentcolon=\begin{bmatrix}x_{0}^{\mbox{\tiny\sf T}}&\cdots&x_{N}^{\mbox{\tiny\sf T}}\end{bmatrix}^{\mbox{\tiny\sf T}}\in\mathbb{R}^{(N+1)n_{x}}$, control vector $U\vcentcolon=U_{N-1}\in\mathbb{R}^{Nn_{u}}$, and disturbance vector $W\vcentcolon=W_{N-1}\in\mathbb{R}^{Nn_{x}}$, the system dynamics are written as the matrix equation

[TABLE]

The matrices $\mathcal{A}^{i},\mathcal{B}^{i},R^{i}$, and $\mathcal{G}^{i}$ are formed by appropriately stacking the terms from (25).

Let $\mathcal{Q}_{x}$ and $\mathcal{Q}_{u}$ be block-diagonal state and control cost weight matrices with entries corresponding to the continuous weights $Q_{x,t}$ and $Q_{u,t}$ from (5):

[TABLE]

[TABLE]

The quadratic state cost at step $N$ is neglected, since the terminal state is fixed. Letting $E_{k}=\begin{bmatrix}0_{n_{x},kn_{x}},I_{n_{x}},0_{n_{x},(N-k)n_{x}}\end{bmatrix}$ and $E^{u}_{k}=\begin{bmatrix}0_{n_{u},kn_{u}},I_{n_{u}},0_{n_{u},(N-k-1)n_{u}}\end{bmatrix}$, such that $x_{k}=E_{k}X$ and $u_{k}=E^{u}_{k}U$, the continuous time cost functional (8) can be rewritten in terms of the linearized system as

[TABLE]

and the boundary conditions (2) and (4) as

[TABLE]

The chance constraints (3), enforced at each time step $k$, are written as

[TABLE]

In summary, we have approximated the continuous time, nonlinear stochastic system (1) by the discrete, linear stochastic system (29). Problem 1 can be accordingly restated in terms of this approximate system as follows.

Problem 2

Find the control sequence $U^{*}$ that minimizes (32) subject to the dynamics (29), boundary conditions (33), and chance constraints (34).

Remark 1

In the discrete-time formulation, the chance constraints are only enforced at the discrete times $\tau_{k}$, and therefore the original constraints (3) may be violated for some $\tau\in(\tau_{k},\tau_{k+1})$. Constraint violation in this interval is likely when the discretization is too coarse.

III COVARIANCE STEERING

For the remainder of this paper, we will restrict the control law to be of the form [9]

[TABLE]

were $v_{k}\in\mathbb{R}^{n_{u}}$ is a feedforward control, $K_{k}\in\mathbb{R}^{n_{u}\times n_{x}}$ is a feedback gain matrix, and $y_{k}\in\mathbb{R}^{n_{x}}$ is a zero-mean random process given by

[TABLE]

In vector notation, we have

[TABLE]

and thus,

[TABLE]

where the block feedback matrix $K\in\mathbb{R}^{Nn_{u}\times(N+1)n_{x}}$ is given by

[TABLE]

Substituting the control into the state equation, we obtain the expressions for the mean and deviation states as

[TABLE]

Similarly for the control, we obtain

[TABLE]

It follows that the state and control covariances, in terms of the covariance of the process $y_{k}$, are given as

[TABLE]

Substituting (45) and (46) into the cost function (32) and simplifying, we obtain

[TABLE]

where

[TABLE]

III-A Endpoint Constraints

Substituting (40) into (33b), we obtain the equality constraint on the final mean state

[TABLE]

Since the equality constraint $E_{N}\mathcal{P}_{x}E_{N}^{\mbox{\tiny\sf T}}=P_{x_{f}}$ is not convex in $K$, and since in practice a smaller than anticipated state covariance is acceptable, we instead enforce the relaxed inequality constraint [14]

[TABLE]

which is convex in $K$. This constraint may be equivalently stated in the more standard form [8]

[TABLE]

III-B Chance Constraints

Assume that at each time step the convex regions $\mathcal{X}_{k}\vcentcolon=\mathcal{X}_{\tau_{k}}$ and $\mathcal{U}_{k}\vcentcolon=\mathcal{U}_{\tau_{k}}$ can be represented by the finite intersection of half spaces

[TABLE]

where $\mathcal{X}_{k,m}\vcentcolon=\{x\in\mathbb{R}^{n_{x}}:\bar{a}_{k,m}^{\mbox{\tiny\sf T}}x\leq\alpha_{k,m}\}$ and $\mathcal{U}_{k,m}\vcentcolon=\{u\in\mathbb{R}^{n_{u}}:\bar{b}_{k,m}^{\mbox{\tiny\sf T}}u\leq\beta_{k,m}\}$ are given in terms of the vectors $\bar{a}_{k,m}\in\mathbb{R}^{n_{x}}$, $\bar{b}_{k,m}\in\mathbb{R}^{n_{u}}$ and scalars $\alpha_{k,m},\beta_{k,m}\in\mathbb{R}$. By subadditivity of probability, we have

[TABLE]

It follows that if $\mathbb{P}(x_{k}\in\mathcal{X}^{c}_{k,m})\leq p^{x}_{k,m}$ for a set of positive numbers $\{p^{x}_{k,m}\}$ that sum over the index $m$ to less than $p_{x,k}$, then $\mathbb{P}(x_{k}\in\mathcal{X}_{k}^{c})\leq p_{x,k}$ [3, 15]. In terms of the concatenated state and control vectors, and since $x_{k}=E_{k}X$ and $u_{k}=E^{u}_{k}U$, the events $\{x_{k}\in\mathcal{X}_{k}\}\subset\mathbb{R}^{n_{x}}$ and $\{E_{k}X\in\mathcal{X}_{k}\}\subset\mathbb{R}^{(N+1)n_{x}}$ have the same probability. Therefore, when relabeling indices of the inequality constraints according to

[TABLE]

if $\{p^{x}_{m,k}\}$ and $\{p^{u}_{m,k}\}$ are given sets of positive numbers that satisfy, for each $k$,

[TABLE]

then, from (53),

[TABLE]

it follows that

[TABLE]

The same construction applies to the control sets $\mathcal{U}_{k}$.

Next, we formulate the chance constraint $\mathbb{P}(a_{m}^{\mbox{\tiny\sf T}}E_{k}X\leq\alpha_{k,m})$ into a deterministic expression of the control variables. From (29) it follows that $X$ is normally distributed and hence $a_{m}^{\mbox{\tiny\sf T}}E_{k}X$ is a scalar normal random variable with mean $a_{m}^{\mbox{\tiny\sf T}}E_{k}\bar{X}$ and covariance $a_{m}^{\mbox{\tiny\sf T}}E_{k}\mathcal{P}_{x}E_{k}^{\mbox{\tiny\sf T}}a_{m}$. It follows that

[TABLE]

where $\operatorname{cdfn}$ is the cumulative normal distribution function. Therefore, the chance constraint $\mathbb{P}(a_{m}^{\mbox{\tiny\sf T}}E_{k}X\leq\alpha_{k,m})\geq 1-p_{x,k}$ can be equivalently written as

[TABLE]

Putting it all together, if (56) holds, and if

[TABLE]

then we ensure that the chance constraints (34) will be satisfied.

Remark 2

This work assumes that $\{p^{x}_{m,k}\}$ and $\{p^{u}_{m,k}\}$ are given sets of positive numbers that satisfy (56). This assumption allows (61) to be convex. Otherwise, (61) becomes non-convex, and we need to consider an optimal risk allocation problem. Several approaches have been proposed, such as [16, 17], to handle the risk allocation problem. In addition, the authors of [18] used a primal-dual interior point method to find an optimal risk allocation.

We are now ready to restate the covariance steering problem as a deterministic, finite dimensional optimization problem.

Problem 3

Linear Covariance Steering. Find $K^{*}$ and $V^{*}$ that minimize the cost (47) subject to the terminal state constraints (49) and (51) and the chance constraints (61).

IV ITERATIVE COVARIANCE STEERING

In the previous sections, we have locally approximated the continuous time, nonlinear system (1) with the discrete linear system (29), and we have restated the cost function and constraints in terms of the discrete linear system as functions of a feedfoward control $V$ and feedback gain $K$. We will search for solutions to the original nonlinear system by successively solving this approximate convex problem, where the optimal controls from each successive problem are used to propagate trajectories of the nonlinear system to obtain references for the next linearization step. This method is referred to in the literature as successive convexification [19, 20].

IV-A Stochastic Trust Region

The linear approximation of the system dynamics is only valid in a neighborhood around the reference trajectory, so care must be taken to ensure that the optimal controls for the linear problem are relevant to the nonlinear problem. For this reason, variations in the state and the control from the previous solution are bounded inside a trust region [19]. In this paper, we are successively approximating a stochastic system, and since the state of a system with Brownian noise is unbounded, we must define a stochastic trust region instead as follows

[TABLE]

where $p^{x}_{\mathrm{tr}}$, $p^{u}_{\mathrm{tr}}$, $\Delta_{x}^{i}$, and $\Delta_{u}^{i}$ are user-defined limits. Constraints in the 1-norm can be represented by $2n_{x}$ or $2n_{u}$ inequality constraints for the state or control, respectively, using (61). As a consequence of these trust region constraints, if the reference trajectory is sufficiently far away from the terminal constraint, then the problem may become infeasible. For these situations, which are most likely encountered when initializing the problem, we relax the hard constraint (49) on the terminal state mean to the soft constraint

[TABLE]

with a corresponding term $\eta_{x_{f}}w_{x_{f}}$ added to the cost, where $w_{x_{f}}$ is a user-defined weight. This constraint may be replaced with the hard constraint (49) when the reference trajectory $\hat{x}^{i}$ is sufficiently close to the terminal constraint. In the case (63) is active, we use the augmented cost function given by

[TABLE]

In summary, we have modified Problem 3 to the following convex optimization problem.

Problem 4

iCS Convex Subproblem. Find $K^{*}$ and $V^{*}$ that minimize the cost (64) subject to the terminal state constraints (61) and (49) (or (63) if the reference trajectory is sufficiently far from the target), the chance constraints (61), and the trust region constraints (62).

This problem is solved successively in order to find a solution to Problem 1 using the iCS algorithm presented in Procedure 1.

Remark 3

The nonlinear mean dynamics can be propagated through Monte Carlo, which can be parallelized. In the case when computational resources are limited, we can approximate $\mathbb{E}[f(x_{t},u_{t},t)]\approx f(\bar{x}_{t},\bar{u}_{t},t)$ so that the mean state evolves according to

[TABLE]

In this case, the mean state can be estimated by integrating a single trajectory.

V NUMERICAL EXAMPLE

In this section we apply the iCS algorithm to control a double integrator subject to a quadratic drag force. Let the position $\xi\in\mathbb{R}^{2}$ and velocity $v\in\mathbb{R}^{2}$ be described by the stochastic system

[TABLE]

where $c_{d}>0$ is the drag coefficient and $\gamma>0$ is a noise scale parameter. In terms of the state $x=(\xi,v)\in\mathbb{R}^{4}$, the dynamics can be written as

[TABLE]

where $f_{d}$ represents the nonlinear drag dynamics and where

[TABLE]

Linearizing about a reference velocity $\hat{v}^{i}_{\tau}$, we obtain

[TABLE]

where $E_{2}^{\mbox{\tiny\sf T}}=\begin{bmatrix}0_{2},I_{2}\end{bmatrix}$, $B^{i}_{\tau}=\sigma B$, and $r^{i}_{\tau}$ is given as in (11). In addition, we enforce the the chance constraint

[TABLE]

where $e_{1}=\begin{bmatrix}1,0\end{bmatrix}$. The initial state is normally distributed with mean and covariance

[TABLE]

and the terminal distribution is constrained by the mean and covariance

[TABLE]

We set the drag coefficient $c_{d}=0.005$ and the noise scale $\gamma=0.01$. For the solution, we let the number of discrete steps $N=25$, terminal mean error weight $w_{x_{f}}=1000$, and time scale $\sigma=15$. The mean cost function was

[TABLE]

the weight matrices were $Q_{x,\tau}\equiv 5I$ and $Q_{u,\tau}\equiv I$, and the algorithm was seeded with the initial guess

[TABLE]

Since the initial guess violates the chance constraint (71), we relax the chance constraint for the first iteration and tighten it to the final constraint over the first several iterations. The algorithm converged in five iterations, and solutions for each iteration are shown in Figure 1. Samples from a 5,000 trial Monte Carlo simulation are shown in Figure 2. The maximum probability of constraint violation was at step $k=11$, with 9.14% of states having $\left\lVert e_{1}\xi_{k}\right\rVert_{1}\geq 6$, which is below the limit of 10% set in (71). Also from the Monte Carlo simulation, the final state mean was

[TABLE]

which is very close to the specified value in (73), and the covariance

[TABLE]

is less than the upper bound specified in (73).

VI CONCLUSION

In this paper we presented an algorithmic solution to the chance constrained nonlinear CS problem. We began by approximating the original nonlinear stochastic system by a linear discrete stochastic system, and then we formulated the linear CS problem as a deterministic optimization problem. Next, in order for the linearized problem formulation to be a reasonable approximation of the original nonlinear problem, we constrained the trajectory at each iteration within a probabilistic trust region about the trajectory from the previous iteration. The size of the trust region depends on the nonlinearity of the dynamics, and therefore convergence properties are problem-specific.

Since the proposed iCS algorithm linearizes the dynamics at each iteration, an initial trajectory must be given for the first iteration of the algorithm. At the same time, the difference in the trajectories between iterations is constrained within a trust region, and so a poor initialization may cause the first step to be infeasible. In the numerical example we addressed this problem by relaxing the chance constraints in the first iterations of the algorithm. Since the chance constrained region is assumed to be a convex polytope, the region can be easily expanded by scaling the inequality constraints. Another solution would be to first solve a deterministic optimization problem with tightened inequality constraints representing a worst-case chance constrained region. In this case, the iCS algorithm would be used to improve the solution from the deterministic problem by adding the closed-loop system statistics to the optimization. The latter approach could be applied to problems such as planetary entry and powered descent by iterating on a given reference trajectory that is to be tracked in the presence of uncertainty.

In future work we plan to apply iCS to problems in entry, descent, and landing (EDL) with nonlinear dynamics, such as entry and powered descent [21]. Another extension to this work would be the addition of time-varying chance constraints that are satisfied for all time, rather than for each time, while not being overly conservative.