A Simulation Study on Turnpikes in Stochastic LQ Optimal Control

We consider stochastic LQ OCPs of the following form

whereby at each time step $k$, $X(k)\in L^{2}(\Omega,\mathcal{F}_{k},\mu;\mathbb{R}^{n_{x}})$ is an $i.i.d.$ random variable on the underlying filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{k})_{k=0,\ldots,N},\mu)$, where $\Omega$ is the set of realizations, $\mathcal{F}$ is a $\sigma$-algebra, $(\mathcal{F}_{k})_{k=0,\ldots,N}$ is a filtration, and $\mu$ is the probability measure.

Note that the $\sigma$-algebras $\mathcal{F}_{k}$ are related by the time evolution of the information, i.e.

We choose $(\mathcal{F}_{k})_{k=0,\ldots,N}$ as the smallest filtration such that $X$ is an adapted process (which results form the evolution of the dynamics (1b)), i.e.

Then, the control at time $k$ is modeled as a stochastic process which is adapted to the filtration $(\mathcal{F}_{k})_{k}$, i.e. $U(k)\in L^{2}(\Omega,\mathcal{F}_{k},\mu;\mathbb{R}^{n_{u}})$.

The concept of stochastic filtrations can be understood as a stochastic causality requirement, i.e. the the stochastic input process $U(k)$ at time step $k$ may only depend on the realization of the random variables $X(0),\dots,X(k)$ up to time step $k$. Note that the influence of the noise, which is also an adapted stochastic process, is handled implicitly via the state recursion. That is, $U(k)$ depends on $X(0),\dots,X(k)$ and thus also on $W(0),\dots,W(k-1)$. For details on stochastic filtrations we refer to (Fristedt and Gray, 2013).

The stage cost is given by

Indeed typical choice for $\ell$ are a combination of expected value and variance of some underlying deterministic stage cost function. The considered dynamics are a time-discrete stochastic system subject to noise $W(\cdot)$ modelled as a stochastic process. The constraints are written as chance constraints for states and inputs, whereby the underlying sets $\mathbb{X}$ and $\mathbb{U}$ are assumed to be closed. Moreover, $\varepsilon_{x}$ and $\varepsilon_{u}$ specify the probabilities with which the chance constraints shall be satisfied.

In order to obtain a tractable reformulation of the stochastic LQ OCP (1) we consider the framework of Polynomial Chaos Expansions (PCE). For an in-depth introduction we refer to (Sullivan, 2015). The underlying idea of PCE is that random variables from some $L^{2}(\Omega,\mathcal{F},\mu;\mathbb{R}^{m})$ with $m\in\{n_{x},n_{u},n_{w}\}$ can be described using an appropriate basis. To this end, we consider an orthogonal polynomial basis $\{\phi_{i}\}_{i=0}^{\infty}$ which spans $L^{2}(\Omega,\mathcal{F},\mu;\mathbb{R}^{m})$.

In order to obtain a computationally tractable formulation, one has to truncate the PCE series after $L+1$ terms

The choice of the basis polynomials can be inferred via the Wiener-Askey scheme (Sullivan, 2015). For example, in case of a standard Gaussian random variable one would consider Hermite polynomials as they allow modeling the Gaussian with the first two terms of the PCE series. Moreover, in case of explicit polynomial maps—e.g., consider state transitions $X(k)=A^{k}X(0)$—of finite degree, one can quantify the truncation errors arising from considering only the first $L+1$ PCE coefficients. One may also infer $L\in\mathbb{N}$ such that no truncation error arises. We refer to Mühlpfordt et al. (2018) for details.

Finally, we remark that whenever the truncated series representation (3) is exact, the first two moments of a random variable $X$ can be computed in terms of PCE coefficients as follows

Consider the stochastic variant of the motivating example

where $X_{0}$ denotes the initial random variable with known probability distribution $p_{X_{0}}$, $W(k)$ denotes system noise modeled as a white Gaussian noise such that all $W(k)$ have an identical known probability distribution $p_{W}$. We arrive at the following stochastic LQ OCP

where ${x}_{min}=-2$ and ${x}_{max}=2$. We approximate the chance constraint as

with $\lambda=\sqrt{(1+\varepsilon_{x})/(1-\varepsilon_{x})}$, for a derivation see (Farina et al., 2013).

We consider $X_{0}$ to follow a uniform distribution with the support $[{x}_{0min},{x}_{0max}]=[0.6,1.4]$. The noise at time $k\in\mathbb{I}_{[0,N-1]}$, $W(k)$, is a Gaussian distribution with mean $\mathbb{E}[W(k)]=0$ and variance $\mathbb{V}[W(k)]=0.5^{2}$. Additionally, $\varepsilon_{x}$ is set to 0.8 and thus we have $\lambda=3$.

Without the noise, considering a first-order PCE ($L=1$) for $X(k)$ and $U(k)$ with identical basis functions, exactness of the PCE representation of $X(k+1)$ is guaranteed since the system dynamic is linear, see Mühlpfordt et al. (2018). Therefore, including noise, the PCE dimension needed for an exact representation is determined by the horizon $N$. More precisely, $L=N+1$, where two PCE terms are induced by the uncertainty of the initial condition and the rest is caused by the noise. The PCE basis and coefficients read

where $\xi$ is a standard uniformly distributed random variable while $\theta_{i}$ are independent Gaussian random variables. Additionally, $\textsf{x}_{00}=({x}_{0min}+{x}_{0max})/2$ and $\textsf{x}_{01}={x}_{0max}-{x}_{0min}$, since $X_{0}$ is uniformly distributed.

