Distribution Modeling and Stabilization Control for Discrete-Time Linear Random Dynamical Systems Using Ensemble Kalman Filter

TL;DR

This paper develops a stabilization control framework for discrete-time linear stochastic systems using an ensemble Kalman filter for distribution modeling and feedback gain design, demonstrated through numerical experiments.

Contribution

It introduces a novel output feedback stabilization method combining EnKF-based distribution modeling with feedback control for stochastic systems.

Findings

Effective stabilization demonstrated through numerical experiments

EnKF successfully models system distribution for control synthesis

Framework advances learning-based control for stochastic systems

Abstract

This paper studies an output feedback stabilization control framework for discrete-time linear systems with stochastic dynamics determined by an independent and identically distributed (i.i.d.) process. The controller is constructed with an ensemble Kalman filter (EnKF) and a feedback gain designed with our earlier result about state feedback control. The EnKF is also used for modeling the distribution behind the system, which is required in the feedback gain synthesis. The effectiveness of our control framework is demonstrated with numerical experiments. This study will become the first step toward the realization of learning type control using our stochastic systems control theory.

Equations

$$x_{k+1}=A({\xi }_k)x_k+B({\xi }_k)u_k$$

$$u_k=F{x}_k$$

$$x_{k+1}=A_{\mathrm{cl}}({\xi }_k)x_k,A_{\mathrm{cl}}({\xi }_k)=A({\xi }_k)+B({\xi }_k)F.$$

$$E[A_{ij}({\xi }_k{)}^2]<\infty ,E[B_{ij}({\xi }_k{)}^2]<\infty ,$$

$$\sqrt{E[||x_k|{|}^2]}\le a||x_0||{\lambda }^k(\forall k\in {\mathbf{N}}_0,\forall x_0\in {\mathbf{R}}^n).$$

$$E[{\lambda }^{2}P-A_{\mathrm{cl}}({\xi }_0{)}^{T}P{A}_{\mathrm{cl}}({\xi }_0)]\ge 0.$$

$$\left[\begin{array}{cc}{\lambda }^{2}X& \ast \\ {\bar{G}}_A^{\prime }X+{\bar{G}}_B^{\prime }Y& X\otimes I_{\bar{n}}\end{array}\right]\ge 0$$

$${\bar{G}}_A^{\prime }:=[{\bar{G}}_{A1}^T,\dots ,{\bar{G}}_{An}^T{]}^T\in {\mathbf{R}}^{n\bar{n}\times n},$$

$${\bar{G}}_B^{\prime }:=[{\bar{G}}_{B1}^T,\dots ,{\bar{G}}_{Bn}^T{]}^T\in {\mathbf{R}}^{n\bar{n}\times m},$$

$$\bar{G}=:\left[{\bar{G}}_{A1},\dots ,{\bar{G}}_{An},{\bar{G}}_{B1},\dots ,{\bar{G}}_{Bn}\right]$$

$$({\bar{G}}_{Ai}\in {\mathbf{R}}^{\bar{n}\times n},{\bar{G}}_{Bi}\in {\mathbf{R}}^{\bar{n}\times m}(i=1,\dots ,n))$$

$$\displaystyle\bar{G}^{T}\bar{G}=E\big{[}[{\rm row}(A(\xi_{0})),{\rm row}(B(\xi_{0}))]^{T}$$

$$\displaystyle\cdot[{\rm row}(A(\xi_{0})),{\rm row}(B(\xi_{0}))]\big{]}.$$

$${\psi }_k=f_k({\psi }_{k-1})+w_k$$

$$z_k=h_k({\psi }_k)+v_k$$

$$w_k\sim \mathcal{N}(0,Q),v_k\sim \mathcal{N}(0,R)$$

$$p({\psi }_k|{\mathbf{y}}_{1:k-1})\simeq \frac{1}{M}\sum _{i=1}^M\delta ({\psi }_k-{\psi }_{k|k-1}^{(i)}),$$

$${\psi }_{k|k-1}=\frac{1}{M}\sum _{i=1}^M{\psi }_{k|k-1}^{(i)},$$

$${\psi }_{k|k-1}^{(i)}=f_k({\psi }_{k-1|k-1}^{(i)})+w_k^{(i)}.$$

$$z_k=\frac{1}{M}\sum _{i=1}^M{z}_k^{(i)},$$

$$z_k^{(i)}=h_k({\psi }_{k|k-1}^{(i)})+v_k^{(i)}.$$

$${\psi }_{k|k}^{(i)}={\psi }_{k|k-1}^{(i)}+K_k(y_k-z_k^{(i)}),$$

$$K_k=U_k(V_k{)}^{-1},$$

