# State Estimation for Nonlinear Discrete-Time Systems with Markov Jumps and Nonhomogeneous Transition Probabilities

Authors: Shunyi Zhao; Zhiguo Wang; Fei Liu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

State estimation problem is addressed for a class of nonlinear discrete-time systems with Markov parameters and nonhomogeneous transition probabilities (TPs). In this paper, the optimal estimation mechanism of transition probability matrix is proposed in the minimum mean square error sense to show some critical points. Based on this mechanism, the extended Kalman filters are employed as the subfilters to obtain the subestimates with corresponding models. A novel operator which fuses the prior knowledge and the posterior information embedded in observations is developed to modify the posterior mode probabilities. A meaningful example is presented to illustrate the effectiveness of our method.


---

## Body

## 1. Introduction

Markov jump systems (MJSs) are special systems which involve both time-evolving and event-driven mechanisms and have been used to model a wide variety of dynamic systems with unpredictable system structures. For MJSs, the state of the systems has two components: the finite valued state (denoted as mode) and the state in a finite dimensional Euclidean space. Motivated by a wide spectrum of applications, state estimation for MJSs has received considerable attention in the literatures [1–7] over the last three decades. It is well known that the optimal estimator for linear MJSs without the information about the underlying mode sequence is obtained by a bank of Kalman filters which requires exponentially increasing memory and computation with time. In order to limit the resource requirement, some suboptimal but executable algorithms have been proposed, including the generalized pseudo-Bayesian (GPB) [2], the interacting multiple model (IMM) [3], the improved IMM method [4, 5], and the particle filters (PF) [6, 7]. Moreover, methods based on guaranteed-cost approach which attempt to guarantee that the covariance of state is upper bounded by a certain constant value are also derived in quantity [8–14].It is worth noting that the all the Markov chains considered in the aforementioned estimators are assumed to be homogeneous. More specifically, the transition probabilities (TPs) are constant and are not allowed to be time-variant during the entire estimation process [15]. However, since the external practical environment (humidity, temperature, etc.) always exhibits complex and rich dynamics, it is impossible to guarantee the TPs of all likely modes to be time-invariant. In other words, nonhomogeneous Markov process is a common phenomenon, not the exception in practice [16, 17]. On the other hand, TPs are introduced to depict the uncertainties of transitions between possible system behavior patterns and the accuracy of TPs is critical to the estimation performance of MM-type estimators [18]. In this sense, it is wise and meaningful to consider the nonhomogeneous Markov chain in systems. By observing this, nonlinear stochastic hybrid systems with state-dependent TPs were considered in [19], where PF was embedded into IMM method to obtain the IMM-PF. With the same motivation, state-dependent TPs were modeled in terms of Gaussian-shaped functions in [20] and the Gaussian hypotheses were used to reduce the computational cost of the developed algorithms. However, the existing literatures seem to lack a complete approach of state estimation for MJSs with uncertain time-varying TPs, until the H-infinity filter for nonhomogeneous linear MJSs was proposed in [21]. In this paper, the character of nonhomogeneous TPs is reflected by a higher-level transition matrix and TPs vary within a finite set in discrete-time domain. Later, the convenience brought by the assumption that the transition probability matrix (TPM) randomly takes value in a finite set was demonstrated analytically in [22] by comparing the optimal Bayesian equation derived under two different conditions. Nevertheless, all these estimators require the information about the TPM. More specifically, all the candidate TPMs and high-level transition matrix are required to be known accurately, making the resulting methods intractable in practical applications.In this paper, a recursive estimator is derived for nonlinear MJSs with nonhomogeneous TPs. Different from [21, 22], the information about the dynamic of TPM including candidate TPMs and high-level transition matrix is not required to be known completely. This is a particular advantage of the framework raised here, and to the best of authors’ knowledge, no similar situation has been considered. The Bayesian equation that exploits the likelihoods of TPM to update the prior PDFs is developed. For the output calculation, we are required to find both the likelihoods and the dynamic of TPM which unfortunately turn out to be impossible. Therefore, instead of considering all the possible values of TPM in an admissible region, a constant diagonally dominant TPM which is rudely determined in advance is employed in the entire estimation process. In this sense, the problems involved in the estimation of TPM are solved naturally. The extended Kalman filter (EKF) [23, 24] is utilized as the subfilter under different models to obtain the subestimates and posterior mode probabilities. The exponentially increased computational problem is solved by the IMM-type approximation. Once the posterior mode probabilities are available, a novel operator is proposed to modify the posterior mode probabilities. Specifically, the probabilities of dominant modes are modified to approach one automatically according to their dominating degree. This operation is reasonable, as if the probability of one model is much larger than the others and we can almost make sure that it is the real one.The content of this paper is organized as follows. In Section2, some definitions and preliminaries are given. Section 3 proposes Bayesian estimation of TPM. The derivation of the suboptimal implementation is given in Section 4, where we propose a novel operator to modify the posterior mode probabilities after time update. The objective of Section 5 is to demonstrate the effectiveness of our method with a simulation example. Finally, Section 6 concludes this paper.


