State estimation for Markovian Jump Linear Systems with bounded disturbances

Automatica 49 (2013) 3292–3303

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

State estimation for Markovian Jump Linear Systems with bounded disturbances✩ Hao Wu a , Wei Wang a , Hao Ye a,1 , Zidong Wang a,b a

Department of Automation, Tsinghua University, Beijing, 100084, China

b

Department of Information Systems and Computing, Brunel University, Uxbridge, Middlesex, UB8 3PH, UK

article

info

Article history: Received 13 June 2012 Received in revised form 16 June 2013 Accepted 14 August 2013 Available online 17 September 2013 Keywords: Markovian Jump Linear Systems Polyhedral disturbances Set-membership State estimation Minkowski sum

abstract In this paper, we investigate the state estimation problem for a class of Markovian Jump Linear Systems (MJLSs) in the presence of bounded polyhedral disturbances. A set-membership estimation algorithm is first proposed to find the smallest consistent set of all possible states, which is shown to be expressed by a union of multiple polytopes. The posterior probabilities of the system jumping modes are then estimated by introducing the Lebesgue measure, based on which the optimal point estimate is further provided. Moreover, a similarity relationship for polytopes is defined and an approximate method is presented to calculate the Minkowski sum of polytopes, which can help reduce the computational complexity of the overall estimation algorithm. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction MJLSs comprise an important class of stochastic dynamic systems as they are often used to model problems with abrupt variations occurring in system structure and parameters (Boukas & Yang, 1995; Costa, Fragoso, & Marques, 2005). Such problems commonly exist in practice and they may be caused by different reasons such as component failures or repairs, sudden environmental disturbances and change of subsystem interconnections (Ho, 2011; Hwang, Balakrishnan, & Tomlin, 2006; Wu & Ye, 2012). Motivated by a wide variety of applications, state estimation for MJLSs has received lots of attention; see for instance Andieu, Davy, and Doucet (2003), Blom and Bar-Shalom (1988), Cinquemani and Micheli (2006), Costa (1994), Doucet, Logothetis, and Krishnamurthy (2000), Eeah and Hwang (2009) and Orguner and Demirekler (2008). The interacting multiple model (IMM) algorithm proposed in Blom and Bar-Shalom (1988) is one of the most representative

✩ This work was supported by the 973 Program of China under grant 2010CB731800, and the National Natural Science Foundation of China under grants 60974059, 60736026, 61021063, 61203068, 61290324. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Tongwen Chen under the direction of Editor Ian R. Petersen. E-mail addresses: [email protected] (H. Wu), [email protected] (W. Wang), [email protected] (H. Ye), [email protected] (Z. Wang). 1 Tel.: +86 10 6279 0497; fax: +86 10 6278 6911.

0005-1098/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.automatica.2013.08.030

methods in this area, with which the estimates of the posterior probabilities of all possible system modes can also be provided in addition to state estimates. It is worth noting that most of the currently available state estimation methods for MJLSs including those in the aforementioned results were developed based on a prior stochastic knowledge of the noises being assumed. This assumption cannot be satisfied in many situations. Alternatively, only the bounding information related to the magnitude or energy of the noises is available. In fact, the noises modeled being unknown but bounded commonly exist in many practical problems such as target tracking (Schweppe, 1968), networked control with quantization (Goodwin, Haimovich, Quevedo, & Welsh, 2004) and robotic control (Jaulin, 2009). However, the results on state estimation for MJLSs with unknown but bounded disturbances are still limited and the methods developed from a stochastic point of view may not be applicable to solve the problem. This is mainly because it will be difficult to define a reasonable optimal estimation directly with the co-existence of stochastic and deterministic signals. In Liu, Ho, and Niu (2010); Liu, Ho, and Sun (2008) and Souza, Trofino, and Barbosa (2006), by employing the bounded real lemma (BRL) given in Pan and Basar (1996), H∞ filters were designed for a class of continuous-time stochastically stable MJLSs in the presence of energy bounded disturbances. In this paper, we shall focus on the state estimation problem for MJLSs with magnitude bounded disturbances, which were regarded as a far more natural model than the energy bounded disturbances in practice; see Bertsekas and Rhodes (1971) and

H. Wu et al. / Automatica 49 (2013) 3292–3303

Milanese and Vicino (1991). To address this problem, an indirect strategy to obtain point estimates will be presented by exploiting optimal estimation techniques based on set-membership filtering. More concretely, set-membership filtering will be first investigated to find the smallest set of all states consistent with the model structure, measurements received and the constraints on disturbances. The results will then serve as a basis to develop optimal point estimators similar to Milanese and Vicino (1991). However, different from single-mode systems considered in Milanese and Vicino (1991), the optimal point estimation for MJLSs with set-membership uncertainty involves an additional step, that is, to estimate the posterior probabilities of all possible modes. The main contributions of the papers are summarized as follows. (i) Considering a class of MJLSs in the presence of bounded polyhedral disturbances modeled similarly to those in Goodwin et al. (2004) and Walter and Piet-lahanier (1989), we propose a recursive algorithm for set membership estimation. It is shown that the smallest consistent state set is expressed by a union of multiple polytopes and the true state is contained in this set. The set of state estimates corresponding to all possible modes may be distributed dispersedly in state space, if we utilize a polytope or ellipsoid to cover the set, large conservativeness will be introduced in the estimation results. Therefore, we utilize a new shape, i.e. a union of polytopes, instead of a single convex polytope or an ellipsoid to represent the outer-bound of the smallest consistent state set. (ii) Inspired by Milanese and Vicino (1991), the minmax criterion is adopted to derive the optimal point estimate based on the set-membership estimation results. The posterior probabilities of the modes need to be estimated firstly. In contrast to the IMM algorithm (Blom & Bar-Shalom, 1988), as no stochastic information about the disturbances is available, the Lebesgue measure is introduced to estimate the posterior probabilities of the modes based on the following concept: the set of state estimates obtained under the mode with larger measure will have larger probability. Moreover, as the resulting set of state estimates from set-membership filtering may be neither convex nor simply connected, the point estimate of the consistent set for each mode is determined as the measures weighted outer Chebyshev center (Boyd & Vandenberghe, 2004). (iii) The Minkowski sum of polytopes is involved in our proposed algorithm, which may largely increase the number of faces and vertices of the resulting polytopes. To tackle this problem, a similarity relationship for polytopes is defined and an approximate calculation method is proposed. With this method, the number of faces and vertices can be effectively decreased. Thus the computational complexity for the subsequent steps of the estimation algorithm can be greatly reduced. The rest of the paper is organized as follows. In Section 2, the problem is formulated. In Section 3, a recursive algorithm is firstly proposed to find the smallest consistent state set. Then the mode posterior probability estimation is presented, based on which the point estimation method is provided. In Section 4, an accurate calculation method for affine transformation of a polytope is given. Then an approximate calculation method is developed for Minkowski sum operation and the analysis on the relative approximation error is provided. The computational complexity of the overall state estimation algorithm is discussed in Section 5, followed by the simulation results presented in Section 6. Finally, we conclude the paper in Section 7. 2. Problem formulation and preliminaries We consider a class of MJLSs as follows, xk = Aθk xk−1 + Bθk wk zk = Cθk xk + vk .

(1)

3293

The initial state satisfies x0 ∈ γ , where γ is a known polytope. xk ∈ ℜn denotes the system state. θk ∈ {1, . . . , N } denotes the modes governed by a Markov chain with the initial distribution following {p1 (0), . . . , pN (0)}, and N is the number of modes. The transition probability matrix

ψ11  ψ21 Ψ =  .. . 

ψ12 ψ22 .. .

ψN1 ψN2 and ψij = Pr ( θk+1

 ψ1N ψ2N  N ×N ..  ∈R . · · · ψNN = j| θk = i) denotes the transition probability ··· ··· .. .

from mode i to mode j. Aθk , Bθk , Cθk are matrices with appropriate dimensions under the mode θk . zk is the system output. wk , vk denote the disturbances satisfying

wk ∈ α vk ∈ β where α, β are known polytopes and 0 ∈ α, β .

(2) (3)

