- Email: [email protected]
Modal Estimation in Hybrid Systems
Journal of Mathematical Analysis and Applications 245, 225᎐247 Ž2000. doi:10.1006rjmaa.2000.6764, available online at http:rrwww.idealibrary.com on
Modal Estimation in Hybrid Systems David D. Sworder Department of Electrical and Computer Engineering, Uni¨ ersity of California, San Diego, California E-mail: [email protected]
John E. Boyd Cubic Defense Systems, Inc., San Diego, California E-mail: [email protected]
and Robert J. Elliott Department of Mathematics, Uni¨ ersity of Alberta, Edmonton, Alberta, Canada E-mail: [email protected] Communicated by C. T. Leondes Received September 13, 1999
In hybrid systems, estimation is made difficult by the nonlinear equations of state evolution. Several algorithms for generating an approximation to the conditional-mean estimator have been proposed. These are satisfactory in many cases, but because of the way the observations are fused, the estimate of the modal component of the state is often less good. This paper proposes a novel algorithm for state estimation and compares it with some in current use. 䊚 2000 Academic Press
1. INTRODUCTION Hybrid models naturally arise in applications in which a nonlinear plant operates sequentially about a collection of set-points. The plant model is phrased in terms of a set of nonlinear stochastic differential equations on the probability space Ž ⍀, F , P . and the time interval w0, T x. There are a right continuous filtration Ft ; 0 F t F T 4 , a right continuous, Ft-adapted 225 0022-247Xr00 $35.00 Copyright 䊚 2000 by Academic Press All rights of reproduction in any form reserved.
226
SWORDER, BOYD, AND ELLIOTT
random processes, ⌽t 4 , a piecewise constant regime process ranging over an index set of size S, and wt 4 , a Brownian motion plant disturbance.1 Subject to initial condition 0 , the plant model is d t s f Ž t , ¨ t , ⌽t . dt q g Ž t , ¨ t , ⌽t . dwt ,
Ž 1.
where ¨ t is an s-dimensional actuating signal and t is the n-dimensional plant-state. In this paper, we will ignore endogenous actuation: ¨ t ' 0. The values of the regime process are associated with stasis condition for the plant. The system operates near one set-point for an unpredictable interval, then suddenly shifts to another. The set of nominal plant-states can be arrayed as an Ž n = S .-matrix , with the nominal associated with ith mode given by i : s w i x. During operation, the system will sojourn in one regime for a time Ž ⌽t s ⌽i for t g w a, b .., and then suddenly shift Ž ⌽ b s ⌽j . to another in response to an external event or change in the surrounding environment. Despite the simple form of the regime process, the equation of plant evolution is intractable. As an approximation, Ž1. is often replaced with a family of local models created by linearizing the plant dynamics about the various set-points. Let t be a pointer to the current regime: The state space of t consists of the S canonical unit vectors in ⺢ s Ž t g e 1 , . . . , e S 4.. The component in t with value one marks the current mode of operation: If ⌽t s ⌽ k , then t s e k . The modal dynamics are usually represented by an exogenous Markov model with generator Q⬘. The base-state is the deviation of the plant-state from the current set-point: x t s t y t . The t 4 process is called the modal-state process to differentiate it from the base-state process. The hybrid-state of the system is the composite of the base- and modal-states. A comprehensive state model requires a representation of evolution, both intramodal and intermodal. During an extended modal sojourn, proper control will place and maintain the plant-state vector near the 1 For notational convenience in what follows: S will designate an integer index set 1, . . . , S4; e i is the ith canonical unit vector in a space whose dimension is obvious from the context; where no confusion will arise, a subscript may identify time or the component of a vector Že.g., 1 4 is a scalar process that is the first component of the vector process t 4.; ) is the Hadamard product ŽŽ x ) y . i s x i yi .; and m indicates Kronecker product. If a process is sampled T times per second, the discrete sequence so generated is written yw k x4 where the index denotes sample number rather than time. A process discontinuity is labeled ⌬ x t s x t y Žtime-continuous. or ⌬ x w k x s x w k x y x w k y 1x Žtime-discrete.. Conditional expectation is xy t denoted with a circumflex with the relevant -field apparent from context. A Gaussian random variable with mean ˆ x t and covariance Px x is labeled x g NŽ ˆ x t , Px x . with the same symbol used for the density function itself where no confusion will arise. Denote the density function of a unit Gaussian random vector by ⌽ Ž u.. If A is a positive symmetric matrix and x is a compatible vector, x⬘ Ax is denoted 5 x 5 2A .
MODAL ESTIMATION IN HYBRID SYSTEMS
227
correct set-point. The natural plant model in this circumstance would be that local model, selected from the family of regime-specific, linear models, associated with the current set point. But there will be a base-state discontinuity when the mode changes. If e i ¬ e l , the base-state reference level changes from i to l : ⌬ x t s y ⌬ t .
Ž 2.
Because of the wide applicability of hybrid systems, considerable effort has been devoted to problems of estimation. Both time-continuous and time-discrete models have been used, the latter formed by sampling the relevant variables of the former at rate 1rT. The time-discrete hybrid model accepts certain temporal coincidences; e.g., the modal state is constant between sample times, but if T is small enough, the resulting errors are small. The time-discrete hybrid model can be developed as is done in w11x, xw k x s
Ý Ž A j x w k y 1x q C j w w k x . eXj w k y 1x y ⌬ w k x
Ž 3.
j
w k x s ⌸ w k y 1x q mw k x ,
Ž 4.
where ww k x4 is a unit Gaussian white sequence. If the local models are controllable from the plant noise in every regime, for i g S, Ci and A i are both nonsingular. Label the local plant-noise covariance with R Ž i . s D Ž i .y1 s Ci CiX ) 0. The time-discrete modal transition probabilities are given by ⌸. In the intersample interval wŽ k y 1.T, kT ., the base-state evolves according to the mode at the beginning of the interval. If there is a modal transition, the base-state will experience the discontinuity y ⌬ w k x. Particularly when the value of t 4 is important, a measurement architecture in which the plant-state and the modal-state are assigned separate sensors has been shown to be effective. The sensors could be time-continuous or time-discrete, but the focus here will be on the latter. Suppose the plant-state sensor gives a linear, albeit mode dependent, measurement of the plant-state vector. Alternatively, the modal sensor is a classifier with quality represented by the S = S discernibility matrix D s wDi j x where Di j is the probability that model-state e i will be selected by the classifier if e j is the true mode at time of measurement. The measurement model is then given by yw k x s
Ý Hi eXi w k x w k x q nw k x
Ž 5.
i
z w k x s D w k x q w k x ,
Ž 6.
228
SWORDER, BOYD, AND ELLIOTT
where nw k x4 is a Gaussian white sequence with positive covariance R x s Dy1 s F⬘F. The columns of D are probability vectors, and the form is x flexible enough to include aliasing and bias. The measurement rate is determined by the sensor and its utilization policy. The measurements may occur at a different rate than the basic clock rateᎏusually much slower. Should the observation rate be slower, the sensor gains can be made time dependent so that at times of no measurement yw k x and z w k x are uninformative. This will be done in what follows without comment. The initial plant state categories are assumed to be independent with probability distributions x w0x ; NŽ ˆ x w0x, Px x w0x. and w0x ; ˆw0x. Denote the filtration generated by ww k x4 and mw k x4 along with the initial conditions on the plant by O w k x. Then ww k x4 and mw k x4 are O w k y 1x-martingale increment sequences. Denote the filtration generated by the plant measurements yw k x4 by Y w k x Žrespectively, that generated by the modal observation by Z w k x.. The observation filtration is G w k x s Y w k x k Z w k x. Let F w k y 1x s G w k y 1x k O w k x. The exogenous processes in Ž5. and Ž6. are F w k y 1x-martingale increments. Other authors study state estimation and regulation with an emphasis on time-continuous hybrid plants and a mix of time-continuous and timediscrete observations. Often, an engineer seeks practical algorithms for approximating the Gt-regime probabilities along with those Gt-moments Žincluding cross moments. important in the application. The polymorphic estimator Žthe PME. is a useful tool w10, 11, 14x. The PME is moment based and does not provide the Gt-distribution function of the hybrid-state. Further, an asymmetry in the data processing reduces the quality of the modal estimate. In this paper, we will provide an explicit processing algorithm for approximating the Gt-probability distribution of the plantstate t and the regime index ⌽t . This algorithm affords a more tightly integrated fusion of Y w k x and Z w k x than is provided by the PME, and it leads to improved identification of the regime process.
2. THE BAYES RECURRENCE FORMULA In contrast to the moment-based approach used to frame the PME, we will seek an implementable approximation to the G w k x-distribution of the hybrid-state. At time t s kT, suppose the G w k x-density is known: pw k x s w pi w k xx where pi w k x Ž z . dz s P Ž x w k x g w z, z q dz x , w k x s e i < G w k x . .
Ž 7.
We seek a mapping from Ž pw k x, yw k q 1x, z w k q 1x. to pw k q 1x. This mapping can be displayed most concisely when phrased in terms of an
229
MODAL ESTIMATION IN HYBRID SYSTEMS
unnormalized density, qw k xŽ z . s w qi w k xŽ z .x, derived using a reference probability. On the original event space and filtration, Ž ⍀, F ; F w k x., consider a reference probability measure P with respect to which mw k x4 and ww k x4 retain their earlier character, but yw k x4 is a unit Gaussian white sequence and z w k x4 is an independent, identically distributed Žiid. sequence that is uniformly distributed across e 1 , . . . , e S 4 . For l g 0, 1, . . . 4 , define
w l x s
S ⌽Ž yw l x.
z w l x ⬘ D w l x < F
= Fy1 y w l x y
ž ž
Ý Hi eXi w l x Ž e w l x q i . i
//
Ž 8.
and let ⌳w k x be a continuing product of the w l x: ⌳w k x s Ł 0k w l x. It is shown in w11x that P and P are related by
⭸P ⭸P
F w k xk G w k x
s ⌳w k x Ž zw k x , yw k x. .
Ž 9.
This change to P reflects no change for the plant disturbances, but the character of the observation is considerably different. Denote expectation with respect to P by E. The conditional Bayes theorem w3, Theorem 3.2, Chap. 2x relates the operator E to the operator E. Let f be a scalar test function. The expectation of eXi ⌽w k q 1x f Ž x w k q 1x. can be computed using either P or P E eXi w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x E ⌳ w k q 1x < G w k q 1x
.
Ž 10 .
Equation Ž10. shows that an expectation with respect to P can be performed as an expectation with respect to the reference measure P if the derivative of the measure is included properly in the calculation. The denominator has nothing to do with f and is simply a normalizing factor. The numerator will be thought of as the G w k q 1x-conditional expectation of f Ž x w k q 1x. with respect to an unnormalized G w k q 1x-conditional probability density, qi w k q 1xŽ x w k q 1x.: E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
H⍀ f Ž z . q w k q 1x Ž z . dz. i
Ž 11 .
230
SWORDER, BOYD, AND ELLIOTT
If qi w k q 1x; i g S4 were known, the joint density = mass function of the hybrid-state would be pi w k q 1 x Ž z . s
qi w k q 1 x Ž z . Ý j H⍀ q j w k q 1 x Ž u . du
,
Ž 12 .
from which ˆ x w k q 1x and ˆw k q 1x along with higher moments could be calculated. To find qi w k q 1xŽ z .; i g S4 , expand Ž11. by replacing w k q 1x, ⌳w k q 1x, and x w k q 1x with their values. E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s E eXi Ž ⌸ w k x q m w k q 1 x . ⌳ w k x z w k q 1 x ⬘De i < F
žÝŽ
S ⌽ Ž y w k q 1x .
A j x w k x q C j w w k q 1 x . eXj w k x y Ž e i y w k x . < G w k q 1 x .
/
j
Ž 13 . It is true that Ý j w k x⬘e j ' 1. Substituting this into Ž13., E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x sE
Ý Ž w k x ⬘e j . ŽeXi ⌸ e j . ⌳ w k x z w k q 1x ⬘De i < F
⭈⌽ Ž Fy1 Ž y w k q 1 x y Hi Ž x w k q 1 x q e i . . . ⭈
S ⌽ Ž y w k q 1x .
f Ž A j x w k x q C j w w k q 1x y Ž e i y e j . . < G w k q 1x .
Under P, G w k x4 is uninformative with respect to the hybrid-state. Since an unnormalized distribution is sought, factors common to all regimes can and will be ignored in what follows. E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
H⍀ Ý ⌸
ij z
w k q 1 x ⬘D .i ⌽
j
⭈ Fy1 Ž y w k q 1 x y Hi Ž A j q C j w y Ž e i y e j . q e i . .
ž
⭈ f Ž A j q C j w y Ž e i y e j . . q j w k x Ž . ⌽ Ž w . d dw.
/
MODAL ESTIMATION IN HYBRID SYSTEMS
231
To simplify this, make the change of variable z s A j q C j w y i q j . Then d dw s < C j
H⍀ Ý ⌸
ij z
w k q 1 x ⬘D .i < C j
j
⭈ f Ž z . q j w k x Ž . ⌽ Ž Cy1 j Ž z y A j q i y j . . d dz s
H⍀ f Ž z . q w k q 1x Ž z . dz. i
Hence Observation update qi w k q 1xŽ z . s Ý ⌸ i j z w k q 1x⬘D.i < C j
⭈
H⍀⌽ Ž C
y1 j
Ž z y A j q i y j . . q j w k xŽ . d .
Ž14.
Equation Ž14. is the recurrence formula for the unnormalized density of the hybrid-state. Variants on Ž14. have been derived for kindred plant representations: If s 0 and the noise in z w k x4 is Gaussian, see w4, Eq. Ž8.x. If s 0 and there is a single measurement g w k x s yw k x q z w k x in Gaussian noise, see w5, Eq. Ž10.x. If s 0 and there is a slight indexing difference in the base-state model, see w5, Eq. Ž4.x. 䢇 䢇
䢇
Unfortunately, Ž14. is not an algorithm of the form we seekᎏalthough recursive, it is infinite dimensional.
3. A GAUSSIAN WAVELET ESTIMATOR We seek a finite dimensional, Gaussian wavelet approximation to qw k x. The number of terms is limited only by the computational complexity permitted in the application. We will assume that N terms will suffice in every regime. The resulting algorithm will be called the Gaussian wavelet estimator ŽGWE.. For notational convenience, let the difference in the ith and the jth set-points be written ji s i y j , and let the observation
232
SWORDER, BOYD, AND ELLIOTT
centered about the ith set-point be written y i w k q 1x s yw k q 1x y Hi i . Denote the inverse of a Žpositive. covariance matrix P by D, e.g., D li w k x s Pli w k xy1 , and the product Dm by d, e.g., d li w k x s D li w k x m il w k x. Suppose qj w k x Ž . s
N
Ý ␣ lj w k x
D lj w k x
1r2
1
exp y
2
ls1
y m lj w k x
2 D lj w k x .
Ž 15 .
In Ž15., the unnormalized G w k x-conditional probability of the event x w k x g w , q d x, w k x s e j 4 is given under P by a sum of N Gaussian pattern functions, N Ž m lj w k x, Pl j w k x.; l g N4 with means m lj w k x and positive covariances Pl j w k x. In this approximation, m lj w k x translates the lth pattern function, Pl j w k x adjusts its shape, and N circumscribes the span of the sum. The pattern functions are weighted by ␣ lj w k x; l g N4 . All of the coefficients are G w k x-adapted. There are NS elements in qw k x though many could be zero. The Appendix presents a recurrence formula for the coefficients ␣ lj w k x4 , m lj w k x4 , Pl j w k x4 Ž j g S, l g N.. The recurrence can be most concisely stated in a mixed covariance-information form also used in the Kalman filter w1x. The information matrix is the inverse of the covariance matrix, and instead of framing the estimator using m il 4 as coordinates, the information filter propagates d li w k x. In this mixed coordinate system, the base-state estimate satisfies base-state recurrence: GWE Extrapolation: m il w k q 1xys A i m il w k x; Pli w k q 1xys
A i Pli w k x AXi
i g S,
lgN
q R Ž i .;
i g S,
Ž16. lgN
Ž17.
Update: d i w k q 1xŽ l; j .qs d lj w k q 1xyq HiX Dx y i w k q 1x yD lj w k q 1xy ji ;
i, j g S,
D i w k q 1xŽ l; j .qs D lj w k q 1xyq HiX Dx Hi ;
Ž18. lgN
i, j g S,
l g N.
Ž19.
As in many multiple model algorithms, Ž16. ᎐ Ž19. require a filter bank of size S for extrapolation and update. But in the GWE, each filter extrapolates N initial conditions. After adjusting for the inclusion of the set points, the form of the GWE is close to that of the Kalman filter.
233
MODAL ESTIMATION IN HYBRID SYSTEMS
Under the nontransition hypothesis, Že i ¬ e i ., ⌬ d li w k q 1xqs HiX Dx y i w k q 1x, the Kalman update using the centered observation. Under the transition hypothesis Že j ¬ e i ., the increment in d i w k q 1xŽ l; j .q contains y i w k q 1x as a factor but the translation of the set-point, Žy i q j ., must be included. The weightings in d i w k q 1xŽ l; j .q balance measurement fidelity and the shape of the pattern functions. The size of jump q
d i w k q 1 x Ž l ; j . y d lj w k q 1 x
y
in the Ž l; j .th subfilter will be increased if the measurement noise is small and if yw k q 1x differs from Hi i . Similarly, the set-point offset Ž j y i . is important when the probability of the lth pattern function is concentrated in a small region, i.e., when D lj w k q 1xy is big. Suppose the modal measurement is quite good and an e j ¬ e i transition occurs at time t s Ž k q 1.T. In the absence of the base-state measurement, the extrapolation in the base-state mean should be q
y
m i w k q 1 x Ž l ; j . f m lj w k q 1 x q j y i . The mean translates by the amount of the base-state discontinuity. The size of the difference, m i w k q 1xŽ l; j .qy m lj w k q 1xyq ji, is a coarse measure of the influence of yw k q 1x since this difference is that between the corrected and uncorrected extrapolation. Modal mixing is distributed in the GWE. It appears in the filter update and also in the weighting coefficients, ␣ i w k q 1xŽ l; j .q; i, j g S, l g N4 . modal mixing: GWE
␣ i w k q 1xŽ l; j .qs ␣ lj w k x⌸ i j z w k q 1x⬘D.i ⭈ D lj w k x Ž D lj w k x q AXj D Ž j . A j . ⭈exp
1 2
½
m i w k q 1 x Ž l; j . y
y1
q
P i w k q 1 x Ž l; j . D Ž j .
1r2
q 2 q D i w kq 1 x Ž l; j .
y m lj w k q 1 x y ji
2 y w kq 1 x y D lj
y i w k q 1x
2 Dx
5.
Ž20.
The weighting coefficient ␣ i w k q 1xŽ l; j .q is proportional to its predecessor Ž ␣ lj w k x., proportional to the probability of an e j ¬ e i transition Ž ⌸ i j ., and proportional to the probability that the modal observation was generated from w k q 1x s e i Žthe factor z w k q 1x⬘D.i ..
234
SWORDER, BOYD, AND ELLIOTT
With coefficients thus determined, the unnormalized G w k q 1x-wavelet distribution of the hybrid-state can be written immediately: q
q
q
qi w k q 1 x s Ý ␣ i w k q 1 x Ž l ; j . N Ž m i w k q 1 x Ž l ; j . , P i w k q 1 x Ž l ; j .
q
..
j, l
Ž 21 . From Ž21., the various moments of interest can be computed; e.g.,
ˆi w k q 1 x s
Ý ␣ i w k q 1x Ž l ; j .
q
j, l
Ý ␣ i w k q 1x Ž l ; j .
q
Ž 22 .
i , j, l
and q
ˆx w k q 1 x s Ý ␣ i w k q 1 x Ž l ; j . m i w k q 1 x Ž l ; j . i , j, l
q
Ý ␣ i w k q 1x Ž l ; j .
q
.
i , j, l
Ž 23 . Unfortunately, there are not N terms in Ž21. but NS. To satisfy the complexity constraint, pruning andror merging must be used to reduce the number of terms. For now we will reduce the terms from each progenitor to one by using the conventional Gaussian sum merging formula given in w1x. For a given i g S, consider the jth progenitor. Define
␣ ji w k q 1 x Ž l . s ␣ i w k q 1 x Ž l ; j .
q
Ý ␣ i w k q 1x Ž l ; j .
q
.
Ž 24 .
l
Then let m ij w k q 1 x s
Ý mi w k q 1x Ž l ; j .
q
␣ ji w k q 1 x Ž l .
Ž 25 .
l
Pji w k q 1 x s
Ý ž P i w k q 1x Ž l ; j .
q
l
q
q m i w k q 1 x Ž l ; j . y m ij w k q 1 x
ž
q
/
= m i w k q 1 x Ž l ; j . y m ij w k q 1 x ⬘ ␣ ji w k q 1 x Ž l . Ž 26 .
ž
␣ ji w k q 1 x s
q Ý ␣ i w k q 1x Ž l ; j . .
l
//
Ž 27 .
235
MODAL ESTIMATION IN HYBRID SYSTEMS
The next iteration of the GWE begins with a set of S terms: N s S. updated distribution qi w k q 1x s
Ý ␣ jiw k q 1xNŽ mij w k q 1x, Pjiw k q 1x..
Ž28.
j
In many applications, there is a common set-point and, though the system dynamics change, the reference level does not. Such systems are called linear-jump-systems Žsee w8x. and are modeled with s 0. The GWE-algorithm is considerably simplified in this case. base-state recurrence: GWE Ž s 0. Extrapolation: m il w k q 1xys A i m il w k x;
i, l g S
Ž29.
Pli w k q 1xys A i Pli w k x AXi q R Ž i .;
i, l g S
Ž30.
Update: ⌬ d i w k q 1xŽ l; j .qs HiX Dx yw k q 1x; ⌬ D i w k q 1xŽ l; j .qs HiX Dx Hi ;
i, j, l g S
Ž31.
i, j, l g S,
Ž32.
where ⌬ d i w k q 1xŽ l; j .qs d i w k q 1xŽ l; j .qy d lj w k q 1xy, and similarly for ⌬ D i w k q 1xŽ l; j .q. The weighting coefficients can be simplified too. For i, j, l g S,
␣ i w k q 1xŽ l; j .qs ␣ lj w k x⌸ i j z w k q 1x⬘D.i < D lj w k xŽ D lj w k x q AXj D Ž j . A j .y1 ⭈P i w k q 1xŽ l; j .q < 1r2 exp 12 ⌬ m i w k q 1xŽ l; j .q 5 2D i w kq1xŽ l; j.q.
4. MULTIPLE MODEL ALGORITHMS The simplest of the GWEs are those for systems with a single set point and N s 1. Then qj w k x Ž . s ␣ j w k x D j w k x
1r2
exp y
1 2
y m jw k x
2 D jw k x .
Ž 33 .
236
SWORDER, BOYD, AND ELLIOTT
The extrapolation equations of GWE Ž s 0, N s 1. are direct restrictions of Ž29. ᎐ Ž32.. The extrapolation uses S parallel filters with some mixing in the update. base-state recurrence: GWE Ž s 0. Extrapolation: m i w k q 1xys A i m i w k x; igS p i w k q 1xys A i P i w k x AXi q R Ž i .;
Ž34. Ž35.
igS
Update: ⌬ d i w k q 1xŽ j .qs HiX Dx yw k q 1x; i, j g S ⌬ D i w k q 1xŽ j .qs HiX Dx Hi ; i, j g S
Ž36. Ž37.
The weighting coefficients can be written q
␣ i w k q 1 x Ž j . s ␣ j w k x ⌸ i j z w k q 1 x ⬘D .i D j w k x Ž D j w k x q AXj D Ž j . A j . ⭈P i w k q 1 x Ž j .
q
1r2
exp 12 ⌬ m i w k q 1 x Ž j .
q 2 D i w kq1 xŽ j .q .
y1
Ž 38 .
The unnormalized modal distribution is
␣ i w k q 1 x s z w k q 1 x ⬘D .i
Ý ␣ i w k q 1x Ž j .
q
.
Ž 39 .
j
For systems with a single set-point, several estimation agorithms have been proposed in the literature. That receiving most attention is the interacting-multiple-model ŽIMM. estimator w2x. As usually implemented, the IMM generates a Gaussian sum approximation to the Y w k x-density of the hybrid-state in the form Ý j eXj ˆw k xNŽ m ˆ j w k q 1x, P j w k q 1x.. The equations of implementation for the IMM are provided in several places Že.g., see w7, Table 1x., and an annotated description of applications of the IMM is found in w9x. The IMM was modified to accept the z w k x4 observations in w6x, and this G w k x-IMM was called the image-enhanced IMM. The steps in the implementation of the IMM are given explicitly in w6, Fig. 3x. This implementation of the IMM is closest to the GWE and will be referred to simply as the IMM. The motivation behind the IMM is the efficient computation of a high quality approximation to ˆ x w k x4 . It is not claimed that the IMM produces the G w k x-joint density of the hybrid-state, nor does it. Approximations created by the merging of the subfilters and the update rule preclude
MODAL ESTIMATION IN HYBRID SYSTEMS
237
accurate calculation of ˆw k x4 . Tests have shown that a high quality estimate of the modal-state is not required for accurate base-state estimation because the IMM is essentially a variable bandwidth filter for the latter. The IMM and the GWE Ž s 0, N s 1. perform common operations on the data. The IMM does hypothesis mixing at the beginning of an iteration while the GWE Ž s 0, N s 1. does it at then end. The factor z w k q 1x⬘D.i appears in both along with factors akin to Ý j ␣ j w k x⌸ i j Žafter adjusting for the time of hypothesis merging.. Both require S Kalmanextrapolations. Moreover, the IMM has S Kalman-updates while the GWE Ž s 0, N s 1. has S 2 . There are some differences in indexing because the IMM supposes that the evolution of x t 4 on wŽ K y 1.T, KT . depends on w k x while the GWE supposes a w k y 1x dependence. For small T, this matters little. The fundamental difference between the algorithms arises from the modal update. The ith filter for either the IMM or the GWE provides m i w k q 1xy as shown in Ž34.. Define the ith residual by r i w k q 1x s yw k q 1x y Hi m i w k q 1xy For a system with a linear-Gauss ᎐Markov ŽLGM. model Ž S s 1., r i w k x is the increment in the Y w k x-innovation processᎏ r i w k x4 is a Y w k x-martingale increment sequence. For the LGM estimation problem, the covariance of r i w k q 1x is y
Pyi y w k q 1 x s Hi P i w k q 1 x HiX q R x Ž i . .
Ž 40 .
The IMM uses multiple models, and the update of the modal probabilities depends upon the behavior of the subfilter residuals. The smaller the residuals in a particular subfilter are, weighted by the associated information matrix, the more likely the associated mode is thought to be. The unnormalized modal update takes the form y
q w k x s ˆ w k x ) Dyi y w k q 1 x
y 1r2
=exp y 12 r i w k q 1 x ⬘Dyi y w k q 1 x r i w k q 1 x 4 . Ž 41 . In contrast, the elements from which the GWE generates the modal estimate Žsee Ž38.. are responsive, not to the residuals themselves, but rather to the effect these residuals have on the update of the mean. The weight ␣ i w k q 1xŽ l; j .q is generated as a continuing exponential product of the quadratic differences of the base-state mean before and after a base-state measurement: ˆw k x depends not on the residuals themselves, but on the influence the residuals have on the base-state estimate.
238
SWORDER, BOYD, AND ELLIOTT
The updates of the IMM and GWE Ž s 0, N s 1. are closely related in certain cases. Denote the update gain of the ith filter in the IMM by W i and suppose the W i 4 are all nonsingularᎏthis unfortunately requires that the Hi 4 are all nonsingular. Then r i w k q 1x s ŽW i .y1 ⌬ m i w k q 1x. Expanding the quadratic form in Ž41. r i w k q 1x
2 D iy y w kq1 x
y
i s ⌬ m i w k q 1 x ⬘D i w k q 1 x Hy1 i Py y w k q 1 x
⭈ HyT D i w k q 1x i
y
⌬ mi w k q 1x y
s ⌬ m i w k q 1 x ⬘D i w k q 1 x Hy1 i y
⭈ Hi P i w k q 1 x HiX q R Ž i .
ž
⭈ HyT D i w k q 1x i
y
/
⌬ mi w k q 1x .
If R Ž i . is small Žas occurs during intermodal transients . and if the base-state observation gain is nonsingular Žan unusual event. r i w k q 1x
2 D iy y w kq1 x
f ⌬ mi w k q 1x
2 D i w kq1 xy
/ ⌬ mi w k q 1x Ž j .
q 2 D i w kq1 xŽ j .q .
Ž 42 .
The IMM and GWE Ž s 0, N s 1. have similar modal updates, but the IMM weights a quadratic form in the increments of the estimate while the GWE weights an increment in the quadratic forms. The update of the conditional mean for GWE Ž s 0, N s 1. involves more bookkeeping than does the similar operation IMM. Modal mixing takes place even in the absence of set-point discontinuities.
5. THE GWE IN ESTIMATION AND IDENTIFICATION To illustrate the performance of the GWE, consider a specific example. On the California desert stands a 10 MWe solar electric generating system. Movable mirrors Žheliostats . are used to focus the sun’s energy on a group of boiler panels atop a central tower. A steam temperature regulator controls the feedwater flow rate to the panels so as to maintain the proper temperature of both the outlet steam and receiver panel itself w13x. A remote observer wishes to estimate the panel temperature and identify the current set-point. On a partly cloudy day the insolation changes suddenly and unpredictably. The panel dynamics change as well: at low insolation, the panel dynamics are slow; at high insolation, the panel dynamics are fast. Let Tm be the panel temperature and let Tc be the
239
MODAL ESTIMATION IN HYBRID SYSTEMS
output of a temperature sensor, with m Žrespectively c . the perturbation variables. On a day where panel operation moves between sun and cloud, the local models of the panel and the sensor are given by Panel model dm s ␣ m dt q  dw1 ,
Ž 43 .
dTc s a Ž Tc y Tm . dt q 0.3  dw 2 .
Ž 44 .
Sensor model
The coefficients of the local panel models are given as a function of insolation in w12x and are displayed in Table I. The nominal panel temperature at high insolation is 17⬚F higher than that for low insolation: s Ž1017, 1000. m 1. The state ordering in what follows is x s vecŽm , c .. To particularize the modal dynamics, Q must be specified. The mean sojourn of bright sunlight is 2 min while clouds persist for an average of 4 min. A simple Markov model matches these mean sojourns. The sensors at the site are of two types: a direct measurement of Tc in noise and an insolation measurement. For the purposes of this example, the former will be sampled every 0.5 min in white Gaussian noise of standard deviation 2⬚F: yw k x s D w k x q nw k x, where D s eX2 and R x s 4. The panel measurement is of absolute temperature, not the perturbation variable, c . To improve performance, the insolation sensor gives a direct measurement of 4 . In the application, the insolation is not uniformly distributed across the panel: the geometry of the sun-mirror-panel is such that the distribution of the solar flux across a panel is irregular with possible hot spots. Insolation sensors measure conditions only in a neighborhood. The application employed several insolation sensors arrayed across the panel, and the highest reading was used as a measure of panel insolation. The insolation sensor is sampled every 0.1 min. In w12x several panel temperature estimators were contrasted. An EKF using pseudonoise to represent the modal transitions was a failure. When the discernibility matrix indi-
TABLE I Local Panel Temperature Dynamics for Two Insolation Conditions Modal-state
Insolation
¨n
n
␣