## 2. Problem Formulation and Preliminaries

We consider the discrete-time nonlinear JMSs described by the following equations:(1a)xk=f(xk-1,rk)+wk,(1b)yk=h(xk,rk)+vk,where f(·) and h(·) are nonlinear transition and measurement functions, respectively; xk∈Rn is the continuous state in Euclidean space with known initial distribution of state x0; yk∈Rp is the measurement; wk and vk are, respectively, uncorrelated white process and measurement noises with known statistics; the sequence {rk} is a continuous-time discrete-valued Markov chain and takes values in a finite state space M={1,2,…,m}. In this paper, the Markov chain governing the system dynamics is nonhomogeneous, which means that the time k exists as a guard condition of the probability of a jump from mode i at time k-1 to mode j at time k as
(2)πijk=Pr(rk=j∣rk-1=i,k),
where i,j=1,2,…,m. Then, the TPM can be given as Πk=[π1k,π2k,…,πmk]T, where notation [·]T denotes the transpose of matrix or vector, πik=[πi1k,πi2k,…,πimk]T, and i=1,2,…,m. Note that if Πk is known as a prior, the existing algorithms such as IMM-PF method [2, 5, 6] can be extended to the nonhomogeneous cases easily by using corresponding TPM at each time step. Unfortunately, TPs always change randomly and it is unrealistic to obtain the dynamic trajectories of real TPM in practical applications. Hence, it is assumed in the paper [21] that for arbitrary k, Πk is unknown but takes its value in a finite set Θ={Π1,Π2,…,Πl}, according a Markov chain χ(k) with a high-level transition matrix θ=[θpq], where p,q=1,2,…,l and ∑q=1lθpq=1. In this sense, the discrete distribution of Πp, where p=1,2,…,l, can be considered as a point-mass approximation of Πk over the simplex of valid TPMs. However, the requirement that both the candidate elements in Θ and the high-level transition matrix θ should be available in advance is difficult to be satisfied in applications. That is to say, the finite-valued hypothesis is also impractical. In fact, a filter designer often uses a diagonally dominant constant matrix π=[π1,π2,…,πm]T to replace the true values in reality, where πi=[πi1,πi2,…,πim]T and ∑j=1mπij=1, for any i∈M. Intuitively, if no further operator is implemented, using such an inaccurate matrix directly may lead to performance degradation, as it failed to capture the real possible system behavior patterns. Therefore, the objective of this paper is to derive a recursive, filtering algorithm that yields satisfactory estimates of states in the cases that only the approximate matrix π is available.


## 3. Bayesian Estimation