Remark 1. Similar expressions to (2)–(3) for the magnitude bounded disturbances can also be found in many existing works in literature, such as the polyhedral sets in Kuntsevich and Lychak (1985) and Ong and Gilbert (2006), zonotopes in Alamo, Bravo, Redondo, and Camacho (2008), bounded 1-norm and ∞-norm models in Goodwin et al. (2004) and Vicino and Zappa (1996). Since the zonotopes, bounded 1-norm and ∞-norm models can all be regarded as special cases of the polyhedral model, a more general model for the disturbances is adopted in our paper. Remark 2. In Liu et al. (2010, 2008) and Souza et al. (2006), based on the bounded real lemma (BRL) for MJLSs presented in Pan and Basar (1996), H∞ filters were designed for a class of continuoustime stochastically stable MJLSs with energy bounded disturbances. In contrast to this, the energy bounded condition cannot be satisfied with the magnitude bounded disturbances modeled by (2) and (3). Moreover, the class of MJLSs (1) considered in this paper is not restricted to be stochastically stable. We now give the definition of ‘‘the smallest consistent set’’ as follows. Definition 1. For the initial state x0 ∈ γ , the smallest consistent set is defined as δf (0) = γ . We denote the smallest consistent set corresponding to the state xk at time k as δf (k), then xk+1 is regarded as a consistent state if there exist possible θk+1 , xk ∈ δf (k), wk+1 ∈ α and vk+1 ∈ β such that xk+1 = Aθk+1 xk + Bθk+1 wk+1 and zk+1 = Cθk+1 xk+1 +vk+1 hold. Based on this, all the consistent states form the smallest consistent set. The word ‘‘smallest’’ indicates that there is no inconsistent state existing in the set. The main objectives in this paper are summarized as follows. (i) Finding the smallest consistent set of all possible states that is consistent with system dynamics (1), disturbance constraints (2) and (3), and collected measurements. (ii) Providing the point estimate based on the consistent state set obtained in (i) and the posterior probability estimate of all possible modes. To achieve these objectives, the following assumption is imposed. Assumption 1. xk ∈ H, where the region H is a given polytope; the matrices Aj , Bj , Cj , j = 1, . . . , N are known. Aj , j = 1, . . . , N are invertible. Remark 3. It can be seen that an additional condition on xk is made, i.e. xk ∈ H. Such a condition was also required in Goodwin et al. (2004) to ensure that the consistent state set obtained at each time instant is bounded. Moreover, Aj for j = 1, . . . , N being invertible is necessary for the implementation of affine transformation, which will be adopted in the later state estimator design.

3294

H. Wu et al. / Automatica 49 (2013) 3292–3303

3. State estimator design In this section, a recursive algorithm is first proposed to find the smallest consistent set of all possible states, which will be shown to be expressed by a union of multiple polytopes and the true state is contained in this set. The posterior probabilities of modes are then estimated by introducing the Lebesgue measure. Based on these, the optimal point estimate is provided. Considering the increasing computational complexity with the recursive algorithm proceeds, a moving horizon strategy is also presented. 3.1. A recursive algorithm for state set estimation We denote the set of state estimates obtained at time k for the case of θk by Xθk . Then a recursive algorithm is proposed as follows.

Fig. 1. A block diagram of Algorithm 1.

• Initialization: Xθ0 =j = γ , for j = 1, . . . , N. From (4), (5) and (8), we have Algorithm 1. • Time update: Xθ−k =i,θk−1 =j = Xθ−k =i

N j =1

=∪



Yk,i = ∪Nj=1 x : x = Aθk =i xk−1 + Bθk =i wk , xk−1



x : x = Aθk =i xk−1 + Bθk =i wk ,

xk−1 ∈ Xθk−1 =j , wk ∈ α



Xθ−k =i,θk−1 =j .

(4)

Xθk =i = Xθ−k =i ∩ Xθ∼k =i .

(6) (7)

Remark 4. In (4), Xθ−k =i,θk−1 =j denotes the time prediction of the state set consistent with the case of a transition from mode θk−1 = j to mode θk = i. In (5), Xθ−k =i denotes the time prediction of the state set consistent with θk = i for all possible θk−1 . Note that for the initial update (i.e. k = 1), because Xθ0 =j = γ for all j ∈ {1, . . . , N }, the obtained Xθ−1 =i,θ0 =j in (5) for any given i are also the same for all j. Besides, Xθ∼k =i in (6) represents the set of the states consistent with the collected measurement and the case that θk = i. For a better understanding, a block diagram of Algorithm 1 is provided in Fig. 1. The following results can be drawn. Theorem 1. Consider the MJLSs in (1) with disturbances modeled as (2) and (3) under Assumption 1. By adopting Algorithm 1, the following conclusions hold. (1) xk ∈ ∪Ni=1 Xθk =i ; (2) ∪Ni=1 Xθk =i is the smallest consistent state set. Proof. We shall use the induction approach to prove Theorem 1. (1) From the initialization, it holds for k = 0. Assume that it holds for k − 1, k ≥ 2. Eqs. (4)–(6) indicate that there is some i, for i ∈ 1, . . . , N, such that xk ∈ Xθ−k =i and xk ∈ Xθ∼k =i . Therefore, we



have xk ∈ ∪Ni=1 Xθ−k =i ∩ Xθ∼k =i = ∪Ni=1 Xθk =i from (7). (2) From the initialization, it holds for k = 0. Assume that ∪Ni=1 Xθk−1 =i is the smallest consistent set of state xk−1 . For each mode i, for i ∈ 1, . . . , N, the set of all possible xk consistent with system dynamics (1) and disturbance constraint (2) can be represented as follows.



= x : x = x¯ k−1 + Bθk =i wk , 

x¯ k−1 ∈ ∪Nj=1 Aθk =i Xθk−1 =j , wk ∈ α

 = ∪Nj=1 x : x = x¯ k−1 + Bθk =i wk ,  x¯ k−1 ∈ Aθk =i Xθk−1 =j , wk ∈ α .

= Xθ−k =i .

(9)

For each mode i, for i ∈ 1, . . . , N, the set of all possible xk consistent with system dynamics (1) and disturbance constraint (3) can be expressed by the set Xθ∼k =i in (6). Therefore, from (7) and (9), the consistent state set for xk can be represented by

∪Ni=1 Yk,i ∩ Xθ∼k =i = ∪Ni=1 Xθ−k =i ∩ Xθ∼k =i = ∪Ni=1 Xθk =i . 

(10)

Observing (4)–(7) in Algorithm 1, the operations of polytopes including affine transformation, intersection and the Minkowski sum are involved. Before we proceed with more details about the geometric shapes corresponding to the state sets derived from Algorithm 1, some properties of the aforementioned operations are first introduced. Lemma 1 (Grunbaum (2003) and Ziegler (1995)). (1) The affine transformation of a convex polytope is a convex polytope. (2) The intersection of two convex polytopes X1 and X2 , i.e., X1 ∩ X2 , is also a convex polytope or null. (3) The Minkowski sum of two convex polytopes X1 and X2 , i.e., X1 + X2 , {y : y = y1 + y2 , y1 ∈ X1 , y2 ∈ X2 }, is also a convex polytope. We then have the following results. Theorem 2. Under Assumption 1, Xθ−k =i,θk−1 =j , Xθ−k =i , and Xθk =i , given by (4), (5) and (7), for k = 1, 2, . . . , i = 1, . . . , N, can be represented as a union of finite convex polytopes. Proof. Similar to the proof of Theorem 1, the induction method is adopted to prove this theorem. From (2) and (3), the regions of wk and vk are polytopes. According to Lemma 1, (4), and Assumption 1, Xθ−1 =i,θ0 =j , i, j = 1, 2, . . . , N are convex polytopes or null. From (5), there is Xθ−1 =i =

Xθ−1 =i,θ0 =j , i, j = 1, . . . , N, therefore, Xθ−1 =i , i = 1, . . . , N are convex polytopes or null. As stated in Goodwin et al. (2004), Xθ∼k =i is a convex polytope. Then according to Lemma 1 and (7), Xθ1 =i , i = 1, . . . , N are convex polytopes or null. We assume that Xθ−k =i,θk−1 =j , i, j = 1, . . . , N can be represented

Yk,i = x : x = Aθk =i xk−1 + Bθk =i wk ,



xk−1 ∈ ∪Nj=1 Xθk−1 =j , wk ∈ α



= ∪Nj=1 Xθ−k =i,θk−1 =j

(5)

• Measurement update:   Xθ∼k =i = x : Cθk =i x + vk = zk , vk ∈ β ∩ H . • Correction:



∈ Xθk−1 =j , wk ∈ α



as a union of finite convex polytopes at time k, i.e. Xθ−k =i,θk−1 =j ,

∪t =i,j,1k Xti,j,k , where Xti,j,k , i, j = 1, 2, . . . , N denote convex polytopes, Ti,j,k , i, j = 1, 2, . . . , N denote the number of the convex polytopes T

(8)

H. Wu et al. / Automatica 49 (2013) 3292–3303

in the union operation. From (5), we have Xθ−k =i = ∪Nj=1 ∪ and

  T ∪Nj=1 ∪t =i,j,1k Xti,j,k ∩ Xθ∼k =i   T = ∪Nj=1 ∪t =i,j,1k Xti,j,k ∩ Xθ∼k =i