Stochastic modeling via adapted filtrations expresses the idea that the noise $W(k)$ at time $k$ influences the next state $X(k+1)$ and input $U(k+1)$ and the subsequent time instances $i>k$ but not at $i\leq k$. In terms of PCE representation this implies that

with

and where, respectively, $Z\in\{X,U\}$ and $\textsf{z}\in\{\textsf{x},\textsf{u}\}$ are placeholders and $M=N$ for $X$ and $M=N-1$ for $U$.

To illustrate turnpike behavior, we solve the stochastic OCP (5) as given above over different optimization horizons $N=3,6,...,24$ using PolyChaos.jl (Mühlpfordt et al., 2020). Fig. 2 shows the trajectories of the optimal solutions for a total of 16 realizations of the uncertainties. At first glance, the realizations of the solutions are noisy and appear not to exhibit the turnpike property.

Hence, in order to uncover the turnpike, we consider the trajectories of the PCE coefficients. Note that each PCE basis $\phi_{i}^{w}$ induced by the noise is a standard Gaussian distributed random variable, the sum $\sum_{i=0}^{k-1}\textsf{z}_{i}^{w}(k)\phi_{i}^{w},\textsf{z}\in\{\textsf{x},\textsf{u}\}$ is equal to a new Gaussian distributed random variable with zero mean and variance. Hence, for the sake of simplified illustration, we consider

Here $\phi_{\Sigma}$ is a random variable with standard Gaussian distribution (zero mean). Therefore, instead of $\textsf{z}_{i}^{w}(k),i\in\mathbb{I}_{[0,N-1]}$, only one PCE coefficient $\textsf{z}_{\Sigma}^{w}(k)$ suffices to represent the uncertainty caused by noise. Note that this transformation is used only for illustration priposes and not in the underlying computation.

Using the PCE reformulation detailed above, Fig. 3 illustrates that the turnpike phenomenon occurs in terms of PCE coefficients. Actually, the turnpike property of PCE coefficients suggests that the optimal steady-state is a random variable with stationary distribution. This distribution can be calculated from the optimal steady-state problem formulated via PCE. Doing so, i.e. solving

we obtain the stationary distribution depicted in Fig. 4. This figure also shows the histogram at $k=25$ obtained from sampling $10^{4}$ realizations of the uncertainty from the optimal PCE solution for $N=50$. As one can see, the behavior in the middle of the horizon corresponds to the solution obtained for the steady state problem.

The time evolution of the state histograms and the distributions obainted via PCE is shown in Fig. 5 for $k=0,10,20,30,40,50$. It is not surprising that the histograms follow the calculated PDF perfectly and the state $X(k)$ keeps the same distribution in the middle of trajectory. This illustrates the turnpike phenomenon in the distributions.

As a second example, we consider a linearized CSTR reactor which appeared in several papers such as (Zanon and Faulwasser, 2018). The expected value quadratic stage cost and the linear discrete-time system with noise read

where $W_{A},W_{B}$ are modeled as independent uniformly distributed random variables with support $[-0.1,0.1]$. The initial state $X(0)$ is a Gaussian distributed random vector with known mean $(0.5,0.8)$ and variance $(0.05^{2},0.08^{2})$. We solve the stochastic OCP via PCE over different horizons $N=10,20,...,80$. We obtain the trajectories of the optimal solutions for a total $16$ realization sequences of the uncertainties, see Fig. 6. As the dimension of PCE is quite large, we plot the first two moments instead of PCE coefficients, i.e. mean and the variance of state and input random variables, see Fig. 7. Similar to the previous example, the trajectories of mean and variance exhibit the turnpike property. Similar to before, Fig. 8 depicts the histograms of the state $X_{A}$ at $k=0,10,20,30,40,50$ for the considered realizations and the PDF obtained via PCE, see Fig. 8.

As an additional means of assessing the turnpike phenomenon via simulation, we consider the following numerical setting:

• •

  Compute a random realization of the disturbance $W(k)$, for $k\in\mathbb{I}_{[0,M]}$ denoted as $\{W(k;\omega_{k})\}$

• •

  Pick horizons $N_{1},...,N_{s}\leq M$ and corresponding realizations of the initial condition $X_{0}(\omega_{1}),\dots,X_{0}(\omega_{s})$.

• •

  For all horizons $N_{1},\dots,N_{s}$ and the initial condition $X_{0}(\omega_{1}),\dots,X_{0}(\omega_{s})$, simulate the response of the dynamics under the optimal input policy, while the disturbance sequence is fixed to $\{W(k;\omega_{k})\}$ (or a subpart thereof).

The results of this numerical experiment are depicted in Figure 9.

As one can see, all the trajectories approach the same solution after some time. One can understand this solution as a time-varying path of a stationary turnpike solution, whose shape is governed by the considered disturbance sequence. Observe the difference to Figure 6, wherein for each trajectory a different disturbance realization sequence is considered.

What could we do if we want to get a optimal steady-state with a narrow distribution, or in other words, with small variance? Involving variance penalization in the objective function is one option. We consider the previous example augmented with a penalty of the variance of the state $X_{A}$ in the stage cost

Here we choose $\gamma=10^{4}$ and solve the stochastic OCP via PCE over different horizons $N=10,20,...,80$. The trajectories of the optimal solutions for a total $16$ realization sequences of the uncertainties are shown in Fig. 10. Observe that the state $x_{A}$ is much less effected by the noise. It can also be seen in Fig. 11 that the variance of state $X_{A}$ is smaller than in Example 2, while the variance of state $X_{B}$ and the variance of the input $U$ increase. Fig. 12 also shows that the optimal steady state $X_{A}$ has a narrower distribution.