Randomized Algorithms for Data-Driven Stabilization of Stochastic Linear Systems

TL;DR

This paper introduces two randomized algorithms for stabilizing stochastic linear systems with unknown parameters, analyzing their performance and showing that sufficient randomizations ensure fast stabilization.

Contribution

It provides a detailed analysis of the stabilization speed and failure probability of two novel data-driven randomized control algorithms.

Findings

Fast stabilization is guaranteed with enough independent randomizations.

Performance depends on the magnitude and frequency of randomizations.

Numerical analysis supports the effectiveness of the proposed methods.

Abstract

Data-driven control strategies for dynamical systems with unknown parameters are popular in theory and applications. An essential problem is to prevent stochastic linear systems becoming destabilized, due to the uncertainty of the decision-maker about the dynamical parameter. Two randomized algorithms are proposed for this problem, but the performance is not sufficiently investigated. Further, the effect of key parameters of the algorithms such as the magnitude and the frequency of applying the randomizations is not currently available. This work studies the stabilization speed and the failure probability of data-driven procedures. We provide numerical analyses for the performance of two methods: stochastic feedback, and stochastic parameter. The presented results imply that as long as the number of statistically independent randomizations is not too small, fast stabilization is guaranteed.

Equations

$$x(t+1)$$

$$K\left(\mathbf{\theta }\right)=Q+A^{\prime }K\left(\mathbf{\theta }\right)A-A^{\prime }K\left(\mathbf{\theta }\right)B{\left(B^{\prime }K\left(\mathbf{\theta }\right)B+R\right)}^{-1}{B}^{\prime }K\left(\mathbf{\theta }\right)A.$$

$$L\left(\mathbf{\theta }\right)$$

$${\hat{D}}_i=\arg \underset{D\in {\mathbb{R}}^{p\times q}}{\min }\sum _{\ell =(i-1)\lfloor T/k\rfloor }^{i\lfloor T/k\rfloor -1}{\left|\left|x(\ell +1)-Dx(\ell )\right|\right|}^2.$$

$${\hat{\mathbf{\theta }}}_T=\arg \underset{\mathbf{\theta }\in {\mathbb{R}}^{p\times q}}{\min }\sum _{i=1}^k{\left|\left|\left|{\hat{D}}_i-\mathbf{\theta }\left[\begin{array}{c}{I}_p\\ L_i\end{array}\right]\right|\right|\right|}^2.$$

$$\mathbb{P}\left(\left|\left|\left|{\hat{\mathbf{\theta }}}_T-{\mathbf{\theta }}_0\right|\right|\right|>{ϵ}_0\right)\le e^{-T{c}_0}.$$

$${\hat{\mathbf{\theta }}}_T=\arg \underset{\mathbf{\theta }\in {\mathbb{R}}^{p\times q}}{\min }\sum _{i=1}^k{\left|\left|\left|{\hat{D}}_i-\mathbf{\theta }\left[\begin{array}{c}{I}_p\\ L\left({\mathbf{\theta }}_i\right)\end{array}\right]\right|\right|\right|}^2.$$

$$A_0=\left[\begin{array}{ccc}1.07& 0& -0.37\\ 0.48& -0.88& 0.85\\ 0& 0.03& -0.92\end{array}\right],$$

$$Q=\left[\begin{array}{ccc}1.31& -0.17& -0.28\\ -0.17& 1.14& 0.51\\ -0.28& 0.51& 5.01\end{array}\right],$$

$$K\left({\mathbf{\theta }}_0\right)=\left[\begin{array}{ccc}2.83& 0& -0.87\\ 0& 2.20& -0.32\\ -0.87& -0.32& 7.31\end{array}\right],$$

Full text

I Introduction

Sequential decision-making strategies for stochastic linear systems is an ubiquitous model extensively used in different fields such as biology, robotics, finance, and cryptography [1, 2, 3, 4, 5]. The setting consists of a state-space dynamical model following the stochastic trajectory

[TABLE]