Ti,j,k t =1

i,j,k

Xt

Xθk =i =

(11)

is followed from (7). For the time instant k + 1, according to (4) and (11), we have Xθ−k+1 =s,θk =u

= x : x = Aθk+1 =s xk + Bθk+1 =s wk+1 ,  xk ∈ Xθk =u , wk ∈ α  = x : x = Aθk+1 =s xk + Bθk+1 =s wk+1 ,    Tu,j,k u,j,k xk ∈ ∪Nj=1 ∪t =1 Xt ∩ Xθ∼k =u , wk ∈ α  T = ∪Nj=1 ∪t =u,1j,k x : x = Aθk+1 =s xk + Bθk+1 =s wk+1 ,  u,j,k xk ∈ X t ∩ Xθ∼k =u , wk ∈ α (12) 

where s, u = 1, 2, . . . , N and each set {x : x = Aθk+1 =s xk +

This is mainly because stochastic uncertainties are involved in the set-membership state estimation results. To be concrete, according to the optimal criterion in Milanese and Vicino (1991), the point estimate for the system under mode θk can be represented by xˆ k,θk = argxˆ k,θ min maxx∈Xθ ∥ˆxk,θk −x∥. Since θk is a stochastic vari-

ˆ k,i pˆ i (k). pˆ i (k) is square error (MMSE), i.e. xˆ k = E xˆ k,θk = i=1 x the estimate of Pr(θi |zk ), which is the posterior probability of the case that the mode θk = i. Then the main challenge is to find pˆ i (k) and xˆ k,i . 

Remark 6. In most of the existing results based on the setmembership estimation technique such as Alamo et al. (2008), Bertsekas and Rhodes (1971), Milanese and Vicino (1991) and Schweppe (1968), the outer-bound of the consistent state set is represented by a single ellipsoid or polytope. Although the storage memory required and the computational complexity of the overall estimation methods may be limited by using such traditional types of outer-bounds, they may not be suitable in this paper. This is because the smallest consistent state set obtained here can be expressed by a union of polytopes. Moreover, the polytopes in this set corresponding to all possible jumping modes may be distributed dispersedly in state space. Therefore, if a single ellipsoid or polytope is utilized to cover the whole set, large conservativeness will be introduced in the final estimation results. To avoid this, the union of multiple polytopes is used to express the outer-bound in this paper. 3.2. Mode probability and point estimation As discussed in Remark 5, we now investigate the problem of finding the optimal point estimate based on the smallest consistent state sets obtained by Algorithm 1. In Milanese and Vicino (1991), the optimal point estimate for set-membership consistent set is determined according to an optimal criterion in the sense of minimizing the maximum possible estimation error. However, it cannot be applied directly to solve the problem in this paper.



N

3.2.1. Mode probability estimation Suppose that at time k, the estimates of the posterior probabilities of the modes pˆ i (k − 1) for i = 1, . . . , N are all available for use. We then aim to derive the estimates pˆ i (k). Based on the Bayesian formula, there is Pr(θk = i|zk ) =

p(zk |θk = i)Pr(θk = i) N



Bθk+1 =s wk+1 , xk ∈ Xt

Remark 5. Note that different from other estimation methods including Kalman filtering, Algorithm 1 based on the setmembership technique only aims to generate the smallest sets of possible states which are consistent with the model structure, the constraints on the uncertainties and the received measurements. However, such a form of state estimation results may bring new challenges in designing novel controllers. This motivates us to further develop an optimal point estimator in the subsequent parts on the basis of the set-membership estimation results, which are sometimes known as ‘‘set-membership uncertainties’’ (Milanese & Vicino, 1991).

k

k

able, xˆ k,θk is also a stochastic vector with respect to θk . Therefore, as shown in Lewis, Xie, and Popa (2008), we use the expectation of xˆ k,θk to be the optimal estimate xˆ k in the sense of minimum mean

u,j,k

∩ Xθ∼k =u , wk ∈ α} is a convex polytope u,j,k according to Lemma 1. Because Xt ∩ Xθ∼k =u may be empty, Ti,j,k ≤ k−1 N . The number of convex polytopes in Xθ−k =i and Xθk =i is N × Ti,j,k and Ni,k × Ti,j,k , respectively, where Ni,k ≤ N. 

3295

(13)

p(zk |θk = j)Pr(θk = j)

j =1

where p(zk |θk = i) represents the likelihood function of zk . Pr(θk = i) is the a priori probability of θk = i, which can be estimated by using pˆ i (k − 1) and the mode transition probability matrix Ψ . For p(zk |θk = i), since no stochastic information about the states and the disturbances is available except their possible regions, all the possible states in the set Xθk−1 =j can be regarded to share the same possibility for driving xk−1 to zk under mode i. Therefore, if there are a larger number of elements in Xθk−1 =j which are consistent with the received measurement zk for the case θk = i, then zk is more likely to be generated under this mode, i.e. p(zk |θk = i) is larger. Thus we use the Lebesgue measures of the consistent state sets to represent p(zk |θk = i), which will be further adopted to update the posterior probability estimates of the modes. The time update result Xθ−k =i,θk−1 =j represents the consistent set driven by xk−1 . Besides, the correction result Xθk =i,θk−1 =j represents the part of Xθ−k =i,θk−1 =j that is consistent with zk . Therefore, the following quantity can be used to express the relative size of the consistent part of Xθk−1 =j under mode i. rji = Mθk =i,θk−1 =j /Mθ−k =i,θk−1 =j where

Mθ−k =i,θk−1 =j

(14)

and Mθk =i,θk−1 =j denote the Lebesgue measures

Xθ−k =i,θk−1 =j

of and Xθk =i,θk−1 =j , respectively. Based on the Chapman–Kolmogorov equation, the estimate of the posterior probability of θk = i for i = 1, . . . , N is derived as pˆ i (k) = c

N 

pˆ j (k − 1)ψji rji

(15)

j =1

where c is the normalizing constant which is selected to ensure N that i=1 pˆ i (k) = 1. Remark 7. If for some i, i ∈ {1, . . . , N }, Xθk =i is not null, whereas all the rest Xθk =s for s ̸= i, s ∈ {1, . . . , N } are null, the actual mode for time k can be identified determinately to be i. Remark 8. An IMM algorithm for MJLSs with Gaussian random noises was proposed in Blom and Bar-Shalom (1988), where the estimates of the posterior probabilities of all possible modes were provided in addition to the state estimates. In fact, for the probability estimation of the modes, our strategy is similar to that of the IMM algorithm. In Blom and Bar-Shalom (1988), the N likelihoods

3296

H. Wu et al. / Automatica 49 (2013) 3292–3303

p(zk |θk = i) are updated from the innovations of the N Kalman filters (the measurement prediction errors). Since no stochastic information about the disturbances is available in this paper, we utilize the measure of the consistent state set to update the likelihoods as shown in (14), (15) for MJLSs with set-membership uncertainties. 3.2.2. Point estimation We now consider generating xˆ k,i from the consistent set Xθk =i . Inspired by Milanese and Vicino (1991), the point estimate for setmembership uncertainties will be derived by solving a minmax optimal problem. However, since the obtained set of state estimates in this paper may be non-convex and not simply connected as discussed in Theorem 2, new challenges will be brought about. Tk,i

From Theorem 2, we set Xθk =i , ∪s=1 Xθsk =i , where each Xθsk =i for s = 1, . . . , Tk,i is a polytope. The measure of Xθsk =i is denoted as

Msi,k , and the related calculation method will be presented in the next subsection. Based on the minmax criterion in Milanese and Vicino (1991), we use the outer Chebyshev center of Xθsk =i as the optimal point estimate for polytope Xθsk =i , and it can be computed as csi,k = arg min max ∥x − xˆ ∥2 . x∈Xθs =i k

xˆ

(16)

Then we use the measure weighting method to determine the point estimate of xk for mode i, i = 1, . . . , N, which is expressed by Ti,k 

xˆ k,i =

csi,k Msi,k

 Ti,k

s=1

Msi,k

(17)

s=1

and the final point estimate xk is xˆ k =

N 

xˆ k,i pˆ i (k).

(18)

i=1