$$U_k=\frac{1}{M-1}\sum _{i=1}^M({\psi }_{k|k-1}^{(i)}-{\psi }_{k|k-1})(z_k^{(i)}-z_k{)}^T,$$

$$V_k=\frac{1}{M-1}\sum _{i=1}^M(z_k^{(i)}-z_k)(z_k^{(i)}-z_k{)}^T.$$

$$p({\psi }_k|{\mathbf{y}}_{1:k})\simeq \frac{1}{M}\sum _{i=1}^M\delta ({\psi }_k-{\psi }_{k|k}^{(i)}).$$

$$y_k=C({\xi }_k)x_k+D({\xi }_k)u_k$$

$${\psi }_k=[x_k^T,{\xi }_k^T{]}^T.$$

$$f_k({\psi }_{k-1})=f({\psi }_{k-1},u_{k-1})$$

$$=\left[\begin{array}{c}A({\xi }_{k-1})x_{k-1}+B({\xi }_{k-1})u_{k-1}\\ {\xi }_{k-1}\end{array}\right],$$

$$h_k({\psi }_k)=h({\psi }_k,u_k)=C({\xi }_k)x_k+D({\xi }_k)u_k$$

$${\xi }_{k|k}=\frac{1}{M}\sum _{i=1}^M{\xi }_{k|k}^{(i)},$$

$${\bar{G}}^T\bar{G}=$$

$$\cdot [\mathrm{row}(A({\xi }_{k|k})),\mathrm{row}(B({\xi }_{k|k}))]).$$

$$x_{k|k-1}=\frac{1}{M}\sum _{i=1}^M{x}_{k|k-1}^{(i)}$$

$${\psi }_{k|k-1}^{(i)}={\left[x_{k|k-1}^{(i)T},{\xi }_{k|k-1}^{(i)T}\right]}^T.$$

Full text

I Introduction

The influence of stochastic noise in modeling and control is more or less inevitable in practice. In modeling of systems, for example, noise may cause a fluctuation of estimated parameters. Such a fluctuation is known to arise, e.g., in flight vehicle systems [1]. If the influence is larger than we can disregard, it would be better to take into account the presence of noise in modeling and control. As a system class handling this issue, this paper deals with discrete-time linear systems whose dynamics are determined by a stochastic process (regarded as noise). Such a system is called in the field of analytical dynamics a discrete-time linear random dynamical system [2].

Markov jump systems [3] are one of the most known special cases of random dynamical systems, whose dynamics are determined by a finite-mode Markov chain. Although the class of stochastic dynamics that can be described with a finite-mode Markov chain is limited from the viewpoint of possible noise, the study [4] extended the associated theory so that a more general Markov process can be dealt with. On the other hand, the authors are also developing another theoretical framework of controlling random dynamical systems whose dynamics are determined by a stochastic process that is independent and identically distributed (i.i.d.) with respect to the discrete time [5, 6]. Since i.i.d. processes are a special case of Markov processes, one might consider that our results could be covered by those for the Markov jump systems. However, this is not true and our studies have an advantage that we can deal with unbounded supports for coefficient random matrices (depending on an i.i.d. process) in the system model, which are difficult to deal with in the case of Markov jump systems. The class of such coefficient random matrices includes those having, e.g., normally distributed entries. If the process determining the system dynamics can be seen as an i.i.d. process, and if we do not have sufficient prior information on its support (i.e., the range of fluctuation), then using our framework may be a good option.

In [5, 6], we have discussed state feedback controller synthesis for stabilization of the systems. However, the situation where all the system states can be directly measured is considered to be unusual in practice. In addition, the distribution of the stochastic process determining system dynamics may be also unknown, which is required in designing the state feedback gain with our theory. The purpose of this paper is to resolve these issues by exploiting a sequential Monte Carlo method called the ensemble Kalman filter (EnKF) [7, 8] in distribution modeling as well as online state observation. The EnKF can be used for solving nonlinear filtering problems. The computational cost of the EnKF is relatively not expensive, and it has been exploited in large-scale problems such as atmospheric data assimilation [9]. We have a hope that statistical methods such as the EnKF would be compatible with our stochastic systems control theory.

This paper is organized as follows. Section II describes discrete-time linear random dynamical systems whose dynamics are determined by an i.i.d. process, and then, briefly reviews the earlier result on synthesis of a stabilizing state feedback gain for the system. Section III introduces a standard usage of the EnKF in parameter estimation, and discusses an EnKF-based method of modeling the distribution behind the system, which is used for the above gain design. Section IV proposes an output feedback control framework using the EnKF and the feedback gain designed with the reviewed result. The effectiveness of such a framework is demonstrated in Section V through numerical experiments.