e1
1.8 = 10 BTUrmin 186 lbrmin
1017⬚F y1.8
e2
2.7 = 10 BTUrmin
1000⬚F y0.36 1.9
5
4
30 lbrmin
a
1.58 0.9 0.25
240
SWORDER, BOYD, AND ELLIOTT
cated high quality modal measurements, PME was effective for both panel temperature estimation and modal identification. The performance of the GWE in this circumstance is comparable. A lower quality D created difficulties for the PME because the hypotheses underlying the algorithm became more difficult to justify. Suppose the less accurate insolation sensor introduced in w12x is used in this application: D1 s
0.65 0.35
0.40 . 0.60
Ž 45 .
The inclination toward high insolation readings is accommodated with this asymmetrical D. The precision of the classifier is such that the modal-state is difficult to identify with confidence. Figure 1 shows sample functions of Tm and Tc for an insolation path: has the value e1 on t g w0, 1. j w6, 9. and the value e2 on w1, 6. j w9, 11x. Both temperatures were initialized at 1017⬚F. After the cloud comes by at t s 1 min, Tm moves rapidly toward 1000⬚F, with Tc moving at a much slower rate. A 6⬚F temperature difference between the sensor and the panel develops in 5 min. When the sun returns, the sensor output has significant lag. A sample observation sequence, yw k x4 , is also shown and is marked on the figure with an ‘‘>.’’
FIG. 1. A sample path of the panel temperature Tm and its measurement Tc . Noisy temperature measurements are shown every 0.5 min.
MODAL ESTIMATION IN HYBRID SYSTEMS
241
An estimator that attempted to determine t 4 from yw k x4 alone would have a difficult time. Even determining Tm from such an abbreviated data set is a difficult task indeed. Figure 2 shows a plot of 1 4 Žsojourns in the sun., the insolation measurements Žlabeled >., and probability of sun as computed by GWE. The poor quality of the insolation sensor is reflected in a slow response to modal transition and a generally low confidence in the modal identification. Even so the modal identification is superior to that achieved by the PME in w12x. The performance of the panel temperature estimator is shown in Fig. 3. True temperature is shown along with the observations. The GWE initializes itself in the first minute. The GWE follows the metal temperature, but with some lag. The GWE is superior to the PME when the insolation sensor is of such low quality. Slight improvement in the modal measurements yields much improved set-point identification. The discernibility matrix D 2 is a higher quality insolation sensor, and the modal estimates shown in Fig. 4 reflect this:
D2 s
0.75 0.25
0.30 . 0.70
Ž 46 .
FIG. 2. The probability of sun as computed by GWE along with the mode. The insolation measurements are denoted with an >.
242
SWORDER, BOYD, AND ELLIOTT
FIG. 3. Estimate of the panel temperature using the GWE and sensor D1 .
FIG. 4. The probability of sun as computed by GWE along with the true mode. The insolation measurements are denoted with an >.
243
MODAL ESTIMATION IN HYBRID SYSTEMS
6. CONCLUSION The GWE is a novel, finite dimensional, approximation to the G w k x density of the hybrid-state. It provides the solution in the form of a Gaussian sum. In contrast to the PME, the observation sequences yw k x4 and z w k x4 are fused in a symmetric manner. When the set-point reference is constant, the GWE is considerably simpler. Indeed, GWE Ž s 0, N s 1. is about as simple as the IMM. The GWE will accommodate variable set-points, a subtle problem illustrated in the example.
APPENDIX At time t s kT, qw k x is a Gaussian sum with N terms in each component: N
qj w k x Ž . s
Ý ␣ lj w k x
D lj w k x
1r2
exp y
y m lj w k x
1 2
2 D lj w k x .
ls1
The recurrence relation for qi w k q 1x is given in Ž14.. qi w k q 1 x Ž z . s
Ý ⌸ i j z w k q 1x ⬘D.i < C j
H⍀ ⌽ ž C Ž z y A q . / q w k x Ž . d .
⭈
y1 j
j
i j
j
Expand the product of Gaussian pattern functions. In what follows, factors common across all modal hypotheses will be ignored without comment; i ⌽ Ž Fy1 Ž y i w k q 1 x y Hi z . . ⌽ Cy1 j Ž z y A j q j . q j w k x Ž .
ž
N
s
Ý ␣ lj w k x
D lj w k x
1r2
/
exp y 12 J1 Ž l . ,
ls1
where J1 Ž l . s y m lj w k x
2 D lj w k x
i q Cy1 j Ž z y A j q j .
q Fy1 Ž y i w k q 1 x y Hi z .
2
.
2
244
SWORDER, BOYD, AND ELLIOTT
The argument of the exponential is a quadratic form. It is an exercise to complete the squares in and z. Define a set of coefficient matrices: F j w k x Ž l ; i . s D Ž j . q HiX Dx Hi H j s AXj D Ž j . A lj w k x s D lj w k x q AXj D Ž j . A j E j w k x Ž l ; i . s D lj w k x m lj w k x q H j ji K j w k x Ž l ; i . s HiX Dx y i w k q 1 x y D Ž j . ji G j w k x Ž l ; i . s m lj w k x
2 D lj w k x
q ji
2
q y i w k q 1x
D Ž j .
2 Dx .
Completing the squares and combining exp y 12 J1 Ž l . s A lj w k x
y Ž1r2 .
⭈ exp y
1 2
D i w k q 1x Ž l ; j .
½G w k x Ž l ; i . y j
y Ž1r2 .
E jw k xŽ l; i.
y d i w k q 1x Ž l ; j . ⭈ N A lj w k x
ž
y1
2 A lj w k xy 1
2 P i w kq1 xŽ l ; j .
Ž E j w k x Ž l ; i . y Hzj . , A lj w k x y1 /
⭈ N z Ž mi w k q 1x Ž l ; j . , P i w k q 1x Ž l ; j . . , where D i w k q 1 x Ž l ; j . s F j w k x Ž l ; i . y H j⬘A lj w k x s P i w k q 1x Ž l ; j .
y1
y1
Hj
and
mi w k q 1x Ž l ; j . s P i w k q 1 x Ž l ; j . Ž H j⬘A lj w k x
y1
E jw k xŽ l; i. q K jw k xŽ l; i. . .
Now integrate
H⍀ ⌽ Ž F Ž y w k q 1x y H z . . ⌽ ž C Ž z y A q . / q w k x Ž . d y1
N
s
i
Ý ␣ lj w k x
i
D lj w k x
1r2
y1 j
A lj w k x
j
y Ž1r2 .
ls1
⭈ F j w k x Ž l ; i . y Ž H j . ⬘ A lj w k x
y1
Hj
y Ž1r2 .
i j
j
245
MODAL ESTIMATION IN HYBRID SYSTEMS
⭈ exp y
1 2
½G w k x Ž l ; i . y j
E jw k xŽ l; i.
y d i w k q 1x Ž l ; j .
2 A lj w k xy 1
2 P i w kq1 xŽ l ; j .
5
⭈ N z Ž mi w k q 1x Ž l ; j . , P i w k q 1x Ž l ; j . . . Then qi w k q 1x is a Gaussian sum with 1r2
␣ i w k q 1 x Ž l ; j . s ⌸ i j z w k q 1 x ⬘D .i < C j
y Ž1r2 .
1 2
D i w k q 1x Ž l ; j .
½G w k x Ž l ; i . y j
y Ž1r2 .
E jw k xŽ l; i.
y mi w k q 1x Ž l ; j .
2 A lj w k xy 1
2 D i w kq1 xŽ l ; j .
5.
The expressions for the coefficients can be simplified. Note first that D i w k q 1 x Ž l ; j . s D Ž j . q HiX Dx Hi y D Ž j . A j Ž D lj w k x q AXj D Ž j . A j . s HiX Dx Hi q Ž R Ž j . q A j Pl j w k x AXj .
y1
y1
AXj D Ž j .
.
This can be written in more conventional form as y
Pl j w k q 1 x s A j Pl j w k x AXj q R Ž j . followed by y
D i w k q 1 x Ž l ; j . s D lj w k q 1 x q HiX Dx Hi . Also d i w k q 1x Ž l ; j . s D Ž j . A j Ž D lj w k x q AXj D Ž j . A j .
y1
d lj w k x
q D Ž j . A j Ž D lj w k x q AXj D Ž j . A j .
ž
q
HiX Dx
y1
AXj D Ž j . y D Ž j . ji
y i w k q 1x .
Write this as y
d lj w k q 1 x s D Ž j . A j Ž D lj w k x q AXj D Ž j . A j . y
y1
d lj w k x y
d i w k q 1 x Ž l ; j . s d lj w k q 1 x q HiX Dx y i w k q 1 x y D lj w k q 1 x ji .
/
246
SWORDER, BOYD, AND ELLIOTT
Note that
Ž Dlj w k x q AXj D Ž j . A j .
y1
yT s AyT I y Ž Ay1 q Pl j w k x . j j R Ž j . A j
ž
y1
yT Ay1 . j R Ž j . A j
/
So y
Pl j w k q 1 x D Ž j . A Ž D lj w k x q AXj D Ž j . A j . yT y Ž Ay1 q Pl j w k x . j R Ž j . A j
y1
y1
R Ž j . AyT j
yT yT s A j Ž Ay1 q Pl j w k x . y A j Ay1 j R Ž j . A j j R Ž j . A j
s A j Pl j w k x . Hence y
y
Pl j w k q 1 x d lj w k q 1 x s A j Pl j w k x d lj w k x s A j m lj w k x . This can be written y
m lj w k q 1 x s A j m lj w k x . The expression for the weighting coefficients is complicated in appearance. The argument in the exponential factor in ␣ i w k q 1xŽ l; j . can be rewritten: d lj w k x
2 P lj w k x
q 5 ji 5 2D Ž j. q y i w k q 1 x
y d lj w k x q AXj D Ž j . ji
2 A lj w k xy1
2 Dx
y d i w k q 1x Ž l ; j .
2 P i w kq1 xŽ l ; j . .
But d lj w k x
2 P lj w k x
y d lj w k x
2 A lj w k xy 1 y
s d lj w k x ⬘Pl j w k x AXj Pl j w k q 1 x A j Pl j w k x d lj w k x s A j m lj w k x
2 D lj w kq1 xy
.
Also 5 ji 5 2D Ž j. y AXj D Ž j . ji
2 A lj w k xy1
s 5 ji 5 2D lj w kq1xy .
MODAL ESTIMATION IN HYBRID SYSTEMS
247
Finally, y2 d lj w k x ⬘ A lj w k x
y1
AXj D Ž j . ji s 2 Ž d lj w k q 1 x
y
. ⬘ ji y y s y2 Ž m lj w k q 1 x . ⬘D lj w k q 1 x ji .
So, with some cancellation of terms, the result follows.
REFERENCES 1. B. D. O. Anderson and J. B. Moore, ‘‘Optimal Filtering,’’ Prentice Hall, New York, 1979. 2. H. A. P. Blom and Y. Bar-Shalom, The interacting multiple model algorithm for systems with Markov switching coefficients, IEEE Trans. Automat. Control 3, No. 8 Ž1988., 780᎐783. 3. R. J. Elliott, L. Aggoun, and J. B. Moore, ‘‘Hidden Markov Models: Estimation and Control,’’ Springer-Verlag, BerlinrNew York, 1995. 4. R. J. Elliott, F. Dufour, and D. D. Sworder, Exact hybrid filters in discrete time, IEEE Trans. Automat. Control. 41, No. 12 Ž1996., 1807᎐1810. 5. R. J. Elliott and J. van der Hoek, A finite dimensional filter for hybrid observations, IEEE Trans. Automat. Control AC-43 Ž1998., 736᎐739. 6. J. S. Evans and R. J. Evans, Image-enhanced multiple model tracking, Automatica 35 Ž1999., 1769᎐1786. 7. X. R. Li and Y. Bar-Shalom, Design of an interacting multiple model algorithm for air traffic control tracking, IEEE Trans. Control Systems Technol. 1, No. 3 Ž1993., 186᎐194. 8. M. Mariton, ‘‘Jump Linear Systems in Automatic Control,’’ Dekker, New York, 1990. 9. E. Mazor, A. Averbuch, Y. Bar-Shalom, and J. Dayan, Interacting multiple model methods in target tracking: a survey, IEEE Trans. Aerospace Electron. Systems 34, No. 1 Ž1998., 103᎐123. 10. D. D. Sworder and J. E. Boyd, Jump-diffusion processes in trackingrrecognition, IEEE Trans. Signal Processing 46, No. 1 Ž1998., 235᎐239. 11. D. D. Sworder and J. E. Boyd, ‘‘Estimation Problems in Hybrid Systems,’’ Cambridge Univ. Press, Cambridge, UK, 1999. 12. D. D. Sworder and J. E. Boyd, State estimation with biased observations, IEEE Trans. Systems Man Cybernet. A 29, No. 6 Ž1999.. 13. D. D. Sworder and R. O. Rogers, An LQ-solution to a control problem associated with a solar thermal central receiver, IEEE Trans. Automat. Control 28, No. 10 Ž1983., 971᎐978. 14. D. D. Sworder and R. Vojak, Hybrid estimation algorithms I, J. Optim. Theory Appl. 81, No. 7 Ž1994., 143᎐167.
Modal Estimation in Hybrid Systems David D. Sworder Department of Electrical and Computer Engineering, Uni¨ ersity of California, San Diego, California E-mail: [email protected]
John E. Boyd Cubic Defense Systems, Inc., San Diego, California E-mail: [email protected]
and Robert J. Elliott Department of Mathematics, Uni¨ ersity of Alberta, Edmonton, Alberta, Canada E-mail: [email protected] Communicated by C. T. Leondes Received September 13, 1999
In hybrid systems, estimation is made difficult by the nonlinear equations of state evolution. Several algorithms for generating an approximation to the conditional-mean estimator have been proposed. These are satisfactory in many cases, but because of the way the observations are fused, the estimate of the modal component of the state is often less good. This paper proposes a novel algorithm for state estimation and compares it with some in current use. 䊚 2000 Academic Press
1. INTRODUCTION Hybrid models naturally arise in applications in which a nonlinear plant operates sequentially about a collection of set-points. The plant model is phrased in terms of a set of nonlinear stochastic differential equations on the probability space Ž ⍀, F , P . and the time interval w0, T x. There are a right continuous filtration Ft ; 0 F t F T 4 , a right continuous, Ft-adapted 225 0022-247Xr00 $35.00 Copyright 䊚 2000 by Academic Press All rights of reproduction in any form reserved.
226
SWORDER, BOYD, AND ELLIOTT
random processes, ⌽t 4 , a piecewise constant regime process ranging over an index set of size S, and wt 4 , a Brownian motion plant disturbance.1 Subject to initial condition 0 , the plant model is d t s f Ž t , ¨ t , ⌽t . dt q g Ž t , ¨ t , ⌽t . dwt ,
Ž 1.
where ¨ t is an s-dimensional actuating signal and t is the n-dimensional plant-state. In this paper, we will ignore endogenous actuation: ¨ t ' 0. The values of the regime process are associated with stasis condition for the plant. The system operates near one set-point for an unpredictable interval, then suddenly shifts to another. The set of nominal plant-states can be arrayed as an Ž n = S .-matrix , with the nominal associated with ith mode given by i : s w i x. During operation, the system will sojourn in one regime for a time Ž ⌽t s ⌽i for t g w a, b .., and then suddenly shift Ž ⌽ b s ⌽j . to another in response to an external event or change in the surrounding environment. Despite the simple form of the regime process, the equation of plant evolution is intractable. As an approximation, Ž1. is often replaced with a family of local models created by linearizing the plant dynamics about the various set-points. Let t be a pointer to the current regime: The state space of t consists of the S canonical unit vectors in ⺢ s Ž t g e 1 , . . . , e S 4.. The component in t with value one marks the current mode of operation: If ⌽t s ⌽ k , then t s e k . The modal dynamics are usually represented by an exogenous Markov model with generator Q⬘. The base-state is the deviation of the plant-state from the current set-point: x t s t y t . The t 4 process is called the modal-state process to differentiate it from the base-state process. The hybrid-state of the system is the composite of the base- and modal-states. A comprehensive state model requires a representation of evolution, both intramodal and intermodal. During an extended modal sojourn, proper control will place and maintain the plant-state vector near the 1 For notational convenience in what follows: S will designate an integer index set 1, . . . , S4; e i is the ith canonical unit vector in a space whose dimension is obvious from the context; where no confusion will arise, a subscript may identify time or the component of a vector Že.g., 1 4 is a scalar process that is the first component of the vector process t 4.; ) is the Hadamard product ŽŽ x ) y . i s x i yi .; and m indicates Kronecker product. If a process is sampled T times per second, the discrete sequence so generated is written yw k x4 where the index denotes sample number rather than time. A process discontinuity is labeled ⌬ x t s x t y Žtime-continuous. or ⌬ x w k x s x w k x y x w k y 1x Žtime-discrete.. Conditional expectation is xy t denoted with a circumflex with the relevant -field apparent from context. A Gaussian random variable with mean ˆ x t and covariance Px x is labeled x g NŽ ˆ x t , Px x . with the same symbol used for the density function itself where no confusion will arise. Denote the density function of a unit Gaussian random vector by ⌽ Ž u.. If A is a positive symmetric matrix and x is a compatible vector, x⬘ Ax is denoted 5 x 5 2A .
MODAL ESTIMATION IN HYBRID SYSTEMS
227
correct set-point. The natural plant model in this circumstance would be that local model, selected from the family of regime-specific, linear models, associated with the current set point. But there will be a base-state discontinuity when the mode changes. If e i ¬ e l , the base-state reference level changes from i to l : ⌬ x t s y ⌬ t .
Ž 2.
Because of the wide applicability of hybrid systems, considerable effort has been devoted to problems of estimation. Both time-continuous and time-discrete models have been used, the latter formed by sampling the relevant variables of the former at rate 1rT. The time-discrete hybrid model accepts certain temporal coincidences; e.g., the modal state is constant between sample times, but if T is small enough, the resulting errors are small. The time-discrete hybrid model can be developed as is done in w11x, xw k x s
Ý Ž A j x w k y 1x q C j w w k x . eXj w k y 1x y ⌬ w k x
Ž 3.
j
w k x s ⌸ w k y 1x q mw k x ,
Ž 4.
where ww k x4 is a unit Gaussian white sequence. If the local models are controllable from the plant noise in every regime, for i g S, Ci and A i are both nonsingular. Label the local plant-noise covariance with R Ž i . s D Ž i .y1 s Ci CiX ) 0. The time-discrete modal transition probabilities are given by ⌸. In the intersample interval wŽ k y 1.T, kT ., the base-state evolves according to the mode at the beginning of the interval. If there is a modal transition, the base-state will experience the discontinuity y ⌬ w k x. Particularly when the value of t 4 is important, a measurement architecture in which the plant-state and the modal-state are assigned separate sensors has been shown to be effective. The sensors could be time-continuous or time-discrete, but the focus here will be on the latter. Suppose the plant-state sensor gives a linear, albeit mode dependent, measurement of the plant-state vector. Alternatively, the modal sensor is a classifier with quality represented by the S = S discernibility matrix D s wDi j x where Di j is the probability that model-state e i will be selected by the classifier if e j is the true mode at time of measurement. The measurement model is then given by yw k x s
Ý Hi eXi w k x w k x q nw k x
Ž 5.
i
z w k x s D w k x q w k x ,
Ž 6.
228
SWORDER, BOYD, AND ELLIOTT
where nw k x4 is a Gaussian white sequence with positive covariance R x s Dy1 s F⬘F. The columns of D are probability vectors, and the form is x flexible enough to include aliasing and bias. The measurement rate is determined by the sensor and its utilization policy. The measurements may occur at a different rate than the basic clock rateᎏusually much slower. Should the observation rate be slower, the sensor gains can be made time dependent so that at times of no measurement yw k x and z w k x are uninformative. This will be done in what follows without comment. The initial plant state categories are assumed to be independent with probability distributions x w0x ; NŽ ˆ x w0x, Px x w0x. and w0x ; ˆw0x. Denote the filtration generated by ww k x4 and mw k x4 along with the initial conditions on the plant by O w k x. Then ww k x4 and mw k x4 are O w k y 1x-martingale increment sequences. Denote the filtration generated by the plant measurements yw k x4 by Y w k x Žrespectively, that generated by the modal observation by Z w k x.. The observation filtration is G w k x s Y w k x k Z w k x. Let F w k y 1x s G w k y 1x k O w k x. The exogenous processes in Ž5. and Ž6. are F w k y 1x-martingale increments. Other authors study state estimation and regulation with an emphasis on time-continuous hybrid plants and a mix of time-continuous and timediscrete observations. Often, an engineer seeks practical algorithms for approximating the Gt-regime probabilities along with those Gt-moments Žincluding cross moments. important in the application. The polymorphic estimator Žthe PME. is a useful tool w10, 11, 14x. The PME is moment based and does not provide the Gt-distribution function of the hybrid-state. Further, an asymmetry in the data processing reduces the quality of the modal estimate. In this paper, we will provide an explicit processing algorithm for approximating the Gt-probability distribution of the plantstate t and the regime index ⌽t . This algorithm affords a more tightly integrated fusion of Y w k x and Z w k x than is provided by the PME, and it leads to improved identification of the regime process.
2. THE BAYES RECURRENCE FORMULA In contrast to the moment-based approach used to frame the PME, we will seek an implementable approximation to the G w k x-distribution of the hybrid-state. At time t s kT, suppose the G w k x-density is known: pw k x s w pi w k xx where pi w k x Ž z . dz s P Ž x w k x g w z, z q dz x , w k x s e i < G w k x . .
Ž 7.
We seek a mapping from Ž pw k x, yw k q 1x, z w k q 1x. to pw k q 1x. This mapping can be displayed most concisely when phrased in terms of an
229
MODAL ESTIMATION IN HYBRID SYSTEMS
unnormalized density, qw k xŽ z . s w qi w k xŽ z .x, derived using a reference probability. On the original event space and filtration, Ž ⍀, F ; F w k x., consider a reference probability measure P with respect to which mw k x4 and ww k x4 retain their earlier character, but yw k x4 is a unit Gaussian white sequence and z w k x4 is an independent, identically distributed Žiid. sequence that is uniformly distributed across e 1 , . . . , e S 4 . For l g 0, 1, . . . 4 , define
w l x s
S ⌽Ž yw l x.
z w l x ⬘ D w l x < F
= Fy1 y w l x y
ž ž
Ý Hi eXi w l x Ž e w l x q i . i
//
Ž 8.
and let ⌳w k x be a continuing product of the w l x: ⌳w k x s Ł 0k w l x. It is shown in w11x that P and P are related by
⭸P ⭸P
F w k xk G w k x
s ⌳w k x Ž zw k x , yw k x. .
Ž 9.
This change to P reflects no change for the plant disturbances, but the character of the observation is considerably different. Denote expectation with respect to P by E. The conditional Bayes theorem w3, Theorem 3.2, Chap. 2x relates the operator E to the operator E. Let f be a scalar test function. The expectation of eXi ⌽w k q 1x f Ž x w k q 1x. can be computed using either P or P E eXi w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x E ⌳ w k q 1x < G w k q 1x
.
Ž 10 .
Equation Ž10. shows that an expectation with respect to P can be performed as an expectation with respect to the reference measure P if the derivative of the measure is included properly in the calculation. The denominator has nothing to do with f and is simply a normalizing factor. The numerator will be thought of as the G w k q 1x-conditional expectation of f Ž x w k q 1x. with respect to an unnormalized G w k q 1x-conditional probability density, qi w k q 1xŽ x w k q 1x.: E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
H⍀ f Ž z . q w k q 1x Ž z . dz. i
Ž 11 .
230
SWORDER, BOYD, AND ELLIOTT
If qi w k q 1x; i g S4 were known, the joint density = mass function of the hybrid-state would be pi w k q 1 x Ž z . s
qi w k q 1 x Ž z . Ý j H⍀ q j w k q 1 x Ž u . du
,
Ž 12 .
from which ˆ x w k q 1x and ˆw k q 1x along with higher moments could be calculated. To find qi w k q 1xŽ z .; i g S4 , expand Ž11. by replacing w k q 1x, ⌳w k q 1x, and x w k q 1x with their values. E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s E eXi Ž ⌸ w k x q m w k q 1 x . ⌳ w k x z w k q 1 x ⬘De i < F
žÝŽ
S ⌽ Ž y w k q 1x .
A j x w k x q C j w w k q 1 x . eXj w k x y Ž e i y w k x . < G w k q 1 x .
/
j
Ž 13 . It is true that Ý j w k x⬘e j ' 1. Substituting this into Ž13., E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x sE
Ý Ž w k x ⬘e j . ŽeXi ⌸ e j . ⌳ w k x z w k q 1x ⬘De i < F
⭈⌽ Ž Fy1 Ž y w k q 1 x y Hi Ž x w k q 1 x q e i . . . ⭈
S ⌽ Ž y w k q 1x .
f Ž A j x w k x q C j w w k q 1x y Ž e i y e j . . < G w k q 1x .
Under P, G w k x4 is uninformative with respect to the hybrid-state. Since an unnormalized distribution is sought, factors common to all regimes can and will be ignored in what follows. E eXi w k q 1 x ⌳ w k q 1 x f Ž x w k q 1 x . < G w k q 1 x s
H⍀ Ý ⌸
ij z
w k q 1 x ⬘D .i ⌽
j
⭈ Fy1 Ž y w k q 1 x y Hi Ž A j q C j w y Ž e i y e j . q e i . .
ž
⭈ f Ž A j q C j w y Ž e i y e j . . q j w k x Ž . ⌽ Ž w . d dw.
/
MODAL ESTIMATION IN HYBRID SYSTEMS
231
To simplify this, make the change of variable z s A j q C j w y i q j . Then d dw s < C j
H⍀ Ý ⌸
ij z
w k q 1 x ⬘D .i < C j
j
⭈ f Ž z . q j w k x Ž . ⌽ Ž Cy1 j Ž z y A j q i y j . . d dz s
H⍀ f Ž z . q w k q 1x Ž z . dz. i
Hence Observation update qi w k q 1xŽ z . s Ý ⌸ i j z w k q 1x⬘D.i < C j
⭈
H⍀⌽ Ž C
y1 j
Ž z y A j q i y j . . q j w k xŽ . d .
Ž14.
Equation Ž14. is the recurrence formula for the unnormalized density of the hybrid-state. Variants on Ž14. have been derived for kindred plant representations: If s 0 and the noise in z w k x4 is Gaussian, see w4, Eq. Ž8.x. If s 0 and there is a single measurement g w k x s yw k x q z w k x in Gaussian noise, see w5, Eq. Ž10.x. If s 0 and there is a slight indexing difference in the base-state model, see w5, Eq. Ž4.x. 䢇 䢇
䢇
Unfortunately, Ž14. is not an algorithm of the form we seekᎏalthough recursive, it is infinite dimensional.
3. A GAUSSIAN WAVELET ESTIMATOR We seek a finite dimensional, Gaussian wavelet approximation to qw k x. The number of terms is limited only by the computational complexity permitted in the application. We will assume that N terms will suffice in every regime. The resulting algorithm will be called the Gaussian wavelet estimator ŽGWE.. For notational convenience, let the difference in the ith and the jth set-points be written ji s i y j , and let the observation
232
SWORDER, BOYD, AND ELLIOTT
centered about the ith set-point be written y i w k q 1x s yw k q 1x y Hi i . Denote the inverse of a Žpositive. covariance matrix P by D, e.g., D li w k x s Pli w k xy1 , and the product Dm by d, e.g., d li w k x s D li w k x m il w k x. Suppose qj w k x Ž . s
N
Ý ␣ lj w k x
D lj w k x
1r2
1
exp y
2
ls1
y m lj w k x
2 D lj w k x .
Ž 15 .
In Ž15., the unnormalized G w k x-conditional probability of the event x w k x g w , q d x, w k x s e j 4 is given under P by a sum of N Gaussian pattern functions, N Ž m lj w k x, Pl j w k x.; l g N4 with means m lj w k x and positive covariances Pl j w k x. In this approximation, m lj w k x translates the lth pattern function, Pl j w k x adjusts its shape, and N circumscribes the span of the sum. The pattern functions are weighted by ␣ lj w k x; l g N4 . All of the coefficients are G w k x-adapted. There are NS elements in qw k x though many could be zero. The Appendix presents a recurrence formula for the coefficients ␣ lj w k x4 , m lj w k x4 , Pl j w k x4 Ž j g S, l g N.. The recurrence can be most concisely stated in a mixed covariance-information form also used in the Kalman filter w1x. The information matrix is the inverse of the covariance matrix, and instead of framing the estimator using m il 4 as coordinates, the information filter propagates d li w k x. In this mixed coordinate system, the base-state estimate satisfies base-state recurrence: GWE Extrapolation: m il w k q 1xys A i m il w k x; Pli w k q 1xys
A i Pli w k x AXi
i g S,
lgN
q R Ž i .;
i g S,
Ž16. lgN
Ž17.
Update: d i w k q 1xŽ l; j .qs d lj w k q 1xyq HiX Dx y i w k q 1x yD lj w k q 1xy ji ;
i, j g S,
D i w k q 1xŽ l; j .qs D lj w k q 1xyq HiX Dx Hi ;
Ž18. lgN
i, j g S,
l g N.
Ž19.
As in many multiple model algorithms, Ž16. ᎐ Ž19. require a filter bank of size S for extrapolation and update. But in the GWE, each filter extrapolates N initial conditions. After adjusting for the inclusion of the set points, the form of the GWE is close to that of the Kalman filter.
233
MODAL ESTIMATION IN HYBRID SYSTEMS
Under the nontransition hypothesis, Že i ¬ e i ., ⌬ d li w k q 1xqs HiX Dx y i w k q 1x, the Kalman update using the centered observation. Under the transition hypothesis Že j ¬ e i ., the increment in d i w k q 1xŽ l; j .q contains y i w k q 1x as a factor but the translation of the set-point, Žy i q j ., must be included. The weightings in d i w k q 1xŽ l; j .q balance measurement fidelity and the shape of the pattern functions. The size of jump q
d i w k q 1 x Ž l ; j . y d lj w k q 1 x
y
in the Ž l; j .th subfilter will be increased if the measurement noise is small and if yw k q 1x differs from Hi i . Similarly, the set-point offset Ž j y i . is important when the probability of the lth pattern function is concentrated in a small region, i.e., when D lj w k q 1xy is big. Suppose the modal measurement is quite good and an e j ¬ e i transition occurs at time t s Ž k q 1.T. In the absence of the base-state measurement, the extrapolation in the base-state mean should be q
y
m i w k q 1 x Ž l ; j . f m lj w k q 1 x q j y i . The mean translates by the amount of the base-state discontinuity. The size of the difference, m i w k q 1xŽ l; j .qy m lj w k q 1xyq ji, is a coarse measure of the influence of yw k q 1x since this difference is that between the corrected and uncorrected extrapolation. Modal mixing is distributed in the GWE. It appears in the filter update and also in the weighting coefficients, ␣ i w k q 1xŽ l; j .q; i, j g S, l g N4 . modal mixing: GWE
␣ i w k q 1xŽ l; j .qs ␣ lj w k x⌸ i j z w k q 1x⬘D.i ⭈ D lj w k x Ž D lj w k x q AXj D Ž j . A j . ⭈exp
1 2
½
m i w k q 1 x Ž l; j . y
y1
q
P i w k q 1 x Ž l; j . D Ž j .
1r2
q 2 q D i w kq 1 x Ž l; j .
y m lj w k q 1 x y ji
2 y w kq 1 x y D lj
y i w k q 1x
2 Dx
5.
Ž20.
The weighting coefficient ␣ i w k q 1xŽ l; j .q is proportional to its predecessor Ž ␣ lj w k x., proportional to the probability of an e j ¬ e i transition Ž ⌸ i j ., and proportional to the probability that the modal observation was generated from w k q 1x s e i Žthe factor z w k q 1x⬘D.i ..
234
SWORDER, BOYD, AND ELLIOTT
With coefficients thus determined, the unnormalized G w k q 1x-wavelet distribution of the hybrid-state can be written immediately: q
q
q
qi w k q 1 x s Ý ␣ i w k q 1 x Ž l ; j . N Ž m i w k q 1 x Ž l ; j . , P i w k q 1 x Ž l ; j .
q
..
j, l
Ž 21 . From Ž21., the various moments of interest can be computed; e.g.,
ˆi w k q 1 x s
Ý ␣ i w k q 1x Ž l ; j .
q
j, l
Ý ␣ i w k q 1x Ž l ; j .
q
Ž 22 .
i , j, l
and q
ˆx w k q 1 x s Ý ␣ i w k q 1 x Ž l ; j . m i w k q 1 x Ž l ; j . i , j, l
q
Ý ␣ i w k q 1x Ž l ; j .
q
.
i , j, l
Ž 23 . Unfortunately, there are not N terms in Ž21. but NS. To satisfy the complexity constraint, pruning andror merging must be used to reduce the number of terms. For now we will reduce the terms from each progenitor to one by using the conventional Gaussian sum merging formula given in w1x. For a given i g S, consider the jth progenitor. Define
␣ ji w k q 1 x Ž l . s ␣ i w k q 1 x Ž l ; j .
q
Ý ␣ i w k q 1x Ž l ; j .
q
.
Ž 24 .
l
Then let m ij w k q 1 x s
Ý mi w k q 1x Ž l ; j .
q
␣ ji w k q 1 x Ž l .
Ž 25 .
l
Pji w k q 1 x s
Ý ž P i w k q 1x Ž l ; j .
q
l
q
q m i w k q 1 x Ž l ; j . y m ij w k q 1 x
ž
q
/
= m i w k q 1 x Ž l ; j . y m ij w k q 1 x ⬘ ␣ ji w k q 1 x Ž l . Ž 26 .
ž
␣ ji w k q 1 x s
q Ý ␣ i w k q 1x Ž l ; j . .
l
//
Ž 27 .
235
MODAL ESTIMATION IN HYBRID SYSTEMS
The next iteration of the GWE begins with a set of S terms: N s S. updated distribution qi w k q 1x s
Ý ␣ jiw k q 1xNŽ mij w k q 1x, Pjiw k q 1x..
Ž28.
j
In many applications, there is a common set-point and, though the system dynamics change, the reference level does not. Such systems are called linear-jump-systems Žsee w8x. and are modeled with s 0. The GWE-algorithm is considerably simplified in this case. base-state recurrence: GWE Ž s 0. Extrapolation: m il w k q 1xys A i m il w k x;
i, l g S
Ž29.
Pli w k q 1xys A i Pli w k x AXi q R Ž i .;
i, l g S
Ž30.
Update: ⌬ d i w k q 1xŽ l; j .qs HiX Dx yw k q 1x; ⌬ D i w k q 1xŽ l; j .qs HiX Dx Hi ;
i, j, l g S
Ž31.
i, j, l g S,
Ž32.
where ⌬ d i w k q 1xŽ l; j .qs d i w k q 1xŽ l; j .qy d lj w k q 1xy, and similarly for ⌬ D i w k q 1xŽ l; j .q. The weighting coefficients can be simplified too. For i, j, l g S,
␣ i w k q 1xŽ l; j .qs ␣ lj w k x⌸ i j z w k q 1x⬘D.i < D lj w k xŽ D lj w k x q AXj D Ž j . A j .y1 ⭈P i w k q 1xŽ l; j .q < 1r2 exp 12 ⌬ m i w k q 1xŽ l; j .q 5 2D i w kq1xŽ l; j.q.
4. MULTIPLE MODEL ALGORITHMS The simplest of the GWEs are those for systems with a single set point and N s 1. Then qj w k x Ž . s ␣ j w k x D j w k x
1r2
exp y
1 2
y m jw k x
2 D jw k x .
Ž 33 .
236
SWORDER, BOYD, AND ELLIOTT
The extrapolation equations of GWE Ž s 0, N s 1. are direct restrictions of Ž29. ᎐ Ž32.. The extrapolation uses S parallel filters with some mixing in the update. base-state recurrence: GWE Ž s 0. Extrapolation: m i w k q 1xys A i m i w k x; igS p i w k q 1xys A i P i w k x AXi q R Ž i .;
Ž34. Ž35.
igS
Update: ⌬ d i w k q 1xŽ j .qs HiX Dx yw k q 1x; i, j g S ⌬ D i w k q 1xŽ j .qs HiX Dx Hi ; i, j g S
Ž36. Ž37.
The weighting coefficients can be written q
␣ i w k q 1 x Ž j . s ␣ j w k x ⌸ i j z w k q 1 x ⬘D .i D j w k x Ž D j w k x q AXj D Ž j . A j . ⭈P i w k q 1 x Ž j .
q
1r2
exp 12 ⌬ m i w k q 1 x Ž j .
q 2 D i w kq1 xŽ j .q .
y1
Ž 38 .
The unnormalized modal distribution is
␣ i w k q 1 x s z w k q 1 x ⬘D .i
Ý ␣ i w k q 1x Ž j .
q
.
Ž 39 .
j
For systems with a single set-point, several estimation agorithms have been proposed in the literature. That receiving most attention is the interacting-multiple-model ŽIMM. estimator w2x. As usually implemented, the IMM generates a Gaussian sum approximation to the Y w k x-density of the hybrid-state in the form Ý j eXj ˆw k xNŽ m ˆ j w k q 1x, P j w k q 1x.. The equations of implementation for the IMM are provided in several places Že.g., see w7, Table 1x., and an annotated description of applications of the IMM is found in w9x. The IMM was modified to accept the z w k x4 observations in w6x, and this G w k x-IMM was called the image-enhanced IMM. The steps in the implementation of the IMM are given explicitly in w6, Fig. 3x. This implementation of the IMM is closest to the GWE and will be referred to simply as the IMM. The motivation behind the IMM is the efficient computation of a high quality approximation to ˆ x w k x4 . It is not claimed that the IMM produces the G w k x-joint density of the hybrid-state, nor does it. Approximations created by the merging of the subfilters and the update rule preclude
MODAL ESTIMATION IN HYBRID SYSTEMS
237
accurate calculation of ˆw k x4 . Tests have shown that a high quality estimate of the modal-state is not required for accurate base-state estimation because the IMM is essentially a variable bandwidth filter for the latter. The IMM and the GWE Ž s 0, N s 1. perform common operations on the data. The IMM does hypothesis mixing at the beginning of an iteration while the GWE Ž s 0, N s 1. does it at then end. The factor z w k q 1x⬘D.i appears in both along with factors akin to Ý j ␣ j w k x⌸ i j Žafter adjusting for the time of hypothesis merging.. Both require S Kalmanextrapolations. Moreover, the IMM has S Kalman-updates while the GWE Ž s 0, N s 1. has S 2 . There are some differences in indexing because the IMM supposes that the evolution of x t 4 on wŽ K y 1.T, KT . depends on w k x while the GWE supposes a w k y 1x dependence. For small T, this matters little. The fundamental difference between the algorithms arises from the modal update. The ith filter for either the IMM or the GWE provides m i w k q 1xy as shown in Ž34.. Define the ith residual by r i w k q 1x s yw k q 1x y Hi m i w k q 1xy For a system with a linear-Gauss ᎐Markov ŽLGM. model Ž S s 1., r i w k x is the increment in the Y w k x-innovation processᎏ r i w k x4 is a Y w k x-martingale increment sequence. For the LGM estimation problem, the covariance of r i w k q 1x is y
Pyi y w k q 1 x s Hi P i w k q 1 x HiX q R x Ž i . .
Ž 40 .
The IMM uses multiple models, and the update of the modal probabilities depends upon the behavior of the subfilter residuals. The smaller the residuals in a particular subfilter are, weighted by the associated information matrix, the more likely the associated mode is thought to be. The unnormalized modal update takes the form y
q w k x s ˆ w k x ) Dyi y w k q 1 x
y 1r2
=exp y 12 r i w k q 1 x ⬘Dyi y w k q 1 x r i w k q 1 x 4 . Ž 41 . In contrast, the elements from which the GWE generates the modal estimate Žsee Ž38.. are responsive, not to the residuals themselves, but rather to the effect these residuals have on the update of the mean. The weight ␣ i w k q 1xŽ l; j .q is generated as a continuing exponential product of the quadratic differences of the base-state mean before and after a base-state measurement: ˆw k x depends not on the residuals themselves, but on the influence the residuals have on the base-state estimate.
238
SWORDER, BOYD, AND ELLIOTT
The updates of the IMM and GWE Ž s 0, N s 1. are closely related in certain cases. Denote the update gain of the ith filter in the IMM by W i and suppose the W i 4 are all nonsingularᎏthis unfortunately requires that the Hi 4 are all nonsingular. Then r i w k q 1x s ŽW i .y1 ⌬ m i w k q 1x. Expanding the quadratic form in Ž41. r i w k q 1x
2 D iy y w kq1 x
y
i s ⌬ m i w k q 1 x ⬘D i w k q 1 x Hy1 i Py y w k q 1 x
⭈ HyT D i w k q 1x i
y
⌬ mi w k q 1x y
s ⌬ m i w k q 1 x ⬘D i w k q 1 x Hy1 i y
⭈ Hi P i w k q 1 x HiX q R Ž i .
ž
⭈ HyT D i w k q 1x i
y
/
⌬ mi w k q 1x .
If R Ž i . is small Žas occurs during intermodal transients . and if the base-state observation gain is nonsingular Žan unusual event. r i w k q 1x
2 D iy y w kq1 x
f ⌬ mi w k q 1x
2 D i w kq1 xy
/ ⌬ mi w k q 1x Ž j .
q 2 D i w kq1 xŽ j .q .
Ž 42 .
The IMM and GWE Ž s 0, N s 1. have similar modal updates, but the IMM weights a quadratic form in the increments of the estimate while the GWE weights an increment in the quadratic forms. The update of the conditional mean for GWE Ž s 0, N s 1. involves more bookkeeping than does the similar operation IMM. Modal mixing takes place even in the absence of set-point discontinuities.
5. THE GWE IN ESTIMATION AND IDENTIFICATION To illustrate the performance of the GWE, consider a specific example. On the California desert stands a 10 MWe solar electric generating system. Movable mirrors Žheliostats . are used to focus the sun’s energy on a group of boiler panels atop a central tower. A steam temperature regulator controls the feedwater flow rate to the panels so as to maintain the proper temperature of both the outlet steam and receiver panel itself w13x. A remote observer wishes to estimate the panel temperature and identify the current set-point. On a partly cloudy day the insolation changes suddenly and unpredictably. The panel dynamics change as well: at low insolation, the panel dynamics are slow; at high insolation, the panel dynamics are fast. Let Tm be the panel temperature and let Tc be the
239
MODAL ESTIMATION IN HYBRID SYSTEMS
output of a temperature sensor, with m Žrespectively c . the perturbation variables. On a day where panel operation moves between sun and cloud, the local models of the panel and the sensor are given by Panel model dm s ␣ m dt q  dw1 ,
Ž 43 .
dTc s a Ž Tc y Tm . dt q 0.3  dw 2 .
Ž 44 .
Sensor model
The coefficients of the local panel models are given as a function of insolation in w12x and are displayed in Table I. The nominal panel temperature at high insolation is 17⬚F higher than that for low insolation: s Ž1017, 1000. m 1. The state ordering in what follows is x s vecŽm , c .. To particularize the modal dynamics, Q must be specified. The mean sojourn of bright sunlight is 2 min while clouds persist for an average of 4 min. A simple Markov model matches these mean sojourns. The sensors at the site are of two types: a direct measurement of Tc in noise and an insolation measurement. For the purposes of this example, the former will be sampled every 0.5 min in white Gaussian noise of standard deviation 2⬚F: yw k x s D w k x q nw k x, where D s eX2 and R x s 4. The panel measurement is of absolute temperature, not the perturbation variable, c . To improve performance, the insolation sensor gives a direct measurement of 4 . In the application, the insolation is not uniformly distributed across the panel: the geometry of the sun-mirror-panel is such that the distribution of the solar flux across a panel is irregular with possible hot spots. Insolation sensors measure conditions only in a neighborhood. The application employed several insolation sensors arrayed across the panel, and the highest reading was used as a measure of panel insolation. The insolation sensor is sampled every 0.1 min. In w12x several panel temperature estimators were contrasted. An EKF using pseudonoise to represent the modal transitions was a failure. When the discernibility matrix indi-
TABLE I Local Panel Temperature Dynamics for Two Insolation Conditions Modal-state
Insolation
¨n
n
␣