3.2.3. Measure calculation for polytope The measures of the polytope and the union of multiple polytopes involved in (14) and (17) are difficult and impractical to be calculated directly through an integration operation. To avoid the complexity in deriving high-dimensional integration, a Monte Carlo numerical method (Robert & Casella, 2004) is adopted. Note that only the case with the union of multiple polytopes is considered here. The method can be directly applied to the case of single polytope. Suppose that there is a union of polytopes X = ∪Tt=1 Xt , where T denotes the number of the polytopes in the union. Let I denote the smallest axis-aligned hypercube that contains X . I can be determined by the maximum and minimum values of X on each axis, which are generated by comparing the coordinates of the vertices of X . The measure of I, which is denoted as M, can be calculated by multiplying the measure of I on each axis. We sample I according to the uniform distribution for Ξ times and assume ¯ samples falling within X . Then the measure of X , that there are Ξ denoted as MX can be approximately calculated by MX ≈

Ξ¯ · M. Ξ

xk,i subspace, the error is ±1 with comparison to Ξi . If we further assume that Ξ1 = Ξ2 = · · · = Ξn , then the relative error of the calculation of the Lebesgue measure can be quantified by

   (Ξ1 ± 1)n   Em = 1 − Ξ1 n  in the case of a uniform sampling. 3.3. Moving horizon strategy for online implementation From (5), (7), and the proof of Theorem 2, the number of union operations in Xθ−k =i will increase exponentially as Algorithm 1 proceeds. Besides, the number of the (n − 1)-dimensional faces as well as the vertices in Xθ−k =i may also be increased due to the intersection operations. To treat these issues, a moving horizon strategy is presented as follows. If either of the following cases occurs at the kth time instant, (a) the number of the union operations in Xθ−k =i exceeds a threshold Tu , (b) the number of the vertices in Xθ−k =i exceeds a threshold Tv or the number of the (n − 1)-dimensional faces in Xθ−k =i exceeds a threshold Tf ,

then we place the current Xθ−k =i with the smallest rectangular box containing Xθ−k =i .

Remark 10. The moving horizon strategy has been proved to be an effective method to tackle the computational complexity issues in many control and estimation problems. For example in Goodwin et al. (2004), the approximation using a rectangular box needs to be performed after a fixed period of time, which is known as ‘horizon length’. Besides, a larger ‘horizon length’ leads to less conservative estimation results, but larger computational complexity. In this paper, the approximation needs be performed once the number of the union operations, the number of the vertices and the number of the (n − 1)-dimensional faces in Xθ−k =i exceed the respective prescribed thresholds. Similar to the ‘horizon length’ in Goodwin et al. (2004), determining the thresholds of the moving horizon strategy in this paper is also a trade-off between the estimation accuracy and scalability of the algorithm. Remark 11. Observing (4)–(7) of Algorithm 1, several basic operations of polytopes including affine transformation, intersection and the Minkowski sum are involved. It should be noted that the precise calculations of Minkowski sum operations will also lead to a greatly increased number of faces and vertices for the resulting polytopes. To deal with this issue, we shall propose an approximated calculation method for the Minkowski sum in Section 4. Since the approximation approach will inevitably lead to approximation errors which may affect the state estimation accuracy, an analysis of the relative approximation errors will also be provided. Besides, some discussions on the computational complexity of the overall state estimation algorithm will be presented in Section 5. 4. Calculation for convex polytopes operations

(19)

Remark 9. Clearly, the Monte Carlo method utilized for approximately calculating the measures of polytopes will inevitably bring approximation errors. However, such errors can be reduced by increasing the number of the samplings Ξ . As mentioned previously, we sample the hypercube I in accordance with the uniform distribution. Assume that the number of the uniform samplings in the xk,i subspace is Ξi , then Ξ = Ξ1 × Ξ2 × · · · × Ξn . Thus for the

As discussed in Remark 11, the operations of the polytopes including affine transformation, intersection, and Minkowski sum are involved in Algorithm 1. In Broman and Shensa (1990), Mo and Norton (1990) and Walter and Piet-lahanier (1989), some algorithms were proposed to precisely or approximately determine the intersection of two polytopes. In this section, an accurate calculation method for affine transformation will be provided. Furthermore, as the existing accurate algorithms for computing the Minkowski sum of two polytopes in Barki, Denis, and Dupont

H. Wu et al. / Automatica 49 (2013) 3292–3303

3297

(2009), Fogel and Halperin (2007), Fukuda (2004) and Fukuda and Weibel (2005) will lead to greatly increased number of faces and vertices, a Minkowski sum approximation scheme will be proposed to treat this problem. As stated in Grunbaum (2003), Yemelichev, Kovalev, and Kravtsov (1984) and Ziegler (1995), a convex polytope can be expressed by a group of vertices or a group of bounded intersections of closed half-spaces. Although these two kinds of representations have a one-to-one relationship, a lot of calculation may be needed to perform a transformation from one to the other. Thus both kinds of representations will be used in this paper for easier implementation of the algorithm. 4.1. Calculation for affine transformation A convex polytope X can normally be expressed in the following two forms. (1) X , (D, d) = {x : Dx ≤ d}, where D is a matrix, d is a vector. Each active inequality denotes an (n − 1)-dimensional face of X . R R (2) X , {x : x = i=1 γi Λi , i=1 γi = 1, γi ≥ 0, i = 1, . . . , R}, v where X = {Λ1 , Λ2 , . . . , ΛR } are the set of vertices of X . We have the following results. The related proofs are easy and omitted here. Lemma 2. Given a polytope X and Y , {y : y = Ax, x ∈ X }, where A is an invertible linear transformation, then the polytope Y can be represented as Y , DA−1 , d = y : DA−1 y ≤ d









(20)

and the set of vertices of Y can be represented as Y

. . . , AΛR }.

v