That is, the current state vector $x(t)\in\mathbb{R}^{p}$ together with the current action $u(t)\in\mathbb{R}^{r}$ lead the system to the next state $x(t+1)$, according to (1). The linear dynamical model is being disturbed by the stochastic noise process $\left\{\xi(t)\right\}_{t=1}^{\infty}$. The state transition matrix $A_{0}$, and the action influence matrix $B_{0}$, both of appropriate dimensions, are unknown.

The objective is to analyze the performance of adaptive stabilization procedures. The focus is on randomized algorithms for stabilizing the system by designing the control actions according to the observed state sequence. Note that the control action can cause instability of the state vector. That is due to the uncertainty of the decision maker about the dynamics matrices $A_{0},B_{0}$. After the stabilization procedure, the control policy for the operating plant under consideration needs to be designed [6, 7, 8]. The real applications determine the goal being for example minimizing the input energy or regulating the state vector. Therefore, a reliable stabilization algorithm is required to be executed prior to the regulation step, and is desired to guarantee the stabilization quite fast. In general, design of effective methods for finite time stabilization has been recognized as a difficult problem [9, 10].

Next, we briefly review various algorithms proposed in the literature for stabilization. Early work consider a restricted setting where the operator norm of the closed-loop matrix is less than one [11]. Then, employing parameter estimation and policy design in an alternate manner, the system is shown to be stabilized [12]. Ensuing work establishes the theoretical guarantee for a randomized algorithm based on stochastic feedbacks (SF) [13]. In case of availability of many independent trajectories, stabilizing feedbacks can be learned using only the last sample of each trajectory [14]. Further, reinforcement learning algorithms relying on other randomization methods such as stochastic parameters (SP) [15], and input perturbation [16, 17] are presented in the literature. The former designs the stabilizing feedback for a randomly generated dynamics parameter, while the latter adds an independent dither signal to the control command. Both procedures then learn the true dynamical parameters through observing the state sequence.

This work studies the performance of SF and SP which compared to the other aforementioned methods, are able to stabilize a larger class of linear systems. In both algorithms, $k$ random control policies are being applied to the system for episodes of equal lengths.111$k\geq 1+\lceil r/p\rceil$ is required, see Section III for details Although some theoretical results exist in the literature [13, 15], practicality of these stabilization algorithms is not studied. Further, analysis of SF and SP is of independent interest because the underlying stochastic systems are unstable and time-varying. This work aims to address the above gap through considering the implementation of the algorithms. We provide numerical analysis to compare the amount of time needed for the stabilization. Further, we study the learning accuracy, the failure probability, and the effect of both the number of episodes, and the standard deviation (of the stochastic feedbacks and the stochastic parameters).

This paper is organized as follows. In Section II we formulate the stabilization problem this work is addressing and discuss the theoretical backgrounds. In Section III, two randomized algorithms to stabilize the system in finite time are discussed. Then, in Section IV the numerical analysis of the presented algorithms are provided by depicting stabilization speeds, learning accuracies, and failure probabilities.

Notation: For square matrix $M$, the largest eigenvalue of $M$ (in magnitude) is denoted by $\lambda_{\max}(M)$. For vector $v$, ${\left|\kern-1.29167pt\left|v\right|\kern-1.29167pt\right|}$ is the Euclidean norm, and the operator norm of matrices is denoted by ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|M\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}=\sup\limits_{{\left|\kern-0.90417pt\left|v\right|\kern-0.90417pt\right|}=1}{\left|\kern-1.29167pt\left|Mv\right|\kern-1.29167pt\right|}$. For the notation $\bm{\theta}$, see Definition 1. Finally, $\mathcal{N}\left(\mu,\Sigma\right)$ denotes the normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$.

II Theoretical Framework

To set the stage, we briefly discuss the theoretical framework of the problem. First, the system needs to be stabilizable:

Assumption 1**.**

There exists a $r\times p$ stabilizer matrix $L$. That is, $\left|\bm{\lambda}_{\max}\left(A_{0}+B_{0}L\right)\right|<1$.

The noise in (1) is a sequence of centered independent random vectors; $\mathbb{E}\left[\xi(t)\right]=0$, satisfying the following.

Assumption 2**.**