We use the following notation in this paper. The set of real numbers, that of positive real numbers and that of non-negative integers are denoted by ${\bf R}$, ${\bf R}_{+}$ and ${\bf N}_{0}$, respectively. The set of $n$-dimensional real column vectors and that of $m\times n$ real matrices are denoted by ${\bf R}^{n}$ and ${\bf R}^{m\times n}$, respectively. The set of $n\times n$ positive definite matrices is denoted by ${\bf S}^{n\times n}_{+}$. The identity matrix of size $n$ is denoted by $I_{n}$. The Euclidean norm is denoted by $||\cdot||$. The vectorization of a matrix in the row direction is denoted by ${\rm row}(\cdot)$, i.e., ${\rm row}(\cdot):=[{\rm row}_{1}(\cdot),\ldots,{\rm row}_{m}(\cdot)]$, where $m$ is the number of rows of the matrix and ${\rm row}_{i}(\cdot)$ denotes the $i$th row. The Kronecker product is denoted by $\otimes$. The (block) diagonal matrix is denoted by ${\rm diag}(\cdot)$. The Dirac delta function is denoted by $\delta(\cdot)$. The expectation (i.e., the expected value) of a random variable is denoted by $E[\cdot]$; this notation is also used for the expectation of a random matrix. If $s$ is a random variable obeying the distribution ${\cal D}$, then we represent it as $s\sim{\cal D}$.

II Discrete-Time Linear Random

Dynamical Systems and State Feedback Stabilization

II-A *Discrete-Time Linear Random Dynamical Systems with

an i.i.d. Process*

Let us consider the $Z$-dimensional discrete-time stochastic process $\xi=\left(\xi_{k}\right)_{k\in{\bf N}_{0}}$ satisfying the following assumption.

Assumption 1

$\xi_{k}$* is independent and identically distributed (i.i.d.) with respect to the discrete time $k\in{\bf N}_{0}$.*

This assumption naturally makes $\xi$ stationary and ergodic [10]. The support of $\xi_{k}$ is denoted by ${\boldsymbol{\mit\Xi}}$. By definition, ${\boldsymbol{\mit\Xi}}\subset{\bf R}^{Z}$, and ${\boldsymbol{\mit\Xi}}$ corresponds to the set of values that $\xi_{k}$ can take (at each $k$).

With such a process $\xi$, consider the discrete-time linear random dynamical system

[TABLE]

where $x_{k}\in{\bf R}^{n}$, $u_{k}\in{\bf R}^{m}$, $A:{\boldsymbol{\mit\Xi}}\rightarrow{\bf R}^{n\times n}$, $B:{\boldsymbol{\mit\Xi}}\rightarrow{\bf R}^{n\times m}$, and the initial state $x_{0}$ is assumed to be deterministic. The representation of the coefficient random matrices in (1) is general in the sense that any random matrices (denoted by $A_{k}$ and $B_{k}$) can be represented in the form with appropriate $\xi$; we can always take $A(\cdot)$, $B(\cdot)$ and $\xi$ such that $A(\xi_{k})=A_{k}$ and $B(\xi_{k})=B_{k}$ (under Assumption 1). Our control approach is developed for such a system.

II-B Stability and State Feedback Stabilization

If full prior information about the distribution of $\xi_{0}$ as well as the structure of $A(\cdot)$ and $B(\cdot)$ (involving constant coefficients) is available, and if the exact value of the system state can be obtained at each time step, then we can readily use our earlier results in [6] for stabilizing system (1). This subsection briefly reviews the earlier results under such an ideal situation, before proceeding to the case with a lack of information.

Let us consider the state feedback

[TABLE]

with the static time-invariant gain $F\in{\bf R}^{m\times n}$, and the associated closed-loop system

[TABLE]

To define second-moment stability for this system, we introduce the following assumption (which is actually a minimal requirement for the definition, although the details are omitted).

Assumption 2

The squares of elements of $A(\xi_{k})$ and $B(\xi_{k})$ are all Lebesgue integrable, i.e.,

[TABLE]

where $A_{ij}(\xi_{k})$ and $B_{ij}(\xi_{k})$ represent the $(i,j)$-entries of $A(\xi_{k})$ and $B(\xi_{k})$, respectively.

For each fixed $F$, this assumption ensures the squares of elements of $A_{\rm cl}(\xi_{k})$ are also Lebesgue integrable. Under Assumptions 1 and 2, we define exponential second-moment stability (i.e., exponential mean square stability) [11] as follows.

Definition 1

The system (3) with a fixed $F$ satisfying Assumptions 1 and 2 is said to be exponentially stable in the second moment if there exist $a\in{\bf R}_{+}$ and $\lambda\in(0,1)$ such that

[TABLE]

This stability notion can be characterized by a Lyapunov inequality as in the following theorem [6].