= {AΛ1 , AΛ2 ,

Lemma 3. Given that Y , {y : y = −x, x ∈ X }, the polytope Y can be represented by Y , (−D, d) = {y : −Dy ≤ d}

(21) v

and the set of vertices of Y can be represented as Y = {−Λ1 , −Λ2 ,

. . . , −ΛR }.

4.2. A Minkowski sum approximation approach Before we present an approximation method for Minkowski sum operation, a similarity relationship between two polytopes is first introduced. Given two polytopes X˜ 1 ⊂ ℜn and X˜ 2 ⊂ ℜn , the ˜ 1 ,1 , Λ ˜ 1,2 , sets of vertices of X˜ 1 and X˜ 2 are denoted by X˜ 1v = Λ

  ˜ 1,T1 and X˜ 2v = Λ ˜ 2,1 , Λ ˜ 2 ,2 , . . . , Λ ˜ 2,T2 , respectively. T1 and ...,Λ T2 are the numbers of the vertices of X˜ 1 and X˜ 2 , respectively. The 

Fig. 2. X1 + X2 and its similar approximation X1′ .

Property 4. If τ → ∞, then X˜ 2 tends to ℜn . Note that the Minkowski sum is involved in (4) in Algorithm 1. We use X1 and X2 to represent the state set and disturbance, resets of vertices are denoted byX1v = spectively. The corresponding   v Λ1,1 , Λ1,2 , . . . , Λ1,T1 and X2 = Λ2,1 , Λ2,2 , . . . , Λ2,T2 . The (n − 1)-dimensional faces of polytope Xs for s = 1, 2 are denoted as Xs = (Ds , ds ), with Ds =

T

DTs,1 , DTs,2 , . . . , DTs,Πs



T

and

ds = ds,1 , ds,2 , . . . , ds,Πs . Πs denotes the number of active inequalities in (Ds , ds ), i.e. the number of (n − 1)-dimensional faces of polytope Xs for s = 1, 2. We now utilize a convex polytope X1′ , which is similar to X1 , to approximate X1 + X2 . It will be seen that by adopting this method, the number of the (n − 1)-dimensional faces and the vertices of the approximated X1 + X2 can be kept the same as that of X1 . Fig. 2 is given to illustrate the basic idea. The convex polytope X1′ with vertices Λ′1 , Λ′2 and Λ′3 is similar to X1 . For a large enough τ , X1 + X2 ⊂ X1′ can always be satisfied. We then aim to find the polytope similar to X1 with the smallest measure to approximate X1 + X2 . As the disturbance sets in (2) and (3) are symmetrical around zero, the disturbance sets can be regarded as X2 . Set the translation vector as ϑ = 0. Then the objective is to find an interior point of X1 , denoted as O, and a constant τ (> 1), such that the measure of the convex polytope X1′ (similar to X1 ), denoted by MX ′ ,



is minimized under the condition that X1 + X2 ⊂ X1′ , i.e., min MX ′ , O,τ

1

s.t. X1 + X2 ⊂ X1′ .

1

(23)

Definition 2. The polytope X˜ 1 is similar to X˜ 2 if T1 = T2 = T and there is an interior point of X˜ 1 (i.e. O), a positive constant τ and a translation vector ϑ , such that for each i = 1, 2, . . . , T , there is

Since it is difficult to obtain the explicit expression of MX ′ , solving 1 the optimization problem (23) will not be easy. It can be seen that if the interior point O is very close to some face, then τ must be very large to ensure that the condition X1 + X2 ⊂ X1′ is satisfied. Then it will result in a large MX ′ . To avoid this, we select the inner 1 Chebyshev center (Boyd & Vandenberghe, 2004) of X1 as O, i.e.,

  ˜ 2 ,i + ϑ − O = τ Λ ˜ 1,si − O Λ

O = arg max dist O1 , ℜn \ X1

following definition is made.

(22)

where si = 1, 2, . . . , T , sa ̸= sb if a ̸= b, ∀a, b ∈ {1, 2, . . . , T }. Some obvious properties of the similarity relationship are then provided. Property 1. If X˜ 1 is similar to X˜ 2 , then X˜ 2 is similar to X˜ 1 . Property 2. If X˜ 1 is similar to X˜ 2 and τ > 1, then X˜ 1 ⊂ X˜ 2 + ϑ ; if τ < 1, then X˜ 1 ⊃ X˜ 2 + ϑ ; if τ = 1, then X˜ 1 = X˜ 2 + ϑ . Property 3. If X˜ 1 is similar to X˜ 2 , then each pair of hyperplanes determined by two groups of vertices of X˜ 1 and X˜ 1 satisfying (22), ˜ 2,i and Λ ˜ 1,si , are in the same direction. i.e. Λ





(24)

O1 ∈X1

where ℜn \ X1 denotes the complementary set of X1 in ℜn , dist (O1 , ℜn \ X1 ) denotes the Euclidean distance from O1 to ℜn \X1 , which can be calculated by dist O1 , ℜn \ X1 =





min

η=1,2,...,Π1

d1,η − D1,η O1

∥D1,η ∥2

.

(25)

Remark 12. It should be noted that finding the inner Chebyshev center as defined in (24) consists in solving a linear programming (LP) problem (Serpedin, Chen, & Rajan, 2011). However, the solution may not be unique. Consider the case if X1 is a special regular

3298

H. Wu et al. / Automatica 49 (2013) 3292–3303

box, the middle point chosen as the interior point is more likely to result in the smallest approximation error. Motivated by this, for the more general case that X1 has an irregular shape, we may choose the middle point of the smallest box that covers X1 as the initial value for solving the LP problem. Then the approximation error is more likely to be made as small as desired. The smallest τ to ensure that X1 + X2 ⊂ X1′ will be found by the following three steps. Step 1: For η = 1, 2, . . . , Π1 , carry on the following operations. (1) Select randomly a vertex of X1 that is located in the hyperplane {Λ : D1,η Λ = d1,η }, denoted as Λ1,tη , and define the point set X¯ ηv = Λ : Λ = Λ2,s + Λ1,tη , Λ2,s ∈ X2v . ¯ η ∈ X¯ ηv by solving (2) Find the point Λ





¯ η = arg max Λ

Fig. 3. Approximate representation of the real Minkowski sum and cutting off the redundant expansion.

D1,η Λs − d1,η

Λs ∈X¯ ηv

(26)

∥D1,η ∥2

and let

¯ η. d¯ η = D1,η Λ

(27)

(3) Solve the equation D1,η τ Λ1,tη − O + O = d¯ η .

 





(28)

Let τη denote the solution of (28).

Step 2: Determine η¯ by solving

η¯ = arg max τη .

(29)

η=1,2,...,Π1

Step 3: For η = 1, 2, . . . , Π1 , let d¯¯ η = D1,η τη¯ Λ1,tη − O + O











(30)



X1′ = D1 , d¯¯

(31)



where d¯¯ , d¯¯ 1 , . . . , d¯¯ η , . . . , d¯¯ Π1

T

. For t1 = 1, 2, . . . , T1 . Let

¯¯ = τ Λ − O + O Λ t1 η¯ 1 ,t1   ¯¯ = Λ ¯¯ , . . . , Λ ¯¯ , . . . , Λ ¯¯ Λ 1 t1 T1 .

(32) (33)

Remark 13. In Step 1, it aims to determine the expansion coefficient of each (n − 1)-dimensional faces of X1 , i.e. τη . In Step 2, the maximum τη¯ for η = 1, . . . , Π1 is selected as the unified expansion coefficient τ to construct a similar approximation X1′ . The (n − 1)-dimensional faces and the vertices of X1′ are computed in Step 3. Proposition 1. X1′ =





D1 , d¯¯

¯¯ is the smallest with vertices set Λ

convex polytope that is similar to X1 with the interior point O selected by (24) and it contains X1 + X2 . Proof. It can be easily seen that X1′ is similar to X1 with τ = τη¯ , interior point O, and translation vector ϑ = 0. Each vertex of X1 + X2    Π ¯ ¯ ̸∈ D1,η , d1,η , Λ ¯ ∈ X¯ ηv . Λ ¯ η , η = 1, 2, belongs to ∪η=1 1 Λ :Λ . . . , Π1 , must belong to the vertiex set of X1 + X2 . From (26)–(31) ¯ η¯ is located in one hyperplane of X1′ , i.e., D1,η¯ Λ ¯ η¯ = and Property 3, Λ d¯¯ η¯ . Combining that X1 + X2 ⊂ X1′ , we have X1′ = smallest convex polytope.





D1 , d¯¯ is the

relative error can be ensured to fall into [0, 1]. As we know, for a polytope with n dimensions, the Minkowski sum will result in an expansion of the polytope from the perspective of (n − i)dimensional faces for i = 1, . . . , n − 1 as well as the vertices, in which the expansions of the (n − 1)-dimensional faces account for the largest proportion in most cases. To simplify the error evaluation, only the expansions of the (n − 1)-dimensional faces are considered to approximately derive the relative error as defined above. The case of xk ∈ ℜ2 is discussed first, then the results are extended to the more general case of xk ∈ ℜn . Fig. 3 is provided as an example with two polytopes in 2-dimensional space. We denote the inner quadrilateral (Λ1 Λ2 Λ3 Λ4 ) as one of the original polytopes, i.e. X1 , whereas the other original polytope is denoted by X2 . Then the shaded region represents the real expansions of 1-dimensional faces (i.e. the edges) of X1 after real X1 + X2 is performed. The outer quadrilateral Λ′1 Λ′2 Λ′3 Λ′4 , denoted as X1′ , is the polytope similar to X1 which covers X1 + X2 with the smallest measure. Thus X1′ is the result we utilize to approximate real X1 + X2 . Besides, the measure of real X1 + X2 is dominated by the sum of the areas of X1 and the shaded region. The area of the blank outer region represents an approximate measure of the absolute error of our proposed Minkowski sum approximation method (Ma ). Before we proceed with the mathematical expression of the relative error, some notations need to be defined. Let the normal distance from the interior point of X1 (O in Fig. 3) to the ηth edge of X1 be δη for η = 1, . . . , Π1 , where Π1 is the number of edges of X1 . Let the perpendicular distance between the ηth edge of X1 and its parallel expanded edge in the real X1 + X2 be φη . From (26)–(28), we have τη = (δη + φη )/δη , where τη is the result obtained by solving (28) which describes the expansion coefficient of the ηth edge. From (29), there exists an edge of overlapping with an edge of real X1 + X2 (see edge Λ′3 –Λ′4 in Fig. 3). The expanding coefficient for this edge is τη¯ = (δη¯ +φη¯ )/δη¯ . Thus the absolute error between the expansions of real X1 + X2 and X1′ corresponding to this edge is zero. From (30) and (33), it can be seen that X1′ is formed by choosing τ = τη¯ as the unified expanding coefficients for the rest of the edges in X1 . Observing Fig. 3, the outer quadrilateral Λ′1 Λ′2 Λ′3 Λ′4 (i.e. X1′ ) can be divided into four parts, i.e. triangles OΛ′1 Λ′2 , OΛ′2 Λ′3 , OΛ′3 Λ′4 and OΛ′1 Λ′4 . According to the definition of relative error given at the beginning of this part, the relative error of the measure of the approximated Minkowski sum corresponding to each triangle can be computed as





δη + φη τ δη

2 =1−

 τ 2 η

.

4.3. Relative error in the measure of Minkowski sum approximation

Eη = 1 −

Here we use Mt and Ma to represent the measure of true Minkowski sum and that of absolute error due to our proposed approximated calculation algorithm, respectively. Then the relative error is formally defined as er = Ma /(Ma + Mt ). Therefore, the

The results presented above for the case xk ∈ ℜ2 can easily be extended to the case xk ∈ ℜn , as we can also divide the similar approximation of the Minkowski sum into Π1 super cubic cones, where Π1 is the number of (n − 1)-dimensional faces. Then the

τ

(34)

H. Wu et al. / Automatica 49 (2013) 3292–3303

error proportion of each part, denoted by Π1 , η = 1, . . . , Π1 , can be expressed by

 Eη = 1 −

δη + φη τ δη

n =1−

 τ n η

τ

.

(35)

From (34) and (35), if X1 is a regular polyhedron similar to X2 which represents the disturbance, all Eη will be equal to zero. If the expanding coefficient corresponding to one (n − 1)-dimensional face of X1 is significantly different from that of another when real X1 + X 2 is performed, Eη will be large. Thus we can adopt maxη Eη or η Eη to judge whether the proposed Minkowski sum approximation method is acceptable in the sense of approximation accuracy. More concretely, if maxη Eη ≤ E0 or η Eη ≤ E0′ holds, where E0 and E0′ are prescribed thresholds, it indicates that the proposed Minkowski sum approximation algorithmis acceptable. We then discuss the case that maxη Eη > E0 or η Eη > E0′ . In fact, although such an error is relatively large in this case, it can still be reduced by introducing an extra modification technique in our proposed approximation approach. As shown in Fig. 3, in order to reduce the error caused by the expansion of the edge Λ1 Λ2 , we can cut off a portion of redundant expansions by drawing a straight line which is parallel to and close to the corresponding edge of real X1 + X2 . If the remaining part of the polytope is treated as the final approximated Minkowski sum, the modified approximation method  can effectively reduce maxη Eη or η Eη . Moreover, the total number of the vertices of the remaining polytope is still smaller than that of the polytope corresponding to the real Minkowski sum. In summary, the Minkowski sum approximation approach in this paper can effectively reduce the numbers of the (n − 1)dimensional faces and the vertices of the resulting polytopes at the cost of acceptable errors. This feature can significantly reduce the computational complexity of the operations involved in the later steps of our state estimation algorithm, which will be seen through the analysis in the next section. 5. Computational complexity analysis In the proposed state estimation method, the main operations occupying relatively large computational complexity include Minkowski sum approximation, intersection operations, Monte Carlo methods for measure calculation of polytopes and determining the outer Chebyshev centers. Before we proceed with the detailed complexity analysis, the notation O(∗) is defined as the upper bound of the running time of a particular algorithm, where ∗ only refers to the term with the highest order. 5.1. About the Minkowski sum approximation

¯ represent the maximum numbers of the (n − 1)Let Π and Π dimensional faces and the vertices of the polytopes for Minkowski sum and intersection. We now discuss the computational complexity of the approximated calculation method developed in Section 4.2 for the Minkowski sum operations involved in (4) of Algorithm 1. For determining the interior point by calculating the inner Chebyshev center in (24) with x ∈ ℜn , it is a linear programming (LP) problem with n + 1 variables as stated in Serpedin et al. (2011). Several effective methods have been proposed to solve the LP problem including the simplex algorithm, the ellipsoid algorithm (Khachiyan, 1979) and the interior algorithm (Karmarkar, 1984). In our paper, the interior algorithm is adopted, of which the time complexity is O((n + Π )3.5 ) according to James (1988). Then it can be seen that the calculations in Steps 1–3 are trivial with comparison to the preceding operation of calculating the inner Chebyshev center. Thus the computational complexity of the Minkowski sum approximation can be regarded as O((n + Π )3.5 ).

3299

However, in Fukuda and Weibel (2005) where a method to precisely calculate the (n − 1)-dimensional faces of Minkowski sum is proposed, the total number of LP problems involved is pΠ . Defining the degree of a vertex as the number of its adjacent vertices, then p represents the sum of the largest degrees of the two polytopes. The corresponding computation complexity is then derived as O((n + Π )3.5 pΠ ). Obviously, our proposed approximation method can reduce the computational complexity in calculating Minkowski sum. 5.2. About the intersection operation In this paper, we adopt the method in Mo and Norton (1990) to calculate the intersection of two polytopes X1 and X2 , which appears in (6) and (7) of the Algorithm. It is executed by following the stepwise procedure below. (i) Select an (n − 1)-dimensional face of X2 , compute the intersection of X1 and the half-space supported by the selected face of X2 . (ii) Select a new (n − 1)-dimensional face of X2 , calculate the intersection of the polytope obtained in last step and the new half-space supported by the selected face of X2 . (iii) Repeat Step (ii) until all the (n − 1)-dimensional faces of X2 have been selected. In each step, four main procedures are (n − 1)-dimensional face updates, vertiex updates, vertex–vertex adjacency list updates and vertex–face list updates. Only some operations with simple judgments in the form of linear equations and inequalities are involved. Thus the computational complexity for the intersection operation of two polytopes can be expressed by ¯ 2 n). O(Π Π 5.3. About the Monte Carlo method for measure calculation of polytopes For the Monte Carlo method in Section 3.2.3 which is used to calculate the measure of a polytope, suppose that there is an average of Ξ1 samples per space dimension (as in Remark 9). We need to judge whether these samples belong to the polytope or not. Thus the time complexity of this method can be computed as O(Ξ1n Π ), where Π is the number of (n − 1)-dimensional faces of the polytope. 5.4. About the outer Chebyshev center calculation Determining the outer Chebyshev centers of polytopes is required in generating the optimal point estimation. In contrast to the inner Chebyshev center, the outer Chebyshev center of a polytope is unique. In this paper, the recursive algorithm proposed in Botkin and Turova-Botkina (1994) with finite steps is utilized to find the outer Chebyshev center of a polytope. Besides, according to Theorems 2 and 3 of Botkin and Turova-Botkina (1994), the total number of the recursive steps is the number of the vertices of ¯ ). In each recursive step, a non-trivial the polytope (denoted as Π operation is involved to find the distance between a point and the convex hull of the point set. Such an operation is actually to solve a convex quadratic programming problem with time complexity O(n3 ) (Goldfarb & Liu, 1991). Thus the maximum time complexity ¯ ). of finding the outer Chebyshev center of a polytope is O(n3 Π It is also worth mentioning that an alternative approximate calculation method is proposed in Milanese and Vicino (1991) to find the outer Chebyshev center, in which only 2n trivial optimization problems need to be solved. Thus considerable computational complexity can be saved with this approximated calculation method, though calculation accuracy will be affected to some extent. Remark 14. In fact, since Minkowski sum approximation approach proposed in Section 4.2 can effectively decrease the numbers of the (n − 1)-dimensional faces and vertices of the resulted

3300

H. Wu et al. / Automatica 49 (2013) 3292–3303

polytopes, the computational complexity of the operations involved in the later steps of our state estimation algorithm can be reduced significantly. Take (7) in Algorithm 1 as an example, the intersection of Xθ−k =i and Xθ∼k =i needs to be calculated. From the aforementioned results, it is known that for fixed state dimensions (n), the computational complexity of such an intersection operation is determined by the numbers of the (n − 1)-dimensional faces and the vertices of Xθ−k =i and Xθ∼k =i . By adopting the Minkowski sum

approximation method, Xθ−k =i,θk−1 =j in (4) and Xθ−k =i in (5) with reduced numbers of faces and vertices can be obtained. Therefore, less memory is required to calculate the intersection in (7) with comparison to the case when the Minkowski sum was explicitly calculated in (4).

Fig. 4. The projection of Xθ7 =1 on x7,1 –x7,2 2D subspace.

corresponding to the two modes are A1 =

0.5 5.5. Computational complexity of the overall state estimation algorithm We proceed with the computational complexity of the overall state estimation algorithm in this paper. The overall algorithm can be roughly divided into three main parts, i.e. state set estimation, mode probability estimation and point estimation, which correspond to Algorithm 1 in Sections 3.1, 3.2.1 and 3.2.2, respectively. The following notations are first defined. N is the number of possible jumping modes, n is the dimension of the state xk . For time instant k, we use U to denote the maximum number of the poly¯ represent topes in the state sets Xθk −1=i for i = 1, . . . , N. Π and Π the maximum numbers of the (n − 1)-dimensional faces and the vertices of the polytopes obtained in (4)–(7). Ξ1 is the number of the Monte Carlo samples per space dimension.

• Part I. State set estimation From (4), (6) and (7), N 2 U Minkowski sum operations and N + NU intersection operations are involved in this part. The worstcase computational complexity for this part can thus be derived as

¯ 2 n). CC1 = O(N 2 U (n + Π )3.5 ) + O(NU Π Π

(36)

• Part II. Mode probability estimation From (14), 2N 2 U measure calculations of the polytopes by the Monte Carlo method are involved in mode probability estimation. Thus the worst-case computational complexity for this part can be derived as CC2 = O(N 2 U Ξ1n Π ).

(37)

• Part III. Point estimation For this part, there are N 2 U outer Chebyshev centers (see (16)) and the measures of N 2 U polytopes requiring to be calculated. Thus the corresponding computational complexity can be expressed by

¯ ) + O(N 2 U Ξ1n Π ). CC3 = O(N 2 Un3 Π

(38)

In summary, the computational complexity CC for the overall state estimation algorithm in this paper can be derived as CC = CC1 + CC2 + CC3

¯ 2 n + N 2 U Ξ1n Π = O(N 2 U (n + Π )3.5 + NU Π Π ¯ ). + N 2 Un3 Π

(39)

Remark 15. By adopting the moving horizon strategy presented in ¯ correspond to the preSection 3.3, the upper bounds of U, Π and Π scribed thresholds Tu , Tf and Tv , respectively. Thus from (39), for fixed n, N and Ξ , the upper bound of the worst-case computational complexity of the overall state estimation algorithm is determined by the values of Tu , Tf and Tv . 6. Simulation In simulation, we consider an MJLS as described in (1) with two possible jumping modes. The states xk ∈ ℜ3 , the matrices

0 0.7

0.5 0.8 0

0.7 0.7 0

0.7 0 −0.6

0 −0.7 0.6

 , A2 =

 0

0.8 0 .7

, respectively. B1 , B2 , C1 and C2 are identity matrices.

The Markov transition probability matrix is Ψ =

0.7 0.5



0.3 0.5



. The

polyhedral sets for the disturbances are α = {wk : ||wk ||υ ≤ 0.1}, β = {vk : ||vk ||υ ≤ 0.5} with υ = ∞. The initial state x0 is chosen to be x0,1 ∈ [1, 3], x0,2 ∈ [1, 3] , x0,3 ∈ [1, 3]. The initial distribution of the modes is {p1 (0), p2 (0)} = {1, 0}. The number of the samples for each dimension of xk is selected to be 20 in the Monte Carlo algorithm. We use a computer with the following configurations for simulation: Windows XP SP3, Intel (R) Core (TM) 2 Duo CPU, frequency 2.20 GHz, memory 1.96 GB. Matlab R2009a is adopted as the simulation software. The average CPU running time for the overall 10 time-steps is 0.14 s. In the following, we first present the detailed results obtained at each step of the state estimation method (i.e. consistent state set estimation, mode probability estimation and point estimations) for some k. As discussed in Remark 7, there exist special cases for some k such that the consistent state sets obtained corresponding to certain mode may be null. Then the implementation of the state estimation can be simplified at these cases as the mode can be identified straightforwardly. To better illustrate the implementation procedure of our proposed algorithm, the detailed results for k = 8 are provided. Note only the projections of the polytopes on x8,1 –x8,2 2D subspace are depicted for easier reading. At k = 8, the projection of set Xθ7 =1 is a polytope and shown in Fig. 4. But the set Xθ7 =2 is null, thus the mode for time k = 7 should be one, i.e. θ7 = 1. By the time update, the projections of the polyhedral sets Xθ−8 =1,θ7 =1 and Xθ−8 =2,θ7 =1 after the approximated Minkowski sum operation are shown in Fig. 5(a) and (b), respectively. By the measurement update, the projections of Xθ8 =1,θ7 =1 and Xθ8 =2,θ7 =1 after the intersection operation are provided in Fig. 6(a) and (b), respectively. We then derive the mode probability estimation and the point estimation based on the set-membership estimation results obtained as above. Note that the real mode for k = 8 is one. From (14) and (15), we have pˆ 1 (8) = 0.6426 and pˆ 2 (8) = 0.3574. The posterior probability estimate for mode 1 is larger than that for mode 2. Clearly, this result is consistent with the real situation. The outer Chebyshev centers of Xθ8 =1,θ7 =1 and Xθ8 =2,θ7 =1 are determined as [0.2532, 0.2219, 0.6949]T and [0.5279, −0.3020, 0.5760]T . Then from (18), the final point estimate for x8 is calculated as [0.2532, 0.2219, 0.6949]T × 0.6426 + [0.5279, −0.3020, 0.5760]T × 0.3574 = [0.3514, 0.0347, 0.6524]T . The simulation involves 24 Minkowski sum operations. The relative errors for the measure of the approximated calculation are presented in Fig. 7. Observing the 11th-time Minkowski sum operation, the relative approximation error is 15.85%. Note the numbers of vertices of X1 and real X1 + X2 are 22 and 88, respectively. However, the number of vertices of the approximated X1 + X2 is limited to 22, i.e. the same as that of X1 . Therefore, three-quarters of the vertices of Minkowski sum have been reduced at the cost of bringing approximately only 15.85% relative error. We then observe the

H. Wu et al. / Automatica 49 (2013) 3292–3303

(a) Xθ−8 =1,θ7 =1 .

3301

(b) Xθ−8 =2,θ7 =1 .

Fig. 5. The projection of time update results in x8,1 –x8,2 2D subspace.

(a) Xθ8 =1,θ7 =1 .

(b) Xθ8 =2,θ7 =1 .

Fig. 6. The projection of measurement update results on x8,1 –x8,2 2D subspace.

Fig. 8. The estimation results of x1 . Fig. 7. The relative errors of the approximate Minkowski sum.

14th Minkowski sum operation. The number of vertices of X1 and real X1 + X2 are 12 and 48, respectively. By adopting our proposed approximation algorithm, the relative approximation error reaches 44.60%. Therefore, modifications are required to reduce the error. After real X1 + X2 is performed, the expanding coefficients (τη ) with respect to the eight faces of X1 are 2.95, 2.95, 2.94, 1.68, 2.33, 1.24, 2.81, 2.49, respectively. Clearly, the largest τη is 2.95, which is adopted as τ to form the approximated polytope similar to X1 + X2 . Therefore, the largest difference between τη and τ lies in the 6th face. Hence we can utilize a face parallel to the 6th face to cut off the redundant expansion. By doing this, the relative error can be reduced from 44.60% to 14.79% without adding new vertices. The point estimates for xk,1 , xk,2 and xk,3 are shown in Figs. 8–10, respectively. The upper and lower bounds of the estimates are also provided, which correspond to the maximum and minimum values of the set-membership consistent state sets for each dimension. It can be observed that the actual state always locates between the upper bounds and lower bounds of the estimates. It indicates that

we have 100% confidence that the true states belong to the estimated region at any time step. Moreover, the actual states can be tracked effectively with the proposed point estimation algorithm. 7. Conclusion In this paper, state estimation is investigated for a class of Markovian Jump Linear Systems (MJLSs) in the presence of bounded polyhedral disturbances. An integrated optimal criterion is developed for the point estimation. In order to implement the optimal point estimate, a set-membership estimation algorithm is first proposed to find the smallest consistent set of all possible states, which is shown to be expressed by a union of multiple polytopes. Then the posterior probabilities of the system jumping modes are estimated by introducing Lebesgue measure, based on which the optimal point estimate is obtained from the probability weighted Chebyshev centers of the consistent polytopes under all jumping modes. Moreover, a similarity relationship for polytopes

3302

H. Wu et al. / Automatica 49 (2013) 3292–3303

Fig. 9. The estimation results of x2 .

Fig. 10. The estimation results of x3 .

is defined and an approximate method is presented to calculate the Minkowski sum of polytopes, which can help reduce the memory requirement and computational complexity. An analysis of the relative errors of the Minkowski sum approximation and the computational complexity of the overall state estimation algorithm are also provided. It is worth mentioning that the approximate Minkowski sum calculation algorithm in this paper is fixed regardless of the polytope shapes. It would be interesting to extend the algorithm such that it can be adapted to different types of polytopes. Besides, it is known that quantization commonly exists in networked control systems and the quantization errors are often regarded as unknown but bounded disturbances. Therefore, the state estimation problem for MJLS with quantized measurements may be considered as a future work. References Alamo, T., Bravo, J. M., Redondo, M. J., & Camacho, E. F. (2008). A set-membership state estimation algorithm based on DC programming. Automatica, 44(1), 216–224. Andieu, C., Davy, M., & Doucet, A. (2003). Efficient particle filtering for jump Markov systems, application to time-varying autoregressions. IEEE Transations on Signal Processing, 51(7), 1762–1770. Barki, H., Denis, F., & Dupont, F. (2009). Contributing vertices-based Minkowski sum computation of convex polyhedral. Computer-Aided Design, 41(7), 525–538. Bertsekas, D. P., & Rhodes, I. B. (1971). Recursive state estimation for a setmembership description of uncertainty. IEEE Transactions on Automatic Control, 16(2), 117–128.

Blom, H. A. P., & Bar-Shalom, Y. (1988). The interacting model algorithm for systems with Markovian switching coefficients. IEEE Transactions on Automatic Control, 33(8), 780–783. Botkin, N., & Turova-Botkina, V. (1994). An algorithm for finding the Chebyshev center of a convex polyhedron. In System modelling and optimization. Springer. Boukas, E. K., & Yang, H. (1995). Stability of discrete-time linear systems with Markovian jumping parameters. Mathematics of Control, Signals and Systems, 8(4), 390–402. Boyd, S., & Vandenberghe, L. (2004). Convex optimization. New York: SpringerVerlag. Broman, V., & Shensa, M. J. (1990). A compact algorithm for the intersection and approximation of N-dimensional polytopes. Mathematics and Computers in Simulation, 32(5–6), 469–480. Cinquemani, E., & Micheli, M. (2006). State estimation in stochastic hybrid systems with sparse observations. IEEE Transations on Automatic Control, 51(8), 1337–1342. Costa, O. L. V. (1994). Linear minimum mean square error estimation for discretetime Markovian jump linear systems. IEEE Transactions on Automatic Control, 39(8), 1685–1689. Costa, O. L. V., Fragoso, M. D., & Marques, R. P. (2005). Discrete-time Markov jump linear systems. Springer. Doucet, A., Logothetis, A., & Krishnamurthy, V. (2000). Stochastic sampling algorithms for state estimation of jump Markov linear systems. IEEE Transactions on Automatic Control, 55(1), 188–202. Eeah, C. E., & Hwang, I. (2009). Stochastic linear hybrid systems: modeling, estimation, and application in air traffic control. IEEE Transactions on Control Systems Technology, 17(3), 563–575. Fogel, E., & Halperin, D. (2007). Exact and efficient construction of Minkowski sums of convex polyhedra with applications. Computer-Aided Design, 39(11), 929–940. Fukuda, K. (2004). From the zonotope construction to the Minkowski addition of convex polytopes. Journal of Symbolic Computation, 38(4), 1261–1272. Fukuda, K., & Weibel, C. (2005). Computing all faces of the Minkowski sum of V -polytopes. In Proceedings of the 17th Canadian conference on computational geometry. Goldfarb, D., & Liu, S. (1991). An O(n3 ) primal interior point algorithm for convex quadratic programming. Mathematical Programming, 49, 325–340. Goodwin, G. C., Haimovich, H., Quevedo, D. E., & Welsh, J. S. (2004). A moving horizon approach to networked control system design. IEEE Transactions on Automatic Control, 49(9), 1427–1445. Grunbaum, B. (2003). Convex polytopes. New York: Springer. Ho, T. J. (2011). A switched IMM-extended Viterbi estimator-based algorithm for maneuvering target tracking. Automatica, 47(1), 92–98. Hwang, I., Balakrishnan, H., & Tomlin, C. (2006). State estimation for hybrid systems: applications to aircraft tracking. IEE Proceedings—Control Theory and Applications, 153(5), 556–566. James, R. (1988). A polynomial-time algorithm, based on Newtons method, for linear programming. Mathematical Programming, 40(1–3), 59–93. Jaulin, L. (2009). Robust set-membership state estimation: application to underwater robotics. Automatica, 45(1), 202–206. Karmarkar, N. (1984). A new polynomial-time algorithm for linear programming. Combinatorica, 4(4), 373–396. Khachiyan, L. (1979). A polynomial algorithm for linear programming. Soviet Mathematics Doklady, 20(1), 191–194. Kuntsevich, V., & Lychak, M. (1985). Synthesis of optimal and adaptive control systems: the game approach. Kiev: Naukova Dumka. Lewis, F. L., Xie, L. H., & Popa, D. (2008). Optimal and robust estimation: with an introduction to stochastic control theory. CRC Press. Liu, H. P., Ho, D. W. C., & Niu, Y. G. (2010). Robust filtering design for stochastic system with mode-dependent output quantization. IEEE Transactions on Signal Processing, 58(12), 6410–6416. Liu, H. P., Ho, D. W. C., & Sun, F. C. (2008). Design of H∞ filter for Markov jumping linear systems with non-accessible mode information. Automatica, 44(10), 2655–2660. Milanese, M., & Vicino, A. (1991). Optimal estimation theory for dynamic systems with set membership uncertainty: an overview. Automatica, 27(6), 997–1009. Mo, S. H., & Norton, J. P. (1990). Fast and robust algorithm to compute exact polytope parameter bounds. Mathematics and Computers in Simulation, 32(5–6), 481–493. Ong, C., & Gilbert, E. G. (2006). The minimal disturbance invariant set: outer approximations via its partial sums. Automatica, 42(9), 1563–1568. Orguner, U., & Demirekler, M. (2008). Risk-sensitive filtering for jump Markov linear systems. Automatica, 44(1), 109–118. Pan, Z., & Basar, T. (1996). H∞ control of Markovian jump systems and solutions to associated piecewise-deterministic differential games. In G. J. Olsder (Ed.), Annals of the international society of dynamic games (pp. 61–94). Boston, MA: Birkhäuser. Robert, C. P., & Casella, G. (2004). Monte Carlo statistical methods. New York: Springer. Schweppe, F. C. (1968). Recursive state estimation: unknown but bounded errors and system inputs. IEEE Transactions on Automatic Control, 13(1), 22–28. Serpedin, E., Chen, T., & Rajan, D. (2011). Mathematical foundations for signal processing, communications, and networking. CRC Press. Souza, D. C. E., Trofino, A., & Barbosa, K. A. (2006). Mode-independent H∞ filters for Markovian jump linear systems. IEEE Transactions on Automatic Control, 51(11), 1837–1841. Vicino, A., & Zappa, G. (1996). Sequential approximation of feasible parameter sets for identification with set membership uncertainty. IEEE Transactions on Automatic Control, 41(6), 774–785.

H. Wu et al. / Automatica 49 (2013) 3292–3303 Walter, E., & Piet-lahanier, H. (1989). Exact recursive polyhedral description of the feasible parameter set for bounded-error models. IEEE Transactions on Automatic Control, 34(8), 911–915. Wu, H., & Ye, H. (2012). State estimation for networked systems: an extended IMM algorithm. International Journal of Systems and Science, 44(7), 1274–1289. Yemelichev, V. A., Kovalev, M. M., & Kravtsov, M. K. (1984). Polytopes, graphs and optimization. Cambridge University Press. Ziegler, G. M. (1995). Lectures on polytopes. New York: Springer-Verlag.

Hao Wu was born in Shandong, China, in 1985. He received his B.S. degree in Mathematics from Shandong University, China, in 2008. He is now pursuing a Ph.D. in Automation from Tsinghua University, China. His research interests are state and parameter estimation.

Wei Wang received her Bachelor’s degree from the Beijing University of Aeronautics and Astronautics, China in July 2005, her Master’s degree from the University of Southampton, UK in December 2006 and her Ph.D. from Nanyang Technological University Singapore, in November 2011. She is currently a Lecturer with the Department of Automation at Tsinghua University, China. Her research interests include adaptive control, fault tolerant control, distributed cooperative control, state observer design and applications related to mechanical systems.

3303 Hao Ye was born in China, in 1969. He received his Bachelor’s and Doctoral degrees in Automation from Tsinghua University, China, in 1992 and 1996, respectively. He has been with the Department of Automation, Tsinghua University since 1996. He is currently a Professor of the Department of Automation, Tsinghua University. He is mainly interested in fault detection and diagnosis of dynamic systems.

Zidong Wang was born in Jiangsu, China, in 1966. He received his B.Sc. degree in mathematics in 1986 from Suzhou University, Suzhou, China, and his M.Sc. degree in applied mathematics in 1990 and his Ph.D. in electrical engineering in 1994, both from Nanjing University of Science and Technology, Nanjing, China. He is currently Professor of Dynamical Systems and Computing in the Department of Information Systems and Computing, Brunel University, UK. From 1990 to 2002, he held teaching and research appointments in universities in China, Germany and the UK. Prof. Wang’s research interests include dynamical systems, signal processing, bioinformatics, control theory and applications. He has published more than 200 papers in refereed international journals. He is a holder of the Alexander von Humboldt Research Fellowship of Germany, the JSPS Research Fellowship of Japan, William Mong Visiting Research Fellowship of Hong Kong. Prof. Wang serves as an Associate Editor for 11 international journals, including IEEE Transactions on Automatic Control, IEEE Transactions on Control Systems Technology, IEEE Transactions on Neural Networks, IEEE Transactions on Signal Processing, and IEEE Transactions on Systems, Man, and Cybernetics—Part C. He is a Senior Member of the IEEE, a Fellow of the Royal Statistical Society and a member of program committee for many international conferences.