The noise covariance matrices are positive definite; $\left|\bm{\lambda}_{\min}\left(\mathbb{E}\left[\xi(t)\xi(t)^{\prime}\right]\right)\right|>0$. Further, the noise is sub-Gaussian: $\sup\limits_{t\geq 1}\mathbb{E}\left[\mathrm{exp}\left({\beta{\left|\kern-1.29167pt\left|\xi(t)\right|\kern-1.29167pt\right|}^{2}}\right)\right]<\infty$, for some $\beta>0$.

Given the time length $T$, the objective is to design the control sequence $\left\{u(t)\right\}_{t=0}^{T-1}$ effectively. That is, one can stabilize the system utilizing the observations $\left\{x(t)\right\}_{t=0}^{T}$ acquired by applying $\left\{u(t)\right\}_{t=0}^{T-1}$. For this purpose, we employ a model-based approach so that the algorithm learns $A_{0},B_{0}$ through observing the history by date: $\mathcal{H}_{t}=\left\{\left\{x(i)\right\}_{i=0}^{t},\left\{u(i)\right\}_{i=0}^{t-1}\right\}$. The subsequent argument shows that a coarse approximation suffices for stabilization. To proceed, we define a notation:

Definition 1**.**

Henceforth, for all $A\in\mathbb{R}^{p\times p}$, $B\in\mathbb{R}^{p\times r}$, we use $\left[A,B\right]$ and $\bm{\theta}\in\mathbb{R}^{p\times q}$ interchangeably, where $q\equiv p+r$.

It is well known that the system can be stabilized by solving a Riccati equation [18]. First, let $Q\in\mathbb{R}^{p\times p},R\in\mathbb{R}^{r\times r}$ be arbitrary positive definite matrices. Then, define $K\left(\bm{\theta}\right)$ as a possible solution of the algebraic Riccati equation

[TABLE]

Accordingly, for an arbitrary dynamics parameter $\bm{\theta}$, define the linear feedback matrix

[TABLE]

The matrix $K\left(\bm{\theta}\right)$ exists for almost all $\bm{\theta}$ (w.r.t. Lebesgue measure [19]), and $L\left(\bm{\theta}\right)$ stabilizes the system:

Lemma 1**.**

[13]** There exists a fixed $\epsilon_{0}>0$, such that an $\epsilon_{0}$-accurate approximation ensures the stabilization. That is, ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widehat{\bm{\theta}}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}\leq\epsilon_{0}$ implies that $\left|\bm{\lambda}_{\max}\left(A_{0}+B_{0}L\left(\widehat{\bm{\theta}}\right)\right)\right|<1$.

III Algorithms

The system can be stabilized by leveraging Lemma 1, together with the recent results on fast identification of unstable systems [20, 17]. In the sequel, we explain two different randomization methods, and compare their performance through numerical simulations. Note that since during the stabilization period the system is presumably unstable, the state vector can grow unbounded.

In order to learn the unstable linear dynamical systems, a sufficient and necessary condition is the closed-loop regularity [21]. In regular matrices, all eigenspaces corresponding to the eigenvalues of magnitude larger than one are one dimensional. To satisfy the regularity condition, it suffices to apply stochastic feedbacks [13]. The stabilizing algorithm proposed in the above reference utilizes a few random feedback matrices to accurately learn $\bm{\theta}_{0}$.

Namely, suppose that we want to stabilize the system in $T$ time steps. The algorithm applies $k$ statistically independent Gaussian feedback matrices $L_{1},\cdots,L_{k}$ in episodes of the equal length $\lfloor T/k\rfloor$. Collecting the history $\mathcal{H}_{k\lfloor T/k\rfloor}$, the algorithm first estimates the closed-loop matrix $A_{0}+B_{0}L_{i}$:

[TABLE]

Then, the closed-loop estimates $\widehat{D}_{1},\cdots,\widehat{D}_{k}$ are being used to estimate the true dynamics parameter $\bm{\theta}_{0}$ according to:

[TABLE]

The corresponding pseudo-code is provided in Algorithm 1. To ensure that the learning procedure leads to an accurate approximation, the number of episodes $k$ satisfies $k\geq 1+\lceil r/p\rceil$ [13]. The following high probability result is established for the performance of the stochastic feedback (SF) method.