e1
1.8 = 10 BTUrmin 186 lbrmin
1017⬚F y1.8
e2
2.7 = 10 BTUrmin
1000⬚F y0.36 1.9
5
4
30 lbrmin
a
1.58 0.9 0.25
240
SWORDER, BOYD, AND ELLIOTT
cated high quality modal measurements, PME was effective for both panel temperature estimation and modal identification. The performance of the GWE in this circumstance is comparable. A lower quality D created difficulties for the PME because the hypotheses underlying the algorithm became more difficult to justify. Suppose the less accurate insolation sensor introduced in w12x is used in this application: D1 s
0.65 0.35
0.40 . 0.60
Ž 45 .
The inclination toward high insolation readings is accommodated with this asymmetrical D. The precision of the classifier is such that the modal-state is difficult to identify with confidence. Figure 1 shows sample functions of Tm and Tc for an insolation path: has the value e1 on t g w0, 1. j w6, 9. and the value e2 on w1, 6. j w9, 11x. Both temperatures were initialized at 1017⬚F. After the cloud comes by at t s 1 min, Tm moves rapidly toward 1000⬚F, with Tc moving at a much slower rate. A 6⬚F temperature difference between the sensor and the panel develops in 5 min. When the sun returns, the sensor output has significant lag. A sample observation sequence, yw k x4 , is also shown and is marked on the figure with an ‘‘>.’’
FIG. 1. A sample path of the panel temperature Tm and its measurement Tc . Noisy temperature measurements are shown every 0.5 min.
MODAL ESTIMATION IN HYBRID SYSTEMS
241
An estimator that attempted to determine t 4 from yw k x4 alone would have a difficult time. Even determining Tm from such an abbreviated data set is a difficult task indeed. Figure 2 shows a plot of 1 4 Žsojourns in the sun., the insolation measurements Žlabeled >., and probability of sun as computed by GWE. The poor quality of the insolation sensor is reflected in a slow response to modal transition and a generally low confidence in the modal identification. Even so the modal identification is superior to that achieved by the PME in w12x. The performance of the panel temperature estimator is shown in Fig. 3. True temperature is shown along with the observations. The GWE initializes itself in the first minute. The GWE follows the metal temperature, but with some lag. The GWE is superior to the PME when the insolation sensor is of such low quality. Slight improvement in the modal measurements yields much improved set-point identification. The discernibility matrix D 2 is a higher quality insolation sensor, and the modal estimates shown in Fig. 4 reflect this:
D2 s
0.75 0.25
0.30 . 0.70
Ž 46 .
FIG. 2. The probability of sun as computed by GWE along with the mode. The insolation measurements are denoted with an >.
242
SWORDER, BOYD, AND ELLIOTT
FIG. 3. Estimate of the panel temperature using the GWE and sensor D1 .
FIG. 4. The probability of sun as computed by GWE along with the true mode. The insolation measurements are denoted with an >.
243
MODAL ESTIMATION IN HYBRID SYSTEMS
6. CONCLUSION The GWE is a novel, finite dimensional, approximation to the G w k x density of the hybrid-state. It provides the solution in the form of a Gaussian sum. In contrast to the PME, the observation sequences yw k x4 and z w k x4 are fused in a symmetric manner. When the set-point reference is constant, the GWE is considerably simpler. Indeed, GWE Ž s 0, N s 1. is about as simple as the IMM. The GWE will accommodate variable set-points, a subtle problem illustrated in the example.
APPENDIX At time t s kT, qw k x is a Gaussian sum with N terms in each component: N
qj w k x Ž . s
Ý ␣ lj w k x
D lj w k x
1r2
exp y
y m lj w k x
1 2
2 D lj w k x .
ls1
The recurrence relation for qi w k q 1x is given in Ž14.. qi w k q 1 x Ž z . s
Ý ⌸ i j z w k q 1x ⬘D.i < C j
H⍀ ⌽ ž C Ž z y A q . / q w k x Ž . d .
⭈
y1 j
j
i j
j
Expand the product of Gaussian pattern functions. In what follows, factors common across all modal hypotheses will be ignored without comment; i ⌽ Ž Fy1 Ž y i w k q 1 x y Hi z . . ⌽ Cy1 j Ž z y A j q j . q j w k x Ž .
ž
N
s
Ý ␣ lj w k x
D lj w k x
1r2
/
exp y 12 J1 Ž l . ,
ls1
where J1 Ž l . s y m lj w k x
2 D lj w k x
i q Cy1 j Ž z y A j q j .
q Fy1 Ž y i w k q 1 x y Hi z .
2
.
2
244
SWORDER, BOYD, AND ELLIOTT
The argument of the exponential is a quadratic form. It is an exercise to complete the squares in and z. Define a set of coefficient matrices: F j w k x Ž l ; i . s D Ž j . q HiX Dx Hi H j s AXj D Ž j . A lj w k x s D lj w k x q AXj D Ž j . A j E j w k x Ž l ; i . s D lj w k x m lj w k x q H j ji K j w k x Ž l ; i . s HiX Dx y i w k q 1 x y D Ž j . ji G j w k x Ž l ; i . s m lj w k x
2 D lj w k x
q ji
2
q y i w k q 1x
D Ž j .
2 Dx .
Completing the squares and combining exp y 12 J1 Ž l . s A lj w k x
y Ž1r2 .
⭈ exp y
1 2
D i w k q 1x Ž l ; j .
½G w k x Ž l ; i . y j
y Ž1r2 .
E jw k xŽ l; i.
y d i w k q 1x Ž l ; j . ⭈ N A lj w k x
ž
y1
2 A lj w k xy 1
2 P i w kq1 xŽ l ; j .
Ž E j w k x Ž l ; i . y Hzj . , A lj w k x y1 /
⭈ N z Ž mi w k q 1x Ž l ; j . , P i w k q 1x Ž l ; j . . , where D i w k q 1 x Ž l ; j . s F j w k x Ž l ; i . y H j⬘A lj w k x s P i w k q 1x Ž l ; j .
y1
y1
Hj
and
mi w k q 1x Ž l ; j . s P i w k q 1 x Ž l ; j . Ž H j⬘A lj w k x
y1
E jw k xŽ l; i. q K jw k xŽ l; i. . .
Now integrate
H⍀ ⌽ Ž F Ž y w k q 1x y H z . . ⌽ ž C Ž z y A q . / q w k x Ž . d y1
N
s
i
Ý ␣ lj w k x
i
D lj w k x
1r2
y1 j
A lj w k x
j
y Ž1r2 .
ls1
⭈ F j w k x Ž l ; i . y Ž H j . ⬘ A lj w k x
y1
Hj
y Ž1r2 .
i j
j
245
MODAL ESTIMATION IN HYBRID SYSTEMS
⭈ exp y
1 2
½G w k x Ž l ; i . y j
E jw k xŽ l; i.
y d i w k q 1x Ž l ; j .
2 A lj w k xy 1
2 P i w kq1 xŽ l ; j .
5
⭈ N z Ž mi w k q 1x Ž l ; j . , P i w k q 1x Ž l ; j . . . Then qi w k q 1x is a Gaussian sum with 1r2
␣ i w k q 1 x Ž l ; j . s ⌸ i j z w k q 1 x ⬘D .i < C j
y Ž1r2 .
1 2
D i w k q 1x Ž l ; j .
½G w k x Ž l ; i . y j
y Ž1r2 .
E jw k xŽ l; i.
y mi w k q 1x Ž l ; j .
2 A lj w k xy 1
2 D i w kq1 xŽ l ; j .
5.
The expressions for the coefficients can be simplified. Note first that D i w k q 1 x Ž l ; j . s D Ž j . q HiX Dx Hi y D Ž j . A j Ž D lj w k x q AXj D Ž j . A j . s HiX Dx Hi q Ž R Ž j . q A j Pl j w k x AXj .
y1
y1
AXj D Ž j .
.
This can be written in more conventional form as y
Pl j w k q 1 x s A j Pl j w k x AXj q R Ž j . followed by y
D i w k q 1 x Ž l ; j . s D lj w k q 1 x q HiX Dx Hi . Also d i w k q 1x Ž l ; j . s D Ž j . A j Ž D lj w k x q AXj D Ž j . A j .
y1
d lj w k x
q D Ž j . A j Ž D lj w k x q AXj D Ž j . A j .
ž
q
HiX Dx
y1
AXj D Ž j . y D Ž j . ji
y i w k q 1x .
Write this as y
d lj w k q 1 x s D Ž j . A j Ž D lj w k x q AXj D Ž j . A j . y
y1
d lj w k x y
d i w k q 1 x Ž l ; j . s d lj w k q 1 x q HiX Dx y i w k q 1 x y D lj w k q 1 x ji .
/
246
SWORDER, BOYD, AND ELLIOTT
Note that
Ž Dlj w k x q AXj D Ž j . A j .
y1
yT s AyT I y Ž Ay1 q Pl j w k x . j j R Ž j . A j
ž
y1
yT Ay1 . j R Ž j . A j
/
So y
Pl j w k q 1 x D Ž j . A Ž D lj w k x q AXj D Ž j . A j . yT y Ž Ay1 q Pl j w k x . j R Ž j . A j
y1
y1
R Ž j . AyT j
yT yT s A j Ž Ay1 q Pl j w k x . y A j Ay1 j R Ž j . A j j R Ž j . A j
s A j Pl j w k x . Hence y
y
Pl j w k q 1 x d lj w k q 1 x s A j Pl j w k x d lj w k x s A j m lj w k x . This can be written y
m lj w k q 1 x s A j m lj w k x . The expression for the weighting coefficients is complicated in appearance. The argument in the exponential factor in ␣ i w k q 1xŽ l; j . can be rewritten: d lj w k x
2 P lj w k x
q 5 ji 5 2D Ž j. q y i w k q 1 x
y d lj w k x q AXj D Ž j . ji
2 A lj w k xy1
2 Dx
y d i w k q 1x Ž l ; j .
2 P i w kq1 xŽ l ; j . .
But d lj w k x
2 P lj w k x
y d lj w k x
2 A lj w k xy 1 y
s d lj w k x ⬘Pl j w k x AXj Pl j w k q 1 x A j Pl j w k x d lj w k x s A j m lj w k x
2 D lj w kq1 xy
.
Also 5 ji 5 2D Ž j. y AXj D Ž j . ji
2 A lj w k xy1
s 5 ji 5 2D lj w kq1xy .
MODAL ESTIMATION IN HYBRID SYSTEMS
247
Finally, y2 d lj w k x ⬘ A lj w k x
y1
AXj D Ž j . ji s 2 Ž d lj w k q 1 x
y
. ⬘ ji y y s y2 Ž m lj w k q 1 x . ⬘D lj w k q 1 x ji .
So, with some cancellation of terms, the result follows.
REFERENCES 1. B. D. O. Anderson and J. B. Moore, ‘‘Optimal Filtering,’’ Prentice Hall, New York, 1979. 2. H. A. P. Blom and Y. Bar-Shalom, The interacting multiple model algorithm for systems with Markov switching coefficients, IEEE Trans. Automat. Control 3, No. 8 Ž1988., 780᎐783. 3. R. J. Elliott, L. Aggoun, and J. B. Moore, ‘‘Hidden Markov Models: Estimation and Control,’’ Springer-Verlag, BerlinrNew York, 1995. 4. R. J. Elliott, F. Dufour, and D. D. Sworder, Exact hybrid filters in discrete time, IEEE Trans. Automat. Control. 41, No. 12 Ž1996., 1807᎐1810. 5. R. J. Elliott and J. van der Hoek, A finite dimensional filter for hybrid observations, IEEE Trans. Automat. Control AC-43 Ž1998., 736᎐739. 6. J. S. Evans and R. J. Evans, Image-enhanced multiple model tracking, Automatica 35 Ž1999., 1769᎐1786. 7. X. R. Li and Y. Bar-Shalom, Design of an interacting multiple model algorithm for air traffic control tracking, IEEE Trans. Control Systems Technol. 1, No. 3 Ž1993., 186᎐194. 8. M. Mariton, ‘‘Jump Linear Systems in Automatic Control,’’ Dekker, New York, 1990. 9. E. Mazor, A. Averbuch, Y. Bar-Shalom, and J. Dayan, Interacting multiple model methods in target tracking: a survey, IEEE Trans. Aerospace Electron. Systems 34, No. 1 Ž1998., 103᎐123. 10. D. D. Sworder and J. E. Boyd, Jump-diffusion processes in trackingrrecognition, IEEE Trans. Signal Processing 46, No. 1 Ž1998., 235᎐239. 11. D. D. Sworder and J. E. Boyd, ‘‘Estimation Problems in Hybrid Systems,’’ Cambridge Univ. Press, Cambridge, UK, 1999. 12. D. D. Sworder and J. E. Boyd, State estimation with biased observations, IEEE Trans. Systems Man Cybernet. A 29, No. 6 Ž1999.. 13. D. D. Sworder and R. O. Rogers, An LQ-solution to a control problem associated with a solar thermal central receiver, IEEE Trans. Automat. Control 28, No. 10 Ž1983., 971᎐978. 14. D. D. Sworder and R. Vojak, Hybrid estimation algorithms I, J. Optim. Theory Appl. 81, No. 7 Ž1994., 143᎐167.