It is well known that the recursive Bayesian equation of the JMSs with accurate TPs is given as(3)p(xk∣Yk)=1ak∑j=1mp(yk∣xk,rkj)∑i=1mPr(rkj,rk-1i∣Yk-1)×∫p(xk∣rkj,xk-1)p(xk-1∣rk-1i,Yk-1)dxk-1,
where ak is a normalization constant, p(yk∣xk,rkj) and p(xk∣rkj,xk-1) denote the likelihoods and system dynamic under mode j, respectively, Pr(rkj,rk-1i∣Yk-1) are the joint probabilities of system modes at current and previous time steps, and p(xk-1∣rk-1i,Yk-1) denote the prior PDFs of state under mode i, where i=1,2,…,m. To see the derivation in detail, one can refer to one of our recent papers [22]. Considering (3), show that the effect of TPs in filtering process is mainly exhibited by the terms Pr(rkj,rk-1i∣Yk-1), which govern the dynamic of modes. Note that Pr(rkj,rk-1i∣Yk-1) can be decomposed in a more compact way as
(4)Pr(rkj,rk-1i∣Yk-1)=Pr(rkj∣rk-1i,k)Pr(rk-1i∣Yk-1)=πijkPr(rk-1i∣Yk-1),
where Pr(rk-1i∣Yk-1) are the mode probabilities at a previous time step. At this point, we can see that using incorrect TPs will lead to the deviation of mode probabilities and reduce the estimation accuracy further due to the following fact:
(5)πijkPr(rk-1i∣Yk-1)≠πijPr(rk-1i∣Yk-1).
Thus, it is highly desirable and really nice to obtain the true TPM Πk at each time step. Unfortunately, Πk is assumed to be unavailable, and a natural solution of this problem is to estimate the unknown TPM recursively. In the following, we thus focus our attention on the problem of finding the posterior PDF p(Πk∣Yk) of the nonhomogeneous TPM in the Bayesian framework. In view of Bayes’ ruler, we have
(6)p(Πk∣Yk)=p(yk∣Πk,Yk-1)p(Πk∣Yk-1)p(yk∣Yk-1),
where p(yk∣Yk-1) is considered as a normalization constant, p(yk∣Πk,Yk-1) is the likelihood, and p(Πk∣Yk-1) denotes the prior PDF of Πk. Using the law of the total probability, the prior PDF p(Πk∣Yk-1) can be equivalently written as
(7)p(Πk∣Yk-1)=∫p(Πk∣Πk-1)p(Πk-1∣Yk-1)dΠk-1,
where p(Πk∣Πk-1) and p(Πk-1∣Yk-1) denote the dynamic of Πk itself and posterior PDF at a previous time step. Using the total probability theorem, the first term in the numerator of the right-hand side of (6) becomes
(8)p(yk∣Πk,Yk-1)=∑j=1mp(yk∣rkj,Πk,Yk-1)Pr(rkj∣Πk,Yk-1)=∑j=1mp(yk∣rkj,Πk,Yk-1)×∑i=1mPr(rkj∣rk-1i,Πk)Pr(rk-1i∣Πk,Yk-1)=∑j=1mp(yk∣rkj,Πk,Yk-1)×∑i=1mπijkPr(rk-1i∣Πk,Yk-1).
Defining the likelihood vector Λk and mode probability vector υk-1 as Λk=[Λk1,Λk2,…,Λkm]T, υk-1=[υk-11,υk-12,…,υk-1m]T, where Λkj=p(yk∣rkj,Πk,Yk-1) and υk-1i=Pr(rk-1i∣Πk,Yk-1), we obtain
(9)p(yk∣Πk,Yk-1)=∑j=1mΛkj∑i=1mπijkυk-1i=υk-1TΠkΛk.
By the Chapman-Kolmogorov equation, the denominator of the right-hand of (6) can be rewritten as
(10)p(yk∣Yk-1)=∫p(yk∣Πk,Yk-1)p(Πk∣Yk-1)dΠk=∫υk-1TΠkΛkp(Πk∣Yk-1)dΠk.
Substituting (9) and (10) into (6) yields
(11)p(Πk∣Yk)=υk-1TΠkΛkp(Πk∣Yk-1)∫υk-1TΠkΛkp(Πk∣Yk-1)dΠk,
where p(Πk∣Yk-1) is given by (7). Note that formula (11) is the key of the entire estimation process, as it tells that the known prior distribution p(Πk∣Yk-1) is updated using the likelihoods of TPM to obtain the posterior conditional PDF p(Πk∣Yk). Unfortunately, calculating (11) analytically is almost impossible as it involves the following problems: (i) the dynamic of Πk in (7) is difficult to obtain; (ii) notice that to get the posterior conditional distribution, one has to calculate the likelihood of Πk. However, since Πk is unknown and is also required to be estimated at time k, evaluating the likelihood directly is impossible. Due to the above problems, the implementation of (11) is intractable and an alternative approach or some necessary approximation will be introduced in the following section.Remark 1.
If we assume that the TPM takes in a finite setΘ, the dynamic of TPM in (7) can be expressed by the high-level transition probabilities θpq, where p,q=1,2,…,l. Moreover, we are not required to consider all possible values in an admissible region and only consider the likelihood of the values in Θ that need to be evaluated, making the problem considerable simpler.


## 4. Suboptimal Estimator