Lemma 2**.**

[13]** Defining $\epsilon_{0}$ similar to Lemma 1, let $\widehat{\bm{\theta}}_{T}$ be the output of Algorithm 1. Then, there is $c_{0}>0$ such that

[TABLE]

Putting Lemma 1 and Lemma 2 together, we conclude that the failure probability of SF decays exponentially with $T$.

Here, we also study another randomization procedure based on stochastic parameters (SP). In this method, the algorithm draws Gaussian independent parameters $\bm{\theta}_{1},\cdots,\bm{\theta}_{k}$. Then, in each episode $i$, SP applies the control $u(t)=L\left(\bm{\theta}_{i}\right)x(t)$, which is computed according to (2). At the end of every episode, the closed-loop matrix is estimated by (3). Finally, putting the estimates $\widehat{D}_{1},\cdots,\widehat{D}_{k}$ together, SP learns $\widehat{\bm{\theta}}_{T}$ by the minimization problem

[TABLE]

Similarly, the lower bound $k\geq 1+\lceil r/p\rceil$ is necessary to ensure identifiability of $\widehat{D}_{T}$. Algorithm 2 depicts the implementation of SP.

IV Numerical Analysis

In this section we simulate SF and SP for the dynamics matrices

[TABLE]

and the cost matrices

[TABLE]

Entries of $A_{0},B_{0},Q,R$ are rounded, after being randomly generated by MATLAB. For this system, solving the Riccati equation, we get

[TABLE]

Then, if one designs the feedback gain according to the true parameter $\bm{\theta}_{0}$, the spectral radius of the closed-loop matrix is $\left|\bm{\lambda}_{\max}\left(A_{0}+B_{0}L\left(\bm{\theta}_{0}\right)\right)\right|=0.51$.

Figure 1 illustrates the statement of Lemma 1. The largest eigenvalue of the closed-loop matrix is reported as a function of the error in the dynamics parameter. The scatter plot consists of different error values ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widehat{\bm{\theta}}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$, and the resulting spectral radius $\left|\bm{\lambda}_{\max}\left(A_{0}+B_{0}L\left(\widehat{\bm{\theta}}\right)\right)\right|$.

Figure 2, and Fig 3 plot the performance of Algorithm 1 for $100$ replications, as the number of time steps $T$ grows. The accuracy of learning the unknown dynamics parameter $\bm{\theta}_{0}$ is reported in Figure 2 for different values of $k,\sigma$. Further, Figure 3 graphs the number of replications in which the system is stabilized. This can be interpreted as a surrogate for the stabilization probability in Lemma 2.

Similarly, Figure 4 and Figure 5 illustrate the performance of Algorithm 2. The results on the learning error ${\left|\kern-1.29167pt\left|\kern-1.29167pt\left|\widehat{\bm{\theta}}_{T}-\bm{\theta}_{0}\right|\kern-1.29167pt\right|\kern-1.29167pt\right|}$ are depicted in Figure 4. Then, Figure 5 plots the percentage of the $100$ replications that the linear feedback designed by the output of Algorithm 2 stabilizes the true dynamics $\bm{\theta}_{0}$. Note that in both Figure 2 and Figure 4, the graph for $k=2$ is omitted since in this case the learning error is much larger than the other cases $k=3,4,5$.

According to the above figures, as long as $k\geq 2$, stabilization will eventually occur. However, the case $k=2$ is significantly slower, and less accurate, compared to the others. All three cases $k=3,4,5$ although, exhibit similar performances. On the other hand, comparing different values of $\sigma$, it is clear that the magnitude of the random matrices ($L_{i}$ in SF and $\bm{\theta}_{i}$ in SP) does not have any remarkable influence on the performance. More precisely, different values of $\sigma$ are indifferent as long as they are not too small, or too large. If $\sigma$ is very small, the randomization will be masked by the randomness of the noise process $\left\{w(t)\right\}_{t=1}^{\infty}$, and is not practical. Further, for large values of $\sigma$, the spectral radius of the resulting closed-loop is wildly large, and leads to an immediate explosion.