Comparative Study of Markov Chain Filtering Schemas for Stabilization of Stochastic Systems under Incomplete Information

Abstract

The object under investigation is a controllable linear stochastic differential system affected by some external statistically uncertain piecewise continuous disturbances. They are directly unobservable but assumed to be a continuous-time Markov chain. The problem is to stabilize the system output concerning a quadratic optimality criterion. As is known, the separation theorem holds for the system. The goal of the paper is performance analysis of various numerical schemes applied to the filtering of the external Markov input for system stabilization purposes. The paper briefly presents the theoretical solution to the considered problem of optimal stabilization for systems with the Markov jump external disturbances: the conditions providing the separation theorem, the equations of optimal control, and the ones defining the Wonham filter. It also contains a complex of the stable numerical approximations of the filter, designed for the time-discretized observations, along with their accuracy characteristics. The approximations of orders $\frac{1}{2}$, 1, and 2 along with the classical Euler–Maruyama scheme are chosen for the comparison of the Wonham filter numerical realization. The filtering estimates are used in the practical stabilization of the various linear systems of the second order. The numerical experiments confirm the significant influence of the filtering precision on the stabilization performance and superiority of the proposed stable schemes of numerical filtering.

1. Introduction

The theoretical solution to the optimal filtering problem of stochastic differential system states with the development of its effective numerical solution is a traditional field of extensive research [1,2,3,4]. The online calculation of the high precise estimates of the system state given the noisy indirect observations is in demand for navigation [5,6,7], telecommunications [8,9,10], finance [11,12,13,14,15], processing in micro- [16] and nano-structures [17], and many other purposes [18,19,20,21,22,23,24]. However, advanced filtering algorithms play a crucial role in the synthesis of the control in the case of incomplete information [25,26,27,28,29]. The spectrum of such applications is rather broad [30,31,32,33].

When the optimal control depends on the available observation through the optimal estimate, fulfillment of the separation theorem is presumed [34,35]. The corresponding “good lucks” are not numerous and primarily include the case of the Linear Quadratic Gaussian (LQG) systems [36]. When the separation theorem fails for some reason, one can use the so-called “separation principle”, i.e., find the conditional optimum within the class of strategies, which are the functions of the state estimate, not of the original observations. So, the availability of a high precise filtering estimate is a significant factor for the successful solution to numerous control problems under incomplete information.

As is known, the optimal in the mean-square sense state estimate coincides with a conditional mathematical expectation of the state given the available observations. The estimates, which are optimal in the sense of criteria other than mean-square one, are the functionals of the corresponding conditional distribution. This, in turn, is a solution to the Kushner–Stratonovich or Zakai partial stochastic differential equations (SDEs) [29,37,38,39,40]. There are only a few cases when the conditional distribution is a solution to some finite-dimensional closed system of SDSs. The Kalman–Bucy filter [41,42], the Wonham one [2,25], and some particular cases [43] belong to this class. Development of the stable numerical realizations of the optimal filters is nontrivial even for the linear Gaussian observation systems [44].

The aim of this paper is an investigation of how the filtering performance impacts the quality of the subsequent optimal control with incomplete information. As a testbed, we choose the stabilization of a linear stochastic differential system affected by some statistically uncertain external disturbance. This represents a continuous-time Markov chain (CTMC). The theoretical solution to the stabilization problem is presented in [33], and the separation theorem is proved. The optimal control strategy needs the optimal estimate, and the last, in turn, is determined by the Wonham filter. The numerical realization of the filter is nontrivial: the classical numerical schemes such as the Euler–Maruyama one [45,46] are unstable in this case. The point is that they do not provide non-negativity and normalization properties for the estimate approximations. The authors of [47] suggest a new concept of the stable numerical approximations of the filtering estimates. It is left to solve the Wonham filter SDE numerically for direct development of the optimal filtering algorithm by the observations discretized by time. The corresponding estimate is calculated recursively by a variant of the Bayes formula. The formula contains the integrals, which are the shift-scale mixtures of Gaussians, and one cannot calculate them analytically. The author of [48] suggests various numerical schemes of their calculation and the accuracy characteristics. This paper presents a comparison of the numerical filtering schemas in light of their influence on the performance of the subsequent stabilization.

The paper is organized as follows. Section 2 presents the investigated control system affected by the external statistically uncertain piece-wise constant Markov disturbances. We state the problem of optimal stabilization and present the equations describing the optimal stabilization strategy. The separation theorem is valid for the system, so the optimal strategy depends on the observation via the optimal filtering state estimate.

Section 3 introduces the mechanical system, an overhead crane, which serves as a prototype of the controlled stochastic system used for the comparative numerical analysis. We present the general properties of the system, its evolution with and without external disturbances, and the influence of the criterion coefficients on the control character. Section 4 is devoted to the numerical realizations of the Wonham filter. Based on the optimal filtering estimates of the CTMC given the continuous-time observations, we introduce the numerical schemas of its approximation of the order $\frac{1}{2}$, 1, and 2 given the time-discretized observations.

Section 5 plays a vital role in the presentation. It contains the results of the numerical experiments, illustrating the influence of the estimation stability on the final performance of the proposed control. We compare the suggested schemas with the classical Euler–Maruyama one. Section 5.1， Section 5.2 and Section 5.3 present the results of the stabilization of the state in the case of the stable, semistable, and unstable control system. All the examples confirm the superiority of the proposed schemas of numerical filtering for subsequent stabilization purposes. Section 6 contains concluding remarks.

2. Problem Formulation and Optimal Control Equations

Below, we briefly present the control problem and its theoretical solution (see [33] for details).

On the canonical space with filtration $\left(\mathsf{\Omega },\mathcal{F},\mathcal{P},{\mathcal{F}}_t\right)$, $t\in \left[0,T\right],$ let us consider the linear SDE

$$d{z}_t=a_t{y}_{t}dt+b_t{z}_{t}dt+c_t{u}_{t}dt+{\sigma }_{t}d{w}_t,z_0=Z.$$

(1)

Here

• $z_t\in {\mathbb{R}}^{n_z}$ is a controllable output stochastic process of the second order;
• $y_t$ is an uncertain stochastic input which represents a CTMC with the finite-state space $\left\{e_1;\dots ,e_{n_y}\right\}$ of the unit vectors in the Euclidean space ${\mathbb{R}}^{n_y}$, with the transition intensity matrix ${\mathsf{\Lambda }}_t$ and initial distribution ${\pi }_0=\mathbb{E}\left\{y_0\right\}$;
• $w_t\in {\mathbb{R}}^{n_w}$ is a standard Wiener process: $w_t,y_t,z_0$ are mutually independent;
• $u_t\in {\mathbb{R}}^{n_u}$ is a control, which is a process with a finite second moment;
• $a_t\in {\mathbb{R}}^{n_z\times n_y},b_t\in {\mathbb{R}}^{n_z\times n_z},c_t\in {\mathbb{R}}^{n_z\times n_u},{\sigma }_t\in {\mathbb{R}}^{n_z\times n_w}$ are known nonrandom matrix-valued functions.

The admissible control $u_t=U_t\left(z_t\right)$ represents any function of the output $z_t$. The optimal control $U_t=U_t\left(z\right),z\in {\mathbb{R}}^{n_z}$ minimizes the objective function

$$J\left(U_0^T\right)=\mathbb{E}\left\{{\int }_{{}_0^{}}^{{}_{}^T}{\left|P_t{y}_t+Q_t{z}_t{+R}_t{u}_t\right|}_{}^{2}dt+{\left|P_T{y}_T+Q_T{z}_T\right|}_{}^2\right\},$$

(2)

where $U_0^T=\{U_t\left(z\right),0\le t\le T\},P_t\in {\mathbb{R}}^{n_J\times n_y},Q_t\in {\mathbb{R}}^{n_J\times n_z},R_t\in {\mathbb{R}}^{n_J\times n_u},0\le t\le T$ are known bounded matrix-valued functions, $x^{\prime }$ is the x transpose, $\left|x\right|=\sqrt{x^{\prime }x}$ is the Euclidean norm.

Under the standard conditions of (1) and (2) (piecewise continuity and uniform degeneracy of $R_t^{{}^{\prime }}{R}_t$ and ${\sigma }_t{\sigma }_t^{{}^{\prime }}$), the solution to the optimization problem

$${\widehat{U}}_0^T=\left\{{\widehat{U}}_t\left(z\right),0\le t\le T\right\}\in \mathrm{argmin}J\left(U_0^T\right)$$

takes the form

$${\widehat{u}}_t=-\frac{1}{2}{\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}\left(c_t^{{}^{\prime }}\left(2{\alpha }_t{\widehat{z}}_t+{\beta }_t{\widehat{y}}_t\right)+2R_t^{{}^{\prime }}\left(P_t{\widehat{y}}_t+Q_t{\widehat{z}}_t\right)\right),$$

(3)

where

$$\frac{d{\alpha }_t}{dt}-\left(M_t^{\alpha }{\alpha }_t+{\alpha }_t^{{}^{\prime }}{\left(M_t^{\alpha }\right)}^{{}^{\prime }}\right)+N_t^{\alpha }-{\alpha }_t^{{}^{\prime }}{c}_t{\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{c}_t^{{}^{\prime }}{\alpha }_t=0,{\alpha }_T=Q_T^{{}^{\prime }}{Q}_T,$$

(4)

$$\frac{d{\beta }_t}{dt}+{\beta }_t{\mathsf{\Lambda }}_t^{{}^{\prime }}+M_t^{\beta }-N_t^{\beta }{\beta }_t=0,{\beta }_T=2Q_T^{{}^{\prime }}{P}_T,$$

(5)

$$\begin{array}{c}d{\widehat{y}}_t={\mathsf{\Lambda }}_t^{{}^{\prime }}{\widehat{y}}_{t}dt+\left(\mathrm{diag}\left({\widehat{y}}_t\right)-{\widehat{y}}_t{\widehat{y}}_t^{{}^{\prime }}\right)a_t^{{}^{\prime }}{\left({\sigma }_t{\sigma }_t^{{}^{\prime }}\right)}^{-1}\times \hfill \\ \hfill \times \left(d{\widehat{z}}_t-a_t{\widehat{y}}_{t}dt-b_t{z}_{t}dt-c_t{\widehat{u}}_{t}dt\right),{\widehat{y}}_0=\mathbb{E}\left\{Y\right\},\end{array}$$

(6)

where ${\widehat{z}}_t$ denotes the optimal trajectory of the state, i.e., the solution to (1) calculated for the control $u_t={\widehat{u}}_t$.

Above, the functions on the right hand side (RHS) of (4) and (5) have the form

$$M_t^{\alpha }={Q_t^{{}^{\prime }}{R}_t\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{c}_t^{{}^{\prime }},$$

$$N_t^{\alpha }=Q_t^{{}^{\prime }}\left(\mathbf{I}-R_t{\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{R}_t^{{}^{\prime }}\right)Q_t,$$

$$M_t^{\beta }=2\left(\left(a_t^{{}^{\prime }}-P_t^{{}^{\prime }}{R}_t{\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{c}_t^{{}^{\prime }}\right){\alpha }_t+P_t^{{}^{\prime }}\left(\mathbf{I}-R_t{\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{R}_t^{{}^{\prime }}\right)Q_t\right),$$

$$N_t^{\beta }=\left(Q_t^{{}^{\prime }}{R}_t+{\alpha }_t{c}_t\right){\left(R_t^{{}^{\prime }}{R}_t\right)}^{-1}{c}_t^{{}^{\prime }},$$

where $\mathbf{I}$ is an identity matrix of suitable dimensionality. So, the optimal control is a linear function of the system state $z_t$ and the optimal filtering estimate ${\widehat{y}}_t$ of the external input $y_t$. Note that ${\widehat{y}}_t$ is defined by the Wonham filter Equation (6) [2,25]. Moreover, Equation (6) determines the conditional mathematical expectation (CME) ${\widehat{y}}_t=\mathbb{E}\left\{y_t|{\mathcal{F}}_t^z\right\}$ for any control $u_t$ other than the optimal one $u_t^{\ast }$. Here, we denote the natural filtration generated by the state process $z_t\left(u\right)$, governed by the control u as ${\mathcal{F}}_t^z=\sigma \left\{z_{\tau }\left(u\right),0\le \tau \le t\right\}$, ${\mathcal{F}}_t^z\subseteq {\mathcal{F}}_t\subseteq \mathcal{F}$.

The above relations allow construction of the following stabilization algorithm

• to solve numerically ordinary differential Equations (4) and (5) (any stable numerical scheme is valid for this);
• to solve numerically SDE (6), defining the Wonham filter.

The separation theorem is valid in control of a Linear Quadratic Gaussian (LQG) system under incomplete information. The optimal control coincides with its analog under complete information, when the unobservable system state is replaced by its Kalman–Bucy filtering estimate [49]. The realization of the optimal control in the LQG case represents the numerical solution to a couple of the Riccati ODEs and a system of linear SDEs. The numerical solution of the Riccati equation can be effectively performed, for example, by the Runge–Cutta scheme, and the realization of the Kalman–Bucy filter is also well developed.

In contrast with the LQG case, the considered optimal stochastic control problem is challenging exactly in the numerical realization stage. The separation theorem is valid for the class of the investigated systems, and the optimal estimate represents the solution to the system of nonlinear SDEs. It is well known that the CME of the MJP is its conditional distribution, so the trajectories satisfy the non-negativity and normalization conditions. Visually, the instability of a numerical scheme of the Wonham filter realization looks as follows: once an estimate approximation violates the non-negativity, the subsequent trajectory of the filter “blows up” in a few steps, which makes it impossible to synthesize the control strategy. We knew this feature, so in the numerical example of [33] we chose one of the stable numerical schemes [43,47] of the Wonham filter. The point was that the example played an illustrative role only. The article can be treated as the second part of [33]: it presents the way for the practical realization of the stated optimal stochastic control problem.

3. Performance Analysis of Mechanical Actuator

Initially, for the comparative numerical study, we choose a controlled overhead crane model. The crane trolley with a hoist travels along the H-beam. The loaded trolley has significant inertia. The goal is to place it above one of the locations (for example, the railway track). The set of possible locations is known and fixed. The trolley state includes the current coordinate $x_t$ and velocity $v_t$ on the beam. The stabilization system includes a feedback velocity control block, actuated by the couple “coordinate–velocity”. In addition, the system also includes some external unobservable input $y_t$ (the signal, specifying the number of the requested railway track) and the control $u_t$ itself. So, the overhead crane model has the form

$$\begin{array}{c}d{x}_t=v_{t}dt,t\in \left(0,T\right],\hfill \\ d{v}_t=a{x}_{t}dt+b{v}_{t}dt-c{y}_{t}dt+h{u}_{t}dt+\sqrt{g}d{w}_t.\hfill \end{array}$$

(7)

The external governing signal $y_t\in \{e_1,e_2,e_3\}$ represents a homogeneous CTMC with the rransition intensity matrix (TIM) $\mathsf{\Lambda }$ and initial condition ${\pi }_0=e_1$. The constant scalars $a,b,h,g$ and row vector $\left(c_1,c_2,c_3\right)$ are known, and $w_t$ is a standard Wiener process. The term $\sqrt{g}d{w}_t$ simulates both an inner inaccuracy of the actuator and the random disturbances contaminating the governing signal. The initial condition ${(x_0,v_0)}^{{}^{\prime }}$ is a two-dimensional standard Gaussian vector.

One can see that the system (7) is stable, when $b<0$ and $b^2+4a<0$, because b and $b^2+4a$ are eigenvalues of the system matrix $\left(\begin{array}{cc}0& 1\\ a& b\end{array}\right)$.

Initially, we choose the following “basic” values for the system and control problem:

$$\begin{array}{c}a=-1,b=-0.5,T=10,h=10,g=0.01,\hfill \\ \left(c_1,c_2,c_3\right)=\left(1,0,-1\right),\mathsf{\Lambda }=\left(\begin{array}{ccc}-0.5& 0.5& 0\\ 0.5& -1& 0.5\\ 0& 0.5& -0.5\end{array}\right).\hfill & \hfill \end{array}$$

(8)

The vector $C=(\frac{c_1}{a},\frac{c_2}{a},\frac{c_3}{a})$ forms the coordinates of the railway tracks, served by the crane, while the process $C{y}_t$ represents the desired “ideal” piece-wise continuous movement of the trolley between the tracks.

We illustrate the crane’s functioning in three stages. To do this, we have to vary the parameters above to highlight various features of the estimation and/or control algorithms.

First, we consider the system without both the external governing signal $y_t$ and control $u_t$. Second, we model the system with CTMC input $y_t$ without the control $u_t$. Third, we consider the system under the action of both the input $y_t$ and control $u_t$.

We solve (7) numerically via the Euler–Maruyama scheme with the time increments $\delta =0.005,0.001,0.0001$; meanwhile, we model the CTMC $y_t$ with the smaller time step, namely ${\delta }^{small}=0.01\delta$.

At the first stage, $c_1=c_2=c_3=h=0.$

Figure 1 contains typical paths of the coordinate $x_t$ and velocity $v_t$ of the crane trolley in the case where the external input $y_t$ is absent.

The figure confirms the stable character of the system.

Second, we consider the system with the parameters (8), keeping $h=0$. The simulation is performed for the CTMCs with the two variants of the TIM:

$\mathsf{\Lambda }=\left(\begin{array}{ccc}-0.05& 0.05& 0\\ 0.005& -0.01& 0.005\\ 0& 0.5& -0.5\end{array}\right),$ and $\mathsf{\Lambda }=\left(\begin{array}{ccc}-5& 5& 0\\ 0.005& -1.005& 1.0\\ 0& 0.005& -0.005\end{array}\right).$

The first TIM provides a long staying in the initial state $e_1$ for the CTMC, while the second TIM generates the rapid transitions $e_1\to e_2$, $e_2\to e_3$ and the long staying in the last state.

The input CTMC $y_t$ influences the output $z_t$ via the process $C{y}_t$, where $C=\left(\frac{c_1}{a},\frac{c_2}{a},\frac{c_3}{a}\right)$ represent possible drift accelerations.

Figure 2 contains typical paths of the coordinate and velocity of the crane trolley affected by the CTMC with the corresponding TIMs compared with the desired switching $C{y}_t$.

One can see some stabilization of the real trolley coordinate near the “ideal” position $C{y}_t$ without additional control $u_t$, but the stabilization performance looks poor.

Third, the stabilization could be made more effective by the control $u_t$, which minimizes the functional

$$J\left(U_0^T\right)=\mathbb{E}\left\{{\int }_{{}_0^{}}^{{}_{}^T}\left({\left|C{y}_t-x_t\right|}^2+R{\left|u_t\right|}^2\right)dt\right\},$$

(9)

where $R=0.01$ characterizes the control unit cost.

Figure 3 demonstrates typical paths of the coordinate and velocity of the crane trolley with and without control $u_t$. The calculations are performed for the parameter set (8): left subplots correspond to the the controlled case with $h=10$, and right subplots stand for the uncontrolled case $h=0$.

Comparing the results demonstrated in Figure 1, Figure 2 and Figure 3, it can be seen that the control remarkably improves the stabilization performance. Without $u_t$ the system only “declares an intension” to pursue the current stabilization position $C{y}_t$. By contrast, the assisting control $u_t$ admits a more sharp reaction to the external drift change.

Fourth, we consider the system with the parameter set (8), reducing the value R, namely $R=0.00001$. The example demonstrates a potential stabilization performance in the case of low control cost.

Figure 4 illustrates the high performance of the coordinate stabilization and high amplitude oscillating character of the corresponding velocity under the optimal control $u_t$.

Note that the value R of the unit control cost does not affect the filtering performance of the external input y and is valuable only in the control procedure. Let us explain this fact from the physical point of view. The trolley has a positive mass. According to Newton’s second law, the faster the actuator must react to the external governing signal, the higher the force applied. However, each resource has its own cost. The coefficient R sets a regulated trade-off between the inaccuracy of the actuator governance to the external input $y_t$ and the total expenses for the control assistance $u_t$. We can track variation of R even visually. When R is small enough, the actuator follows the external signal relatively accurately, whereas the applied assisting control demonstrates high values and oscillating nature (see, e.g., Figure 3 (left)). The high unit cost R leads to small assisting control $u_t$ values and low quality concerning how the trolley coordinate follows the input $y_t$. Figure 3 (right) can serve as an illustration of this fact.

The provided calculations disclose the main issue of the control realization. The point is that the Euler–Maruyama scheme for the numerical solution to the SDE demonstrates an instability, which leads to divergence of some estimate trajectories. Regardless of the step $\delta$, almost each trajectory ${\widehat{y}}_t$ contains the points with violated conditions concerning the non-negativity ${\left({\widehat{y}}_t\right)}_1\ge 0,$ ${\left({\widehat{y}}_t\right)}_2\ge 0,$ ${\left({\widehat{y}}_t\right)}_3\ge 0$ or normalization ${\left({\widehat{y}}_t\right)}_1+{\left({\widehat{y}}_t\right)}_2+{\left({\widehat{y}}_t\right)}_3=1$. The number of such points varies from single to several thousands, which leads to the filtering “blowing up” and inability of the control synthesis. The more such defective estimates there are in the trajectory, the more likely the trajectory diverges rapidly.

If the defects appear rarely, they can be neutralized by a heuristic modification of the Euler–Maruyama scheme. Once the filtering estimate violates either non-negativity or normalization conditions ${\widehat{y}}_{t_i}$, it is replaced by the CTMC stationary distribution ${\pi }_{\infty }$, or by the estimation at the previous time instant ${\widehat{y}}_{t_{i-1}}$. These modifications help in some cases but are not a panacea, and this fact is shown in Section 5. Further in the presentation we refer to the introduced scheme modifications as ${\widehat{y}}_{t_i}^{lim}$ and ${\widehat{y}}_{t_i}^{del}$, respectively.

The next section presents a stable realization of the Wonham filter [43,47] required for the synthesis of the optimal stabilization strategy ${\widehat{u}}_t$.

4. Stable Filtering Algorithms by Discretized Observations

We realize both the filtering and stabilization algorithms with the same constant time increment $\delta$, such that $t_0=0,t_i=i\delta$, and $\frac{T}{\delta }$ is an integer.

Furthermore, without loss of generality, we suppose for model (1) $a_t\equiv \mathrm{const}=a_{t_i},b_t\equiv \mathrm{const}=b_{t_i}$ and ${\sigma }_t\equiv \mathrm{const}={\sigma }_{t_i}$ over the discretization intervals $\left(t_{i-1},t_i\right]$; otherwise, the functions $a_t,b_t$ and ${\sigma }_t$ should be approximated by suitable piece-wise constant functions.

Let us consider $z_t^0$, which represents a non-anticipated transformation of the observable input $z_t$:

$$z_t^0={\int }_0^t\left(d{z}_{\tau }-{(b}_{\tau }{z}_{\tau }+c_{\tau }{u}_{\tau })d\tau \right)={\int }_0^t\left(a_{\tau }{y}_{\tau }d\tau +{\sigma }_{\tau }d{w}_{\tau }\right).$$

The process $z_t^0$ does not depend on the control $u_t$, and due to the identity

$${\mathcal{F}}_t^z=\sigma \left\{z_{\tau },0\le \tau \le t\right\}\equiv \sigma \left\{z_{\tau }^0,0\le \tau \le t\right\}={\mathcal{F}}_t^{z^0}$$

the equality ${\widehat{y}}_t=\mathbb{E}\left\{y_t|{\mathcal{F}}_t^z\right\}=\mathbb{E}\left\{y_t|{\mathcal{F}}_t^{z^0}\right\}$ holds.

Since the numerical realization of the control strategy $u_t$ presumes its action at the time instants $t_i=i\delta$, we have to estimate the Markovial drift $y_t$ at the same moments. We use the transformed observations $z_t^0$ discretized by time with the step $\delta$:

$${\Delta z}_{t_i}^0={\int }_{t_{i-1}}^{t_i}\left(a_{\tau }{y}_{\tau }d\tau +{\sigma }_{\tau }d{w}_{\tau }\right).$$

The discretized observations generate the new family of $\sigma$ algebras

$${\mathcal{F}}_{t_i}^{\Delta z^0}=\sigma \left\{{\Delta z}_{t_j}^0,1\le j\le i\right\}.$$

If ${\mu }_i={\int }_{t_{i-1}}^{t_i}{y}_{\tau }d\tau ={\left({\mu }_i^1,\dots ,{\mu }_i^{n_y}\right)}^{{}^{\prime }}$ is a random vector composed of the occupation times of $y_t$ in each state during the interval $(t_{i-1},t_i]$, and $\mathcal{N}(z;m,{\sigma }^2)$ is a Gaussian probability density function with the mean m and variance ${\sigma }^2$, then the estimate ${\widehat{y}}_{t_i}=\mathbb{E}\left\{y_t|{\mathcal{F}}_{t_i}^{\Delta z^0}\right\}$ can be calculated by the following recursive procedure [43,47]

$${\widehat{y}}_{t_i}={\left(\mathbf{1}{\widehat{q}}_{t_i}^{{}^{\prime }}{\widehat{y}}_{t_{i-1}}\right)}^{-1}\left({\widehat{q}}_{t_i}^{{}^{\prime }}{\widehat{y}}_{t_{i-1}}\right),{\widehat{y}}_0={\pi }_0,$$

(10)

where $\mathbf{1}={(1,\dots ,1)}^{{}^{\prime }}\in {\mathbb{R}}^{n_y}$ is a vector formed by units, and the matrix ${\widehat{q}}_{t_i}={∥{\widehat{q}}_{t_i}^{k,j}∥}_{k,j=1}^{n_y}$ consists of the random values

$${\widehat{q}}_{t_i}^{k,j}=\mathbb{E}\left\{\mathcal{N}\left({\Delta z}_{t_i}^0;a_{t_i}{\mu }_i,\delta {\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)y_{t_i}^j|y_{t_{i-1}}=e_k\right\}.$$

(11)

The CMEs ${\widehat{q}}_{t_i}^{k,j}$ represent the integrals—the shift/scale mixtures of some Gaussians with the distribution of ${\mu }_i$ as the mixing one. The issue is that this distribution is not continuous with respect to the Lebesgue measure, so the integral (11) cannot be calculated analytically: we have to do it numerically. In this paper, we compare the following numerical schemes for the calculation of ${\widehat{q}}_{t_i}^{k,j}$:

$${\widehat{q}}_{t_i}^{k,j}\approx {\left({\breve{q}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)}_{kj}=\mathcal{N}\left({\Delta z}_{t_i}^0;\delta a_{t_i}{e}_k,{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)\left({\Delta }_{kj}+\delta {\lambda }_{kj}\right)$$

which is a scheme of the “left” rectangles with the accuracy order $\frac{1}{2}$,

$${\widehat{q}}_{t_i}^{k,j}\approx {\left({\breve{q}}_{t_i}^{\delta }\right)}_{kj}={\Delta }_{kj}{e}^{{\lambda }_{kk}\delta }\mathcal{N}\left({\Delta z}_{t_i}^0,\delta a_{t_i}{e}_k,{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)+$$

$$\left(1-{\Delta }_{kj}\right)\delta {\lambda }_{kj}{e}^{\frac{\delta }{2}\left({\lambda }_{kk}+{\lambda }_{jj}\right)}\mathcal{N}\left({\Delta z}_{t_i}^0,\frac{\delta }{2}{a}_{t_i}{(e}_k+e_j),{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)$$

which is a scheme of the “midpoint” rectangles with the accuracy order 1,

$${\widehat{q}}_{t_i}^{k,j}\approx {\left({\breve{q}}_{t_i}^{{\delta }^2}\right)}_{kj}={\Delta }_{kj}{e}^{{\lambda }_{kk}\delta }\mathcal{N}\left(\Delta z_{t_i}^0,\delta a_{t_i}{e}_k,{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)+$$

$${(1-\Delta }_{kj})\frac{\delta }{2}\left(e^{\left({\lambda }_{kk}-{\lambda }_{jj}\right)\frac{(\sqrt{3}-1)\delta }{2\sqrt{3}}}\mathcal{N}\left({\Delta z}_{t_i}^0,\frac{(\sqrt{3}-1)\delta }{2\sqrt{3}}{a}_{t_i}{e}_k+\frac{(\sqrt{3}+1)\delta }{2\sqrt{3}}{a}_{t_i}{e}_j,{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)\right.+$$

$$\left.e^{\left({\lambda }_{kk}-{\lambda }_{jj}\right)\frac{(\sqrt{3}+1)\delta }{2\sqrt{3}}}\mathcal{N}\left({\Delta z}_{t_i}^0,\frac{(\sqrt{3}+1)\delta }{2\sqrt{3}}{a}_{t_i}{e}_k+\frac{(\sqrt{3}-1)\delta }{2\sqrt{3}}{a}_{t_i}{e}_j,{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)\right)+$$

$$\sum _{i:i\ne j,i\ne k}^{}\frac{{\delta }^2}{6}\left(\left(e^{\frac{\delta }{6}\left({\lambda }_{kk}-{\lambda }_{ii}\right)+\frac{\delta }{6}\left({\lambda }_{ii}-{\lambda }_{jj}\right)}\mathcal{N}\left({\Delta z}_{t_i}^0,\frac{\delta }{6}{a}_{t_i}{(e}_k+e_i+{4e}_j),{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)\right.\right.+$$

$$e^{\frac{2\delta }{3}\left({\lambda }_{kk}-{\lambda }_{ii}\right)+\frac{\delta }{6}\left({\lambda }_{ii}-{\lambda }_{jj}\right)}\mathcal{N}\left(\Delta z_{t_i}^0,\frac{\delta }{6}{a}_{t_i}{(e}_k+{4e}_i+e_j),{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)+$$

$$\left.e^{\frac{\delta }{6}\left({\lambda }_{kk}-{\lambda }_{ii}\right)+\frac{2\delta }{3}\left({\lambda }_{ii}-{\lambda }_{jj}\right)}\mathcal{N}\left({\Delta z}_{t_i}^0,\frac{\delta }{6}{a}_{t_i}{(4e}_k+e_i+e_j),{\sigma }_{t_i}{\sigma }_{t_i}^{{}^{\prime }}\right)\right)$$

which is a Gaussian quadrature scheme with the accuracy order 2 (here and below ${\Delta }_{kj}$ stands for the Kronecker symbol).

In the system under investigation ${\left({\breve{q}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)}_{13}={\left({\breve{q}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)}_{31}={\left({\breve{q}}_{t_i}^{\delta }\right)}_{13}={\left({\breve{q}}_{t_i}^{\delta }\right)}_{31}={\left({\breve{q}}_{t_i}^{{\delta }^2}\right)}_{13}={\left({\breve{q}}_{t_i}^{{\delta }^2}\right)}_{31}=0$ since ${\lambda }_{13}={\lambda }_{31}=0$, and this fact simplifies calculations significantly.

For the comparative study we choose the approximations ${\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{y}}_{t_i}^{\delta }$, and ${\breve{y}}_{t_i}^{{\delta }^2}$, calculated by (10) and (11) and the schemes ${\breve{q}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{q}}_{t_i}^{\delta }$, and ${\breve{q}}_{t_i}^{{\delta }^2}$, respectively.

5. Comparative Numerical Study

Through a series of numerical experiments, we try to answer the question of how the accuracy of the filtering approximation affects the total performance of the control of a linear system (1) or, more precisely (7). The first set of tests is devoted to the control of the physical system introduced in Section 3. The investigation object of the second set represents a semi-stable linear system. The third set returns to a consideration of the trolley functioning under various complications. Finally, the fourth set is directed to the unstable linear system control.

We perform all experiments with the time increments $\delta =0.005,0.001,0.0001$. On the one hand, the step $\delta$ is inversely proportional to the computational costs of the input $y_t$ filtering. On the other hand, the decrease of $\delta$ improves the filtering accuracy and raises the control performance. Hence, the value of $\delta$ is a trade-off between the computational costs and reasonable losses of the control quality. However, in some applied problems, the step $\delta$ is constrained from below by the frequency of the actual measuring sensors.

We simulate 1000 random trajectories of the Wiener process $w_t$ and CTMC $y_t$, which are the same for all experiments. Then we synthesize the optimal stabilization strategy ${\widehat{u}}_t$ by use of various approximations of the Wonham filter estimates and various time steps $\delta$. Finally, the control performance index is a result of the Monte-Carlo sampling over the simulated trajectories.

We present the obtained results both in the tables and figures. For convenient analysis, the headers of all tables contain the system parameters, and we mark the important results in bold font. To characterize the quality of the estimates ${\widehat{y}}_{t_i}$, ${\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{y}}_{t_i}^{\delta }$, ${\breve{y}}_{t_i}^{{\delta }^2}$, we calculate the values

$$\widehat{D}\left({\tilde{y}}_{t_i}={\widehat{y}}_{t_i}\right),\widehat{D}\left({\tilde{y}}_{t_i}={\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right),\widehat{D}\left({\tilde{y}}_{t_i}={\breve{y}}_{t_i}^{\delta }\right),\widehat{D}\left({\tilde{y}}_{t_i}={\breve{y}}_{t_i}^{{\delta }^2}\right),$$

$$\widehat{D}\left({\tilde{y}}_{t_i}\right)=\widehat{\mathbb{E}}\left\{\frac{\delta }{T}\sum _{i=1}^{\frac{T}{\delta }}{\left(c{y}_{t_i}-c{\tilde{y}}_{t_i}\right)}^2\right\},$$

where $\widehat{\mathbb{E}}$ denotes averaging over the simulated trajectories.

In the case of the Euler–Maruyama scheme divergence for ${\widehat{y}}_{t_i}$, we use its stable versions ${\widehat{y}}_{t_i}^{lim}$ and ${\widehat{y}}_{t_i}^{del}$, suggested in Section 3.

To analyse the stabilization quality, we calculate the performance functionals

$$\widehat{J}\left({\widehat{u}}_{t_i}\right),\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right),\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right),\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right),$$

where ${\widehat{u}}_{t_i}$, ${\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{u}}_{t_i}^{\delta }$, ${\breve{u}}_{t_i}^{{\delta }^2}$ are optimal control strategies calculated with the use of the estimates ${\widehat{y}}_{t_i}$, ${\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{y}}_{t_i}^{\delta }$, ${\breve{y}}_{t_i}^{{\delta }^2},$ respectively, and

$$\widehat{J}\left(u_{t_i}\right)=\widehat{\mathbb{E}}\left\{\delta \sum _{i=1}^{\frac{T}{\delta }}\left({\left(C{y}_{t_i}-x_{t_i}\right)}^2+R{\left(u_{t_i}\right)}^2\right)\right\}.$$

5.1. Stable System

We perform the first numerical experiment with the trolley model, which has the parameters (8), to compare the modified Euler–Maruyama scheme and the stable ones (10) and (11). The calculation results are accumulated in

Table 1. The results of the first numerical experiment.

$a=-1,b=-0.5,T=10,h=10,g=0.01,R=0.01,$$(c_1,c_2,c_3)=(1,0,-1),\mathsf{\Lambda }=\left[\begin{array}{ccc}-0.5& 0.5& 0\\ 0.5& -1& 0.5\\ 0& 0.5& -0.5\end{array}\right]$
	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$				$\widehat{J}\left({\widehat{u}}_{t_i}\right)$
$\delta$	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{lim}$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{lim}$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{del}$				${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{del}$
	∞				∞
$0.005$	$\mathbf{0}.\mathbf{2269}$	$\mathbf{0}.\mathbf{0496}$	$0.0496$	$0.0496$	$\mathbf{3}.\mathbf{1876}$	$\mathbf{1}.\mathbf{5909}$	$1.5913$	$1.5912$
	$\mathbf{0}.\mathbf{3954}$				$\mathbf{4}.\mathbf{5370}$
	∞				∞
$0.001$	$\mathbf{0}.\mathbf{0539}$	$\mathbf{0}.\mathbf{0485}$	$0.0485$	$0.0485$	$\mathbf{1}.\mathbf{6035}$	$\mathbf{1}.\mathbf{5576}$	$1.5576$	$1.5575$
	$\mathbf{0}.\mathbf{0492}$				$\mathbf{1}.\mathbf{5631}$
$0.0001$	$\mathbf{0}.\mathbf{0493}$	$\mathbf{0}.\mathbf{0492}$	$0.0492$	$0.0492$	$\mathbf{1}.\mathbf{5432}$	$\mathbf{1}.\mathbf{5431}$	$1.5431$	$1.5431$

.

We do not supply the example with an extra illustration, because Figure 3 can serve this purpose. The results of the first experiment can lead to the following conclusions.

• The utilization of the suggested numerical filtering schemes, based on the time discretization of the observations, provides stability to the filtering procedure for all experimental conditions. The superiority is obvious for large steps $\delta$. The offered stable modifications of the Euler–Maruyama scheme demonstrate the expectable behavior. There are many “divergent” filtering trajectories; they need correction for large $\delta$, so both modifications provide a low filtering quality. When $\delta$ decreases, the number of the “divergent” trajectories also decreases, and both modifications provide the acceptable filtering quality. They lose to the suggested stable schemes. When $\delta$ decreases, the probability that the “divergent” filtering trajectory vanishes too and the quality of all considered numerical filtering schemes are identical.
• In general, one “divergent” trajectory is enough to make the value $\widehat{D}\left({\widehat{y}}_{t_i}\right)$ arbitrarily huge. The absence of such trajectories in the sample of 1000 trajectories for $\delta =0.0001$ does not guarantee their absence in the large samples. We detect the “divergent” trajectories in the samples of the size ${10}^7$. Note that the modified Euler–Maruyama schemes provide the acceptable filtering quality in the case of the small $\delta$ values.
• On the set of 1000 simulated trajectories, all suggested stable numerical schemes demonstrate very close performance.
• The absolute values of all $\widehat{D}$ and $\widehat{J}$ are small enough. The reason might be the small noise intensity $g=0.01$. The calculation error and moderate sample size seem to make the main contribution to the performance index comparing with the theoretical properties of all considered estimates.

The second numerical experiment differs from the first one only in that the parameter g value is increased by 10 times: $g=0.1$. The calculation results are collected in

Table 2. The results of the second numerical experiment.

$a=-1,b=-0.5,T=10,h=10,g=0.1,R=0.01,$$(c_1,c_2,c_3)=(1,0,-1),\mathsf{\Lambda }=\left[\begin{array}{ccc}-0.5& 0.5& 0\\ 0.5& -1& 0.5\\ 0& 0.5& -0.5\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	$0.1857$	$0.1861$	$0.1861$	$0.0496$	$2.7371$	$2.7409$	$2.7411$	$2.7411$
$0.001$	$0.1827$	$0.1828$	$0.1828$	$0.1828$	$2.6836$	$2.6845$	$2.6845$	$2.6845$
$0.0001$	$0.1860$	$0.1861$	$0.1861$	$0.1861$	$2.7126$	$2.7128$	$2.7128$	$2.7128$

.

The results of the second experiment can lead to the following conclusions.

• There are no divergent filtering trajectories ${\widehat{y}}_{t_i}$ among the simulated sample for any considered step $\delta$. This does not imply their absence in the samples of the greater size; however, the modifications ${\widehat{y}}_{t_i}^{lim}$ and ${\widehat{y}}_{t_i}^{del}$ give an opportunity to neutralize divergence without significant loss of the estimate accuracy.
• All estimates have moderate quality because of relatively high noise intensity g. We can conclude that under a low signal–noise ratio the performance of all considered estimates becomes similar for any considered step $\delta$.

The numerical experiments above indicate that the reasonability of the stable numerical filtering scheme usage is questionable. To highlight their utility, we have to involve the models with more delicate examples, which help to expose the value of the stable numerical schemes.

5.2. Semi-Stable System

We impair the model quality and consider a semi-stable dynamic system. Let us remind the reader that trolley model (7) keeps its stability for any $a,b>0$. The equalities $a=b=0$ determine the stability bound of the system, and we investigate this case. Obviously, the formula $C=\left(\frac{c_1}{a},\frac{c_2}{a},\frac{c_3}{a}\right)$ for the stable states is not valid in this situation. We also make the following changes to the parameters: $c=\left(-1.5,-0.5,0.5\right)$, $C=\left(-1,0,1\right)$, $g=0.01$. The calculation results are collected in

Table 3. The results of the third numerical experiment.

$a=0,b=0,T=10,h=10,g=0.01,R=0.01,$$(c_1,c_2,c_3)=(-1.5,-0.5,0.5),C=\left(-1,0,1\right),\mathsf{\Lambda }=\left[\begin{array}{ccc}-0.5& 0.5& 0\\ 0.5& -1& 0.5\\ 0& 0.5& -0.5\end{array}\right]$
	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$				$\widehat{J}\left({\widehat{u}}_{t_i}\right)$
$\delta$	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{lim}$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{lim}$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
	${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{del}$				${\widehat{y}}_{t_i}={\widehat{y}}_{t_i}^{del}$
	∞				∞
$0.005$	$\mathbf{0}.\mathbf{6776}$	$0.0557$	$0.0558$	$0.0557$	$\mathbf{4}.\mathbf{4154}$	$1.6670$	$1.6672$	$1.6671$
	$\mathbf{1}.\mathbf{1360}$				$\mathbf{5}.\mathbf{1079}$
	∞				∞
$0.001$	$\mathbf{0}.\mathbf{3500}$	$\mathbf{0}.\mathbf{0528}$	$0.0528$	$0.0528$	$\mathbf{1}.\mathbf{6035}$	$1.6135$	$1.6135$	$1.6135$
	$\mathbf{0}.\mathbf{0742}$				$\mathbf{1}.\mathbf{5631}$
$0.0001$	$\mathbf{0}.\mathbf{0527}$	$\mathbf{0}.\mathbf{0523}$	$0.0523$	$0.0523$	$\mathbf{1}.\mathbf{6168}$	$\mathbf{1}.\mathbf{6158}$	$1.6158$	$1.6158$

. The results of the third experiment can lead to the following conclusions.

• The proposed stabilization strategy demonstrates good performance even in the semi-stable case.
• The proposed stable numerical schemes of the state filtering demonstrate remarkable superiority compared with the Euler–Maruyama scheme and its modifications. This situation can be seen for all considered time steps $\delta$.
• There were five divergent filtering trajectories in the simulated sample, which did not affect the filtering performance dramatically in the case $\delta =0.0001$. Nevertheless, even the fact of their existence justifies the usage of the stable filtering scheme.
• The experiment demonstrates no significant difference in the performance of the stable numerical schemes ${\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, ${\breve{y}}_{t_i}^{\delta }$ and ${\breve{y}}_{t_i}^{{\delta }^2}$.

Figure 5 contains typical state trajectories for the semi-stable case, governed by the same control law ${\widehat{u}}_{t_i}$, but calculated by the various drift estimates ${\widehat{y}}_{t_i}^{lim}$, ${\widehat{y}}_{t_i}^{del},$ and control law ${\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}$.

Figure 6 contains the uncontrolled variant of the same trajectory.

The numerical experiments with the chosen parameters make it possibile to expose the behaviour of the offered heuristic modernizations of the Euler–Maruyama scheme ${\widehat{y}}_{t_i}^{lim}$ and ${\widehat{y}}_{t_i}^{del}$ in comparison with the stable scheme ${\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, when the original Euler–Maruyama scheme diverges. For the case $\delta =0.001$, Figure 7 makes it possibile to compare the estimates $c{\widehat{y}}_{t_i}^{lim}$, $c{\widehat{y}}_{t_i}^{del}$, $c{\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}$ with the true values $c{y}_{t_i}$ visually. Figure 8 contains corresponding state trajectories governed by the controls ${\widehat{u}}_{t_i}\left({\widehat{y}}_{t_i}^{lim}\right)$, ${\widehat{u}}_{t_i}\left({\widehat{y}}_{t_i}^{del}\right)$ and ${\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}$.

The numerical experiments with the semi-stable system allow us to conclude that the proposed stable filtering estimates are very effective for stabilization control synthesis. The differences between the proposed schemes are still minor in the semi-stable case. The subsection below discusses the situation when the choice of the stable numerical scheme makes some sense.

5.3. Stable System with High-Frequency Changing Drift

The numerical experiments of the last subsection demonstrate that the modified estimates ${\widehat{y}}_{t_i}^{lim}$ and ${\widehat{y}}_{t_i}^{del}$ do not “blow up” at any trajectory. Nevertheless, in some cases, they “lose” the CTMC $y_t$, as we can see in Figure 7. Hence, we hold for the comparative study the original Euler–Maruyama scheme, keeping in mind the possibility of its “blowing up” divergence.

The experiments above expose several realistic cases, for which the proposed stable numerical estimates demonstrate remarkable superiority for the stabilization strategy synthesis compared with the Euler–Maruyama scheme and its modifications. Nevertheless, we still do not present a situation when the suggested stable schemes deliver different estimation quality. To do this, we complicate the conditions of the trolley functioning. Namely, we increase the intensity of the CTMC transitions. We do this to confirm the hypothesis that the proposed schemes show different precision when the probability that the CTMC $y_t$ has more than one transition during $(t_{i-1},t_i]$ is significant. To meet this condition, we have to increase the elements of TIM $\mathsf{\Lambda }$. The frequent input signal transitions lead to more “active” control. If the control cost R is rather high, then the optimal control is minor and does not affect the system state. Hence, in this experiment, we set $R=0.00001$, as in Section 3 (see Figure 4).

Table 4. The results with transition intensities 10.

$a=-1,b=-0.5,T=10,h=10,g=0.1,R=0.00001,$$(c_1,c_2,c_3)=(-1.5,-0.5,0.5),\mathsf{\Lambda }=\left[\begin{array}{ccc}-10& 10& 0\\ 10& -20& 10\\ 0& 10& -10\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	∞	$\mathbf{0}.\mathbf{3744}$	$\mathbf{0}.\mathbf{3735}$	$\mathbf{0}.\mathbf{3733}$	∞	$\mathbf{7}.\mathbf{1539}$	$\mathbf{7}.\mathbf{1169}$	$\mathbf{7}.\mathbf{1146}$
$0.001$	∞	$\mathbf{0}.\mathbf{3430}$	$\mathbf{0}.\mathbf{3430}$	$\mathbf{0}.\mathbf{3430}$	∞	$\mathbf{6}.\mathbf{4908}$	$\mathbf{6}.\mathbf{4894}$	$\mathbf{6}.\mathbf{4894}$
$0.0001$	$0.3354$	$0.3357$	$0.3357$	$0.3357$	$6.3079$	$6.3095$	$6.3095$	$6.3095$

and

Table 5. The results with transition intensities 50.

$a=-1,b=-0.5,T=10,h=10,g=0.1,R=0.00001,$$(c_1,c_2,c_3)=(-1.5,-0.5,0.5),\mathsf{\Lambda }=\left[\begin{array}{ccc}-50& 50& 0\\ 50& -100& 50\\ 0& 50& -50\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	∞	$\mathbf{0}.\mathbf{7582}$	$\mathbf{0}.\mathbf{7470}$	$\mathbf{0}.\mathbf{7461}$	∞	$\mathbf{10}.\mathbf{7081}$	$\mathbf{10}.\mathbf{6590}$	$\mathbf{10}.\mathbf{6558}$
$0.001$	∞	$\mathbf{0}.\mathbf{6474}$	$\mathbf{0}.\mathbf{6470}$	$\mathbf{0}.\mathbf{6470}$	∞	$\mathbf{10}.\mathbf{3166}$	$\mathbf{10}.\mathbf{3149}$	$\mathbf{10}.\mathbf{3149}$
$0.0001$	$0.6260$	$0.6277$	$0.6277$	$0.6277$	$10.2395$	$10.2415$	$10.2415$	$10.2415$

contain the results of the numerical experiments, performed with the transition intensities 10 and 50, respectively.

The results of Table 4 correspond to expectations. For $\delta =0.005$ and $\delta =0.001$, the quality of estimates of discretized filters is built in accordance with their theoretical convergence rate. It should be noted that the distinction is quite small and can be seen only in the third or fourth significant digit. The case with $\delta =0.0001$ can already be considered marginal. Assuming that the difference may show an upward trend with an increase in the intensity of the jumps, we will continue to complicate the operating conditions of the actuator.

Figure 9 illustrates the results corresponding to the transition intensities 10. It contains the typical state trajectory governed by the strategy ${\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}$, indistinguishable from ${\breve{u}}_{t_i}^{\delta }$ and ${\breve{u}}_{t_i}^{{\delta }^2}$. We compare it with the trajectory governed by the “ideal” control $u_{t_i}^{\ast }$ synthesized with the complete information about the input $y_t$ [50]. The figure shows that the potential enhancement of the stabilization performance is possible only through the increasing of the filtering accuracy.

The plots corresponding to the system with the transition intensities 50 are given in Figure 10.

The calculation presented in Table 5 meets the expected results to the same extent as the previous one. For $\delta =0.005$, the distinction in the quality of the estimates of the discretized filters and the corresponding stabilization strategies can be seen in the second or third significant digit. For $\delta =0.001$, there is the same difference and the same hierarchy of algorithm quality, but already by the third or fourth significant digit. This means that the order of the transition intensities of the input $y_t$ for such sampling steps provides quite a lot of implementations of more than one jump in the interval $\delta$, which gives an advantage to higher-order filters. Having achieved this result, it should be noted that according to Figure 10, the effectiveness of the stabilization strategy itself, even in the case of complete information about the state of the input $y_t$, is extremely low. This is the result of the deterioration of the operating conditions of the actuator to completely unrealistic parameters. The slight superiority of the Euler–Maruyama scheme ${\widehat{y}}_{t_i}$ in the case $\delta =0.0001$ is illusory: in fact, the sample of the trajectories does not contain the divergent ones.

5.4. Influence of MJP Dimensionality on Control Performance

The dimensionality of the CTMC state space can also indirectly affect the filtering and stabilization performance. Let us repeat the experiments of the last subsection but increase the CTMC dimensionality: $n_y=4$. As in the subsection above, we consider values of the transition rate: 10 and 50. The numerical results of the experiments are given in

Table 6. The results with the CTMC dimensionality 4 and transition intensities 10.

$a=-1,b=-0.5,T=10,h=10,g=0.01,R=0.00001,$$(c_1,c_2,c_3,c_4)=(-1.5,-0.5,0.5,1.5),\mathsf{\Lambda }=\left[\begin{array}{cccc}-10& 10& 0& 0\\ 10& -20& 10& 0\\ 0& 10& -20& 10\\ 0& 0& 10& -10\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	∞	$\mathbf{0}.\mathbf{3349}$	$\mathbf{0}.\mathbf{3345}$	$\mathbf{0}.\mathbf{3345}$	∞	$\mathbf{6}.\mathbf{3259}$	$\mathbf{6}.\mathbf{2977}$	$\mathbf{6}.\mathbf{2952}$
$0.001$	∞	$0.3126$	$0.3125$	$0.3126$	∞	$5.7229$	$5.7219$	$5.7219$
$0.0001$	∞	$0.3063$	$0.3063$	$0.3063$	∞	$5.5687$	$5.5687$	$5.5687$

and

Table 7. The results with the CTMC dimensionality 4 and transition intensities 50.

$a=-1,b=-0.5,T=10,h=10,g=0.01,R=0.00001,$$(c_1,c_2,c_3,c_4)=(-1.5,-0.5,0.5,1.5),\mathsf{\Lambda }=\left[\begin{array}{cccc}-50& 50& 0& 0\\ 50& -100& 50& 0\\ 0& 50& -100& 50\\ 0& 0& 50& -50\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	∞	$0.7011$	$0.6973$	$0.7011$	∞	$11.9747$	$11.8980$	$11.9165$
$0.001$	∞	$\mathbf{0}.\mathbf{6431}$	$\mathbf{0}.\mathbf{6425}$	$\mathbf{0}.\mathbf{6429}$	∞	$\mathbf{11}.\mathbf{1556}$	$\mathbf{11}.\mathbf{1533}$	$\mathbf{11}.\mathbf{1544}$
$0.0001$	∞	$0.6258$	$0.6258$	$0.6258$	∞	$11.0219$	$11.0219$	$11.0219$

, respectively.

Figure 11 and Figure 12 contain the corresponding typical trajectories.

The obtained results emphasize the necessity of using the stable numerical filtering schemes: even for the minimal time step $\delta =0.0001$, the Euler–Maruyama approximation of the Wonham filter “blows up”. The differences between the various stable numerical filtering schemes are still minor.

5.5. Unstable System

To analyze the bundle “stable numerical filtering scheme–optimal stabilization strategy” exhaustively, we finally omit the requirement that system (7) be stable. This means that inequalities $b<0$ and $b^2+4a<0$ are not valid. The numerical experiments in this subsection have no physical sense; they just illustrate the applicability of the stabilization strategy calculated with the aid of the stable numerical filtering scheme.

Note that the linearity of (7) admits the control synthesis without requiring system stability. If $v_t$ is an optimal strategy for the stable system with $a<0$ and $b<0$, then $-v_t$ is optimal for the unstable system with “the mirror parameters” $\overline{a}=-a>0$ and $\overline{b}=-b>0$, and the minimal value of the criterion remains the same as for the stable system (see Table 1, Table 2 and Table 3). Hence, we have to consider the case when the parameters a and b have different signs. We perform a numerical experiment with the parameters $a=-1$, $b=5$, $(c_1,c_2,c_3)=(-2.5,-1.5,1)$; the rest of parameters coincide with ones in (8).

Table 8. The stabilization of the unstable system.

$a=-1,b=5,T=10,h=10,g=0.01,R=0.01,$$(c_1,c_2,c_3)=(-2.5,-1.5,1),\mathsf{\Lambda }=\left[\begin{array}{ccc}-0.5& 0.5& 0\\ 0.5& -1& 0.5\\ 0& 0.5& -0.5\end{array}\right]$
$\delta$	$\widehat{D}\left({\widehat{y}}_{t_i}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{\delta }\right)$	$\widehat{D}\left({\breve{y}}_{t_i}^{{\delta }^2}\right)$	$\widehat{J}\left({\widehat{u}}_{t_i}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^{\frac{1}{2}}}\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{\delta }\right)$	$\widehat{J}\left({\breve{u}}_{t_i}^{{\delta }^2}\right)$
$0.005$	∞	$0.0645$	$0.0646$	$0.0646$	∞	$4.6579$	$4.6602$	$4.6595$
$0.001$	∞	$0.0604$	$0.0604$	$0.0604$	∞	$4.4521$	$4.4525$	$4.4525$
$0.0001$	∞	$0.0594$	$0.0594$	$0.0594$	∞	$4.4753$	$4.4754$	$4.4754$

contains the filtering and stabilization results for various time-discretization steps.

Figure 13 contains the typical system trajectories both in the controlled and uncontrolled cases.

The calculations confirm the effectiveness of the proposed stabilization strategy in the case of an unstable system. When the system is uncontrolled, its state diverges rapidly: the absolute values of the state attain ${10}^{19}$–${10}^{20}$ at the termination moment $T=10$. The Euler–Maruyama scheme is also unstable in this system. All the proposed stable schemes of numerical filtering demonstrate the same performance.

6. Conclusions

The object of investigation in the paper is a stabilization of an overhead crane affected by the uncertain external Markov piece-wise constant disturbances. The goal is to analyze the potential application of the stable approximations of the Wonham filter [2,25] to the system state stabilization in the case of incomplete information. Their utilization should have neutralized the problem of the instability of the Wonham filter approximations delivered by the Euler–Maruyama numerical scheme. A large number of the fulfilled numerical experiments lead to the following conclusions.

• The rational choice of the stable numerical scheme preventing filter divergence is a natural trade-off between the accuracy requirements and the level of statistical uncertainty. In the case of small noise and high-frequency time discretization of the observations, the computational cost of the high-order schemas could exceed its accuracy gain.
• One should not underestimate the divergence of the unstable filtering schemas. There is a multiplicity of mathematical models for which any Euler–Maruyama filtering trajectory diverges.
• All stable schemas demonstrate consistency, i.e., converge to the theoretic Wonham filter as the time-discretization step infinitely vanishes.
• The complex numerical experiments suggest the schema of the “middle” rectangles of the order 1 as the preferable one.