Theorem 1

Suppose the open-loop system (1) satisfies Assumptions 1 and 2. For given $F\in{\bf R}^{m\times n}$, the following two conditions are equivalent.

The closed-loop system (3) is exponentially stable in the second moment. 2. 2.

There exist $P\in{\bf S}^{n\times n}_{+}$ and $\lambda\in(0,1)$ such that

[TABLE]

Based on the Lyapunov inequality (6), we can further obtain the following theorem [6], which gives an inequality condition for designing a stabilizing state feedback gain.

Theorem 2

Suppose the open-loop system (1) satisfies Assumptions 1 and 2. There exists a gain $F$ such that the closed-loop system (3) is exponentially stable in the second moment if and only if there exist $X\in{\bf S}^{n\times n}_{+}$, $Y\in{\bf R}^{m\times n}$ and $\lambda\in(0,1)$ satisfying

[TABLE]

for $\bar{G}^{\prime}_{A}$ and $\bar{G}^{\prime}_{B}$ given by

[TABLE]

with a matrix $\bar{G}\in{\bf R}^{\bar{n}\times(n+m)n}\ (\bar{n}\leq(n+m)n)$ satisfying

[TABLE]

In particular, $F=YX^{-1}$ is one such stabilizing gain.

While the Lyapunov inequality (6) involved decision variables contained in the expectation operation, the decision variables in (7) are all uncontained in the expectation operation. Hence, once $\bar{G}^{\prime}_{A}$ and $\bar{G}^{\prime}_{B}$ are calculated, (7) can be solved as a standard linear matrix inequality (LMI) (for each fixed $\lambda$). By minimizing $\lambda$ through a bisection with respect to $\lambda^{2}$ under such an LMI, we can obtain a stabilizing gain that is optimal in the sense of exponential second-moment stability; the minimal $\lambda$ corresponds to that in (5), which corresponds to the convergence rate of $\sqrt{E[||x_{k}||^{2}]}$ with respect to $k$.

As reviewed above, our earlier results can be readily used when full information of system (1) (except the value of $\xi_{k}$) is available. However, it is not always possible to access full information in practical control problems, and we may have to achieve some required performance only with limited available information. To exploit our control theory even in such a situation, the following sections consider using the ensemble Kalman filter [7, 8] in modeling of unknown system parts as well as online estimation of the system state.

III Distribution Modeling Using Ensemble Kalman Filter

III-A Ensemble Kalman Filter

The ensemble Kalman filter (EnKF) is a sequential Monte Carlo method that solves nonlinear filtering problems. The problem of estimating system states and coefficient matrix parameters for deterministic linear systems can be seen a special case of nonlinear filtering problems, and the EnKF is known to be useful for such estimation. Before proceeding to the arguments for our random dynamical systems, let us briefly review a standard usage of the EnKF in this subsection (see, e.g., [8] for further details).

As a simulation model in the EnKF, we consider the state equation

[TABLE]

and the observation equation

[TABLE]

with given nonlinear functions $f_{k}(\cdot)$ and $h_{k}(\cdot)$, where $\psi_{k}$, $z_{k}$, $v_{k}$ and $w_{k}$ are the state, the observation, the system noise and the observation noise, respectively. The noise inputs are assumed to be independent of each other and obey the following multi-dimensional zero-mean normal distributions with given covariance matrices $Q$ and $R$.

[TABLE]

In the filtering problem, $\psi_{k}$ is regarded as a random vector, which corresponds to an estimate of the system state and unknown parameters111Unknown parameters are also dealt with as the state of the simulation model in the parameter estimation case. of a plant. In particular, the EnKF uses a set $\{\psi^{(i)}_{k|k}\}_{i=1,\ldots,M}$ of sample values called an ensemble for approximating the distribution of $\psi_{k}$; each ensemble member is simulated with the above simulation model. In this paper, we consider approximating the density function $p(\psi_{k}|{\bf y}_{1:k-1})$ of the prior probability distribution as

[TABLE]

where ${\bf y}_{1:k}:=(y_{i})_{i=1,\ldots,k}$ is the measurement (i.e., the output of the real plant) up to $k$, and $\psi^{(i)}_{k|k-1}$ is the forecast of state $\psi_{k}$ calculated with the (filtered) $i$th ensemble member $\psi^{(i)}_{k-1|k-1}$ at $k-1$. Under the approximation (15), an EnKF is given as follows, where $w_{k}^{(i)}$ and $v_{k}^{(i)}\ (i=1,\ldots,M;k\in{\bf N})$ are independent and their distributions are given in a fashion similar to (14) (with common $Q$ and $R$ for $i$).

State prediction (state update):

[TABLE]

Output prediction (output update):

[TABLE]

Filtering:

[TABLE]

In the above EnKF, $K_{k}$ corresponds to a Kalman gain. With this filter, the ensemble $\{\psi^{(i)}_{k|k}\}_{i=1,\ldots,M}$ approximates the density function $p(\psi_{k}|{\bf y}_{1:k})$ of the posterior probability distribution as

[TABLE]

Through a sequential use of (16)–(23) with respect to $k$, the above posterior distribution gives an estimate of the system state and unknown parameters.

Remark 1

In this subsection, we introduced a standard usage of the EnKF. However, this is not a unique usage, and there are varieties; using EnKFs is also not a unique option for nonlinear filtering. For more details about the varieties and other approaches, see [7, 8, 12] and other sophisticated articles. Since discussing the differences among them is beyond the scope of this paper, we only deal with the above type of EnKF.

III-B Distribution Modeling

In most practical control problems, the state of a plant cannot be directly measured, and some observation method is needed when we use state feedback controllers. The same comment also applies to our random dynamical systems, as stated in Section II. Hence, we consider the following output equation in addition to the state equation (1) for the system model222We consider two models of a plant in this paper. One is the simulation model used in the EnKF and the other is the system model for representing a plant in controller synthesis..

[TABLE]

We assume that $x_{k}$ (and $\xi_{k}$) cannot be measured directly, and only $y_{k}$ is available at each $k$. The goal of this paper is to control the plant modeled by (1) and (25) with our control approach reviewed in Section II and an EnKF. In particular, we consider the situation where the prior information on the distribution of $\xi_{k}$ is unavailable, and regard $\xi_{k}$ as a part of the state in the simulation model (12) and (13), i.e.,

[TABLE]

This is reasonable since the prior information on the accurate distributions of fluctuating parameters is less usual to be provided than that on constant parameters in practical problems. To design a gain $F$ in (2) in such a situation, this subsection first uses the EnKF for modeling the distribution of $\xi_{k}$.

In this subsection, we temporarily see $u_{k}$ as an input for identification. To reflect the influence of $u_{k}$ also in the simulation model of the EnKF, we consider in (17) and (19) $f_{k}(\cdot)$ and $h_{k}(\cdot)$ such that

[TABLE]

under (26). Obviously, these functions are nonlinear in $\psi_{k}$ (even when $A$, $B$, $C$ and $D$ are linear in $\xi_{k}$). Here, if $\xi_{k}$ is supposed to be an unknown deterministic parameter vector, the usual deterministic parameter estimation can be carried out through the corresponding EnKF. In practical parameter estimation problems, however, the estimated value of a parameter more or less fluctuates even after sufficient time has elapsed. If the fluctuation is larger than we can disregard, it might be better to take account of it in modeling and control. This is why we deal with representation (1) and (25) with random $\xi_{k}$ as the system model.

Under Assumption 1 on $\xi$, the time sequence of the estimated value of $\xi_{k}$ obtained by the EnKF could be seen as giving an empirical distribution of $\xi_{k}$. This is our basic idea of modeling random $\xi_{k}$. To make the point clearer, we introduce the notation $\xi_{k|k}^{(i)}$ for representing the parameter part of the filtered $i$th ensemble member $\psi_{k|k}^{(i)}$ defined in a way similar to (26). By inputting a signal for identification such as maximal length sequences to a plant (modeled by system (1) and (25)) as $u_{k}$, the EnKF connected to the plant provides us with the time sequence of $\xi_{k|k}^{(i)}$ (for $i=1,\ldots,M$). Then, one of the reasonable estimates of $\xi_{k}$ can be given by

[TABLE]

which corresponds to the mean of the parameter-related part of the posterior distribution (24). If $\xi_{k|k}$ can be seen as being i.i.d. with respect to $k$ in association with Assumption 1, the set $\{\xi_{k|k}\}_{k=K_{0},\ldots,K_{0}+K_{1}}$ with sufficiently large integers $K_{0}$ and $K_{1}$ corresponds to an empirical distribution of the estimate of $\xi_{k}$. Although obtaining the exact value of $\xi_{k}$ at each $k$ is generally difficult even with EnKFs, we have a hope that the obtained empirical distribution would not be totally different from the original distribution (if such an original actually exists), and decided to employ it as the modeling result of the distribution of $\xi_{k}$. If we have the prior information about the class of the distribution of $\xi_{k}$ (such as being a normal distribution), the empirical distribution could be further exploited for determining its hyperparameters.

Remark 2

Not only the information about the distribution with respect to $k$ but also that with respect to $i$ for $\xi_{k|k}^{(i)}$ might be useful in modeling the distribution of $\xi_{k}$. However, using only the information about $i$ at some fixed $k$ is unreasonable since $\xi_{k|k}^{(i)}\ (i=1,\ldots,M)$ are distributed mainly because of the artificial noises $w_{k}^{(i)}$ and $v_{k}^{(i)}$ in the EnKF.

Remark 3

Parameter estimation is usually carried out for stable systems. Hence, the plant in this section basically has to be stable (in a stochastic sense). One might consider this contradicts our tackling a stabilization problem. However, this is a sort of dilemma that is in common with conventional deterministic systems modeling and control. In most practical control problems, we usually have to first stabilize plants without their models by using, e.g., PID controllers, if they seem to be unstable. After such stabilization, we identify the systems, and design controllers with sophisticated model-based control theories to achieve required performance. A similar comment also applies to our problem. Although this paper is the first attempt at combining a nonlinear filtering approach with our recent results, and hence, only the convergence rate $\lambda$ can be dealt with as the control performance index at this moment, other performance indices such as $H_{2}$ norm may be available in the future, as is the case with deterministic systems control [13].

IV Closed-Loop System with Ensemble Kalman Filter and

Stabilization Control

With the modeling result of the distribution of $\xi_{k}$ in the preceding section, the expectation in (11) can be calculated. That is, we have

[TABLE]

Hence, a state feedback gain $F$ that minimizes $\lambda$ in (5) for the ideal closed-loop system (3) with such modeled $\xi_{k}$ can be obtained by minimizing $\lambda$ with respect to the inequality condition in Theorem 2. Although the state feedback gain cannot be directly used for the state-unmeasurable system (1) and (25), it may work well if the state is observed (i.e., sequentially estimated) from the measurable output with satisfactory accuracy333The principle of separation of estimation and control has not been confirmed to hold for our problem at this moment.. The EnKF can be used also as a state observer, and this section discusses the closed-loop control framework using it.

The algorithm of the EnKF for state observation is similar to that described in the preceding section. Only an essential difference is that we have to carry out online estimation of the state in the closed-loop control, which was needless in the preceding offline modeling. In the modeling, we used (the $\xi$-related part of) the filtered ensemble members $\psi_{k|k}^{(i)}$ for obtaining the distribution model of $\xi_{k}$. The filtered ensemble members are, however, calculated with $y_{k}$ as in (20), which is the output of the plant at $k$. Since this $y_{k}$ is the result of inputting $u_{k}$ on the plant at $k$, and since $u_{k}$ has to be determined with (an estimate of) the state at $k$, using (the $x$-related part of) $\psi_{k|k}^{(i)}$ for online estimation is generally difficult. Hence, we use the forecasts $\psi_{k|k-1}^{(i)}\ (i=1,\ldots,M)$ for online estimation of $x_{k}$. More precisely, we regard

[TABLE]

as the estimate of $x_{k}$, where $x_{k|k-1}^{(i)}$ denotes the $x$-related part of $\psi_{k|k-1}^{(i)}$, i.e.,

[TABLE]

Here, since we already obtained a model of the distribution of $\xi_{k}$, one might consider it needless to estimate $\xi_{k}$ simultaneously with $x_{k}$ as in (32). However, we deal with this $\psi_{k|k-1}^{(i)}$ even in the present state observation. This is mainly because we cannot obtain the exact value of $\xi_{k}$ in the simulation model (12) and (13) even if we take a sample from the distribution model at each $k$; what is worse, the sample is generally independent of true $\xi_{k}$, which might cause an unnecessary deterioration of the accuracy of the estimation of $x_{k}$. Another important reason is the hedge against possible modeling errors.

With the above idea of state observation, our control framework with an EnKF can be given as follows.

State prediction: (16) and (17) with $f_{k}(\cdot)=f(\cdot,u_{k-1})$ such that (26) and (27).

Control input decision: (31), (32) and

[TABLE]

Output prediction: (18) and (19) with $h_{k}(\cdot)=h(\cdot,u_{k})$ such that (26) and (28).

Filtering: (20)–(23).

To facilitate further understanding, we also provide Fig. 1 for showing the signal flows between these four functions (constituting the whole controller) as well as their connection with the plant. Since the computational cost of the EnKF is not expensive compared to other heuristic nonlinear filtering approaches, and since the internal control law (33) is very simple, this controller would work with a practical sampling rate in various control problems. In addition, if the EnKF provides a satisfactorily accurate estimate of the state online, the closed-loop system is expected to be stabilized by Theorem 2.

V Numerical Experiments

V-A Application to Partially Unknown System

Consider the three subsystems described by

[TABLE]

where $s_{1}$ and $s_{2}$ are i.i.d. processes that are independent of each other and distributed as

[TABLE]

for each $k\in{\bf N}_{0}$. Let us further consider the networked system described in Fig. 2 consisting of the above subsystems, where $u$ is the input and $y_{k}=[y^{\prime}_{1k},y^{\prime}_{2k},y^{\prime}_{3k}]^{T}$ is the measurable output. This networked system is dealt with as the plant to be controlled in this section. With $q_{k}=[x^{\prime T}_{1k},x^{\prime T}_{2k},x^{\prime T}_{3k}]^{T}$, the plant can be described as follows.

[TABLE]

For the above plant, we assume the situation where its prior information is partially unknown and the system model is given by (1) and (25) with coefficients

[TABLE]

where $\xi_{k}=[\xi_{1k},\xi_{2k},\xi_{3k}]^{T}$ represents the unknown part of the model; $\xi_{1k}$, $\xi_{2k}$ and $\xi_{3k}$ are time-varying parameters introduced correspondingly to $0.7s_{1k}$, $s_{2k}$ and 0.6, respectively. With such a system model, we consider stabilizing the plant (37) and (38) through our approach.

We first use our modeling method discussed in Subsection III-B for determining the distribution of $\xi_{k}$ in the system model. The simulation model for the EnKF is given by (26)–(28) with the coefficients (39) of the system model. We temporarily regard $u$ as the input for identification at this modeling stage, and use a maximal length sequence for the input. In addition, we take $M=300$, $R=0.01$ and $Q=0.01{\rm diag}(1,1,1,1,1,1,2,2,2)$ for EnKF parameters. Then, in our computation, we obtained the response shown in Fig. 3 for a sample path of $s_{1}$ and $s_{2}$, where $q_{0}=[1,1,1,1,1,1]^{T}$ and $\psi_{0}=[x_{0}^{T},\xi_{0}^{T}]^{T}=0$ were used as the initial values. Since the response of the states did not diverge, the plant seems not to be unstable. With the data of the filtered response from $k=20$ through $k=100$ (i.e., $K_{0}=20$ and $K_{1}=80$), we obtained an empirical distribution of $\xi_{k}$. The corresponding sample mean was

[TABLE]

and the covariance was

[TABLE]

We can see that the above sample mean is not very far from the expected value $[0.56,0.2,0.6]^{T}$ of $[0.7s_{1k},s_{2k},0.6]^{T}$, and a similar comment also applies to the sample variance of $\xi_{2k}$, which corresponds to $s_{2k}$ ($0.0750$ is not very far from $0.3^{2}=0.09$). Compared to these values, the sample covariance of $[\xi_{1k},\xi_{3k}]^{T}$ seems not close to the covariance of $[0.7s_{1k},0.6]^{T}$. This might imply that the information obtained from the output of the plant is not sufficient for distinguishing the values of $\xi_{1k}$ and $\xi_{3k}$ in the simulation model, which may be coupled in a single subsystem part (i.e., the $G_{2}$ part). Our purpose of obtaining the distribution of $\xi_{k}$ is, however, not to obtain the exact model of the plant but to achieve high performance in control by taking account of randomness behind the plant. Fortunately, the predicted/filtered $x_{k}$ in Fig. 3 seems to follow $q_{k}$ with a high degree of accuracy, and the above modeling error might not be a problem in characterizing the randomness. Indeed, the system model with the above mean and covariance leads us to a successful result at the next synthesis stage.

With the above system model, we calculated $\bar{G}^{\prime}_{A}$ and $\bar{G}^{\prime}_{B}$ in (8)–(11) (recall (30)), and solved (7) so that the obtained state feedback gain $F=YX^{-1}$ can achieve a minimal $\lambda$ satisfying (5) for the ideal closed-loop system (3). Then, we obtained ($\bar{n}=4$,) $\lambda=0.8432$ and

[TABLE]

With this gain, we constructed the controller described in Fig. 1, and applied it to the present plant. Then, for the external input

[TABLE]

added on $u_{k}$, we obtained the response of the state $q_{k}$ of the plant as shown in Fig. 4, where this response was calculated through approximating the expectation by the sample mean obtained with 200 sample paths of $s_{1}$ and $s_{2}$ (as well as those of $v^{(i)}$ and $w^{(i)}$ in the EnKF). The EnKF parameters and the initial values were the same as those in the above modeling. As we can see, the state after $k=50$ successfully decreases. Although the decay rate after $k=90$ is not low compared to that before $k=90$, this would be an effect of the artificially introduced additive noises $v^{(i)}$ and $w^{(i)}$ in the EnKF. If we do not use any feedback control for the present plant, the response of $q_{k}$ is as shown in Fig. 5, where the same environment and condition (including the external input and samples) as the above closed-loop case was used for this computation. It is obvious from these results that our control approach successfully achieved the improvement of the degree of stability (in the sense of the convergence rate of $\sqrt{E[||q_{k}||^{2}]}$). The convergence rate of $\sqrt{E[||q_{k}||^{2}]}$ calculated with the date at $k=60$ and $k=80$ in Fig. 4 was $0.8473$, which is very close to our theoretical result of $\lambda=0.8432$; possible reasons for the small gap between these values are the influences of modeling errors, our use of an EnKF in control, and the sample-base approximation in response computation.

Remark 4

The computation of a single response (i.e., a single sample path) of $q_{k}$ in our closed-loop system from $k=0$ to $k=100$ took from 0.38 s to 0.51 s with MATLAB R2017b running on a laptop equipped with 8.00 GB RAM and Intel(R) Core(TM) i7-5600U CPU @ 2.60 GHz. Hence, the average computation time for one cycle of our feedback control was less than 10 ms for the present example. This implies that our controller with an EnKF could work in various actual problems without high-spec computers. (Although we used $M=300$ as the size of ensemble, this size is actually conservative and can be reduced for the present example.)

V-B *Comparison with Deterministic Systems Control Approaches

Combined with EnKF*

In the preceding subsection, we confirmed that our control approach can improve the control performance (although only stability degree was dealt with at this moment). This was actually achieved by our appropriately taking account of randomness behind the plant in our synthesis using Theorem 2. To further confirm this, we next compare our control approach with those using the EnKF and conventional deterministic systems control methods. That is, we consider the cases where the feedback gain in (33) is designed with other conventional methods.

If we discard the information on the covariance (41) and take the standpoint that the mean (40) is the true constant value of $\xi_{k}$ in (39), then the system model with the corresponding coefficient matrices becomes a standard linear time-invariant (LTI) system. Hence, we can use various conventional design methods for the LTI system. As an example, we first consider using the pole placement technique. The design result of achieving the poles $(0.9,0.8,0.7,0.6,0.5,0.4)$ is given by

[TABLE]

It is obvious from the poles that this gain can stabilize the ideal closed-loop system (3) with the above LTI system. However, it failed to stabilize the actual closed-loop system with the structure in Fig. 1 (despite the stable behavior of the open-loop plant in Fig. 5). A similar comment also applies to the cases with other poles in our computation. Hence, the pole placement technique cannot be used for the present example. One might consider this would be because of our using not Theorem 2 (i.e., Theorem 1) but other methods for the LTI system; indeed, Theorem 2 can be actually used also for LTI systems as a special case. However, the influence of discarding the information on the covariance is fatal even in this direction; we obtained with Theorem 2 a feedback gain with $\lambda=0.1143$ (which is optimal for the ideal LTI closed-loop system) that failed to stabilize the actual closed-loop system. As is confirmed from these results, using only the mean values of the parameters in the system model may lead us to false results for systems with stochastic dynamics.

To circumvent this kind of issues, robust control methods [14, 15] have been conventionally studied, in which variations of parameters such as those in Fig. 3 are viewed as deterministic uncertainties. However, the usual use of such methods in the present example is actually impossible because the variation range of the parameters is too large; this makes it difficult to obtain a gain that is theoretically guaranteed to stabilize even the ideal closed-loop system. If the variation range of the parameters used for the synthesis is narrowed, one could obtain some feedback gain. However, such synthesis has an issue similar to that in the above nominal stabilization case (the closed-loop system is unstable or the convergence rate is not improved). This difficulty in robust control becomes serious as the number of the parameters to be estimated increases.

VI Conclusions

In this paper, we extended our state feedback control approach for random dynamical systems toward output feedback control using an EnKF. The EnKF was used also for modeling the distributions behind the systems. The results of our numerical experiments demonstrated the potentiality of our approach in improving the control performance through taking account of the randomness behind the systems. In particular, it might be surprising that the theoretical guarantee on the control performance obtained for the ideal closed-loop system (i.e., without an EnKF) provided us with a good estimate of that for the actual closed-loop system using the EnKF even in the situation where the given model of the plant was partially unknown. Although we directly dealt with a random dynamical system as the plant in numerical experiments, our approach is expected to be widely effective in practical problems in which the estimates of parameters fluctuate more than we can disregard.

Since the distribution modeling using the EnKF was successful (in the sense of achieved control performance), and since an estimate of the unknown parameters can also be obtained during the feedback control (recall (32)), we might be able to exploit the estimate to adjust the feedback gain online. That is, a sort of learning type control could be realized in which the controller estimates the distribution behind the system online and determines the internal feedback gain with the estimate by itself. Our earlier results [5, 16] using random polytopes [17] would be useful for this direction of studies; random polytopes can describe, e.g., the variation in mean and variance of random matrices. The present paper actually also has the role of the first step toward such advanced control.