Since the system under consideration is nonlinear, in this paper, EKFs are utilized as subfilters. It is well known that EKF is one of the most widely used algorithms in the framework of nonlinear filtering. Letxk-1∣k-1i, Δk-1∣k-1i denote the estimate of state, covariance matrix computed by the i EKF at time k-1, respectively. Let υk-1i denote the mode probability of mode i at time k-1. Then, the PDF of state at time k-1 conditioned on rk-1i is denoted as
(12)p(xk-1∣rk-1i,Yk-1)=N(xk-1;xk-1∣k-1i,Δk-1∣k-1i),
where N(·) denotes the normal distribution. As done in IMM method, the optimal Bayesian equation is decomposed into a series of basic transitions. In this sense, the PDFs p(xk-1∣rkj,Yk-1) can be equivalently given as follows:
(13)p(xk-1∣rkj,Yk-1)=∑i=1mp(xk-1∣rk-1i,Yk-1)Pr(rk-1i∣rkj,Yk-1).
Note that
(14)υ-k∣k-1ij≜Pr(rk-1i∣rkj,Yk-1)=Pr(rkj∣rk-1i)Pr(rk-1i∣Yk-1)Pr(rkj∣Yk-1)=πijυk-1i∑i=1mπijυk-1i
we have
(15)p(xk-1∣rkj,Yk-1)=∑i=1mυ-k∣k-1ijN(xk-1;xk-1∣k-1i,Δk-1∣k-1i).
Here, the IMM approximation is employed. More specifically, m-component Gaussian mixture is replaced by a single Gaussian density to solve the computational sensitive problem as
(16)p(xk-1rkj,Yk-1)≈N(xk-1;x-k-1∣k-1j,Δ-k-1∣k-1j),
where j=1,2,…,m,
(17)x-k-1∣k-1j=∑i=1mυ-k∣k-1ijxk-1∣k-1i,Δ-k-1∣k-1j=∑i=1mυ-k∣k-1ij[(xk-1∣k-1i-x-k-1∣k-1j)ΤΔk-1∣k-1i+(xk-1∣k-1i-x-k-1∣k-1j)(xk-1∣k-1i-x-k-1∣k-1j)Τ].
Once x-k-1∣k-1j and Δ-k-1∣k-1j are available, m EKFs are executed to obtain the posterior filtered states xk∣kj and covariance Δk∣kj as follows:
(18)xk∣k-1j=f(x-k-1|k-1j,rkj),Δk∣k-1j=FkjΔ-k-1∣k-1j(Fkj)T+Qk,xk∣kj=xk∣k-1j+Kkj(yk-h(xk∣k-1j,rkj)),Δk∣kj=Δk∣k-1j-KkjSkj(Kkj)T,
where Kkj is the filter gain given as Kkj=Δk∣k-1j(Hkj)T(Skj)-1,
(19)Fkj=∂f(xk,rkj)∂xk|xk=x-k-1∣k-1j,(20)Hkj=∂h(xk,rkj)∂xk|xk=x-k-1∣k-1j,(21)Skj=HkjΔk∣k-1j(Hkj)T+Rk,
where Qk and Rk are the covariance matrices of the process noise and measurement noise, respectively. At this point, the posterior mode probabilities can be updated as
(22)υkj=p(yk∣rkj,xk)∑i=1mπijυk-1i∑j=1mp(yk∣rkj,xk)∑i=1mπijυk-1i,
where j=1,2,…,m. Considering (14) and (22), show that only the approximate TPs πij are employed to compute the mode probabilities which are used as “weights” of subfilters in the filtering process. That is to say, the estimation performance will be unsatisfied, as the mode probabilities deviate from their real trajectories. On the other hand, it is known that the ideal case of estimation of MJSs is that the underlying mode sequence is known. In such case, only one subfilter is required to calculate the estimate at every time step and the multiple-model estimation algorithm will degenerate to be the Kalman filter with corresponding model. In this sense, the information of TPs is not required any more in filtering process. By observing these, a novel operator is developed to modify the dominating mode probabilities after time update to make them approach to one according to their dominating degree. Without loss of generality, we present ρ(υkj) at time k as
(23)υ~kj=ρ(υkj)=(υkj)χk(υk1)χk+⋯+(υkm)χk,
where j=1,2,…,m; the index parameter χk is given as follows. For the cases with m>2 and βk+βk′>2/m, we have
(24)χk=βk-1/m1/m-βk′
else χk=1, where βk and βk′ are the maximum mode probability and minimum mode probability at time k, respectively. Note that it is not difficult to verify that β′≠1/m in the cases βk+βk′>2/m. Therefore, the denominator of the right-hand side of (23) will not be zero and can be used safely in computer simulation. For the cases with m=2 and βk-εk>0.1, the index parameter can be evaluated by
(25)χk=10(βk-εk)
else χk=1. At this point, we can see that the index parameter χk which governs the modification degree in our method is adaptive and varies automatically according to the dominance of posterior mode probabilities. Specifically, a large parameter is evaluated when the dominance is evident, otherwise either 1 or another relatively small value is employed. It is also necessary to point out that the assumption that the real TPs are constants is not required, as the compensation parameter is updated using online information at each sampling time step. The proposed method works effectively as long as the dominance is obvious. This is an intuitively meaningful property and can be considered a particular advantage of the method introduced here. After modifying the posterior mode probabilities, the final estimate is calculated by merging the posterior estimates xk|kj and covariance Δk|kj using the modified mode probability υ~kj as
(26)xk∣k=∑j=1mυ~kjxk∣kj,(27)Δk∣k=∑j=1mυ~kj[Δk∣kj+(xk∣kj-xk∣k)(xk∣kj-xk∣k)Τ].
The algorithm based on the above illustrations is detailed below.Step 1.
Initializations: given x0∣0i, Δ0∣0i and υ0i, set k=1.Step 2.
Interacting: compute mixed states x-k-1∣k-1j and Δ-k-1∣k-1j according to (16).Step 3.
Extended Kalman filtering: in accordance with the steps of EKF, obtain the xk∣kj, Δk∣kj, and posterior mode probabilities υkj according to (18)–(22).Step 4.
Modification: modify the posterior mode probabilities υkj according to (23) and (24) to obtain υ~kj.Step 5.
Output: calculate the overall estimate as (26), set k=k+1, and go back to Step 2.


## 5. Numerical Example

To illustrate the performance of our method, a nonlinear dynamic system modeled which is a slightly modified version of the one given in [25] was carried out. Consider
(28)xk=0.5xk-1+axk-11+xk-12+wk-1,yk=xk220+vk,
where wk~N(0,1) and vk~N(0,1). The mode-dependent parameters a are aggregated into three possible modes shown as a1=30 and a2=5. The nonhomogeneous TPM Πk takes its value in a finite set Θ={Π1,Π2,Π3}, where Π1, Π2, and Π3 are presented as follows:
(29)Π1=[0.80.20.20.8],Π2=[0.70.30.30.7],Π3=[0.60.40.40.6].
The high-level transition matrix θ is equal to
(30)θ=[0.90.050.050.050.90.050.050.050.9].
It is necessary to point out that all the candidate TPMs Πi, where i=1,2,3 and high-level transition matrix θ are unknown in advance and only an approximate TPM π which is chosen by designer crudely is available as
(31)π=[0.50.50.50.5].
Two different algorithms were performed:(i)
the IMM-EKF withπ;(ii)
the proposed algorithms (denoted as OA) withπ.To assess the performance of algorithms, the average root mean square (RMS) error based on 5000 times Monte-Carlo simulations is employed. All the filters are initiated with the prior distributionp(x0)~N(0,2). Results from 5000 Monte Carlo runs are shown in Figures 1 and 2 and Table 1.Table 1
Average results of different algorithms.
Average RMS error
IMM
OA
Position
0.907
0.864
Average computational time (sec.)
0.068
0.069Figure 1
Average RMS error of state.Figure 2
Variation of the index parameter in a typical simulation.In Figure1, the average RMS error of state based on 5000 Monte Carlo runs of different algorithms is shown. Together with Table 1 and Figure 1, we can clearly see that OA yields improvements in estimation performance significantly over IMM-EKF as expected. This is the result of using more reasonable mode probabilities, which are mainly caused by the operator developed. That is, the information embedded in the measurements is exploited to counteract the negative effect brought by the wrong TPs. At each time step, the posterior mode probability of dominating mode is modified to approach one according to the dominating degree. The index parameter in a typical run is given in Figure 2. We see that the index parameter varies with the time instants and constantly adjusts and realigns itself according to the dominance of the posterior mode probabilities. In this sense, the operator can be considered as an instrument that can increase the probabilities of dominating modes automatically. Another important issue in practical application is the computation time. Note that in OA, the index parameter is required to be calculated at every time step. However, the computational time of OA is almost the same as IMM-EKF. This is mainly because the computation of OA is quite simple, no circulation or inversion is required, making the proposed method more competitive in online applications.


## 6. Conclusions

This paper proposed a novel state estimation method for a class of nonlinear JMSs with nonhomogeneous TPs. In our work, the optimal Bayesian estimation of TPM is proposed and EKF is employed to handle the nonlinear problem. Moreover, an operator is developed to modify the posterior mode probabilities to compensate the effect caused by wrong TPs used. Simulation implicates that our method yields better estimation performance than IMM-EKF and shows more robustness against the wrong TPs.



---

*Source: 293456-2013-12-24.xml*