CC BY
12
A MONTE CARLO-BASED TECHNIQUE FOR ESTIMATING THE OPERATION MODES OF HYBRID DYNAMIC SYSTEMS
F. Cadini, D. Avram, E. Zio. •
Politecnico di Milano - Dipartimento di Energia Milan, Italy
e-mail: francesco.cadini@polimi. it
ABSTRACT
Many real systems are characterized by a hybrid dynamics of transitions among discrete modes of operation, each one giving rise to a specific continuous dynamics of evolution. The estimation of the state of these hybrid dynamic systems is difficult because it requires keeping track of the transitions among the multiple modes of system dynamics corresponding to the different modes of operation. A Monte Carlo-based estimation method is here illustrated through an application to a case study of literature.
1 INTRODUCTION
Diagnosis and prognosis of system faults rely on the knowledge or anticipation of the system state to provide advanced warning and lead time for preparing the necessary corrective actions to maintain the system in safe operation.
The related state estimation task becomes quite challenging for systems with a hybrid dynamic behavior characterized by continuous states and discrete modes. Sudden transitions of the discrete modes, often autonomously triggered by the continuous dynamics, affect the system evolution and a large computational effort is required to keep track of the multiple models of the discrete system modes and the autonomous transitions between them (Koutsoukos et al. 2002). Since the dynamic states cannot be directly observed, the problem becomes that of inferring the system state from related measured parameters.
The soundest model-based approaches to the estimation of the state of a dynamic system or component build a posterior probability distribution of the unknown states by combining the probability distribution assigned a priori to the possible states with the likelihood of the observations of the measurements actually collected (Doucet 1998, Doucet et al. 2001). In this Bayesian setting, the estimation method most frequently used in practice is the Kalman filter, which is optimal for linear state space models and independent, additive Gaussian noises. In practice, however, the dynamic evolution of many systems and components is non-linear and the associated noises are non-Gaussian (Kitagawa 1987). In these cases, one may resort to Monte Carlo sampling methods also known as particle filtering methods, which are capable of approximating the continuous distributions of interest by a discrete set of weighed 'particles' representing random trajectories of system evolution in the state space and whose weights are estimates of the probabilities of the trajectories (Doucet et al. 2000, Djuric et al. 2003, Cadini et al. 2009a, b).
In this paper, particle filtering is applied for the estimation of the state of a hybrid system of literature often taken as a benchmark for dynamic reliability estimation and fault diagnosis/prognosis methods (Aldemir et al. 1994, Marseguerra et al. 1996, Labeau et al. 1998,
Wang et al. 2002). The system consists of a tank filled with a liquid whose level is autonomously maintained between two thresholds by actuators driving three filling and emptying flows triggered by the actual liquid level. The actuators discrete mode is estimated by the particle filter on the basis of noisy level and temperature measurements.
2 PARTICLE FILTERING FOR OPERATION MODE ESTIMATION
2.1 General model-based framework for state estimation
Let us consider a continuous system whose evolution can be described by:
df> o
where x is the system state vector, f: R"x x R"" ^ R"x is possibly non-linear and « is an independent identically distributed (i.i.d.) state noise vector of known distribution.
The state x cannot in general be directly observed; rather, information about x can be inferred from the observation of a related variable z whose relation to the state x is described in general terms by the equation:
z = h(x, u) (2)
where h: R"x x R"" ^ R"x is possibly non-linear and u is an i.i.d. measurement noise vector sequence of known distribution. The measurements z are, thus, assumed to be conditionally independent given the state process x.
The practical implementation of computational tools for state estimation requires that the continuous system dynamics be discretized appropriately. Regardless of the discretisation method adopted, the system state dynamics can be represented by an unobserved (hidden) Markov process of order one:
x k = f (x t-i.ra k-i) (3)
where fk : R"x x R"" ^ R"x is possibly non-linear and {wk, ke N} is an independent identically distributed (i.i.d.) state noise vector sequence of known distribution.
The transition probability distribution p(xk\xk-i) is defined by the system Equation (3) and the known distribution of the noise vector wk. The initial distribution of the system state p(x0) is assumed known.
A sequence of measurements {zk, ke N} is assumed to be collected at the successive time steps tk. The sequence of measurement values is described by the measurement (observation) equation:
z k = h k (x k, U ) (4)
where hk : R"x x R"" ^ R"x is possibly non-linear and { uk, ke N} is an i.i.d. measurement noise
vector sequence of known distribution. The measurements {zk, ke N} are, thus, assumed to be conditionally independent given the state process {xk, ke N}.
Within a Bayesian framework, the filtered posterior distribution p(xk|z0:k) can be recursively computed in two stages: prediction and update (Doucet 1998, Arulampalam et al. 2002). Given the probability distribution p(xk-1|z0:k-1) at time k-1, the prediction stage involves using the system model (3) to obtain the prior probability distribution of the system state xk at time k via the Chapman-Kolmogorov equation:
p(xk Iz 0k-1 )= j p(xk k-1 Iz 0*-1 )p(xk-1 Iz 0k-1 )^xk-1 = J p(xk |xk-1 )p(xk-1 Iz 0:k-1 )^xk-1 (5)
where the Markovian assumption underpinning the system model (3) has been used.
At time k, a new measurement zk is collected and used to update the prior distribution via Bayes rule, so as to obtain the required posterior distribution of the current state xk (Arulampalam et al. 2002):
i I \ P(xk N 0:k-1 )P(zk Ixk) P(x k z 0:k )= -T"!-r--(6)
P(z k z 0:k-1 )
where the normalizing constant is
p(z
k z 0:k-1 )=i P(x k z 0:k-1 )p(z k\xk Vxk (7)
The recurrence relations (5) and (6) form the basis for the exact Bayesian solution. Unfortunately, except for a few cases, including linear Gaussian state space models (Kalman filter) and hidden finite-state space Markov chains (Wohnam filter), it is not possible to evaluate analytically these distributions, since they require the evaluation of complex high-dimensional integrals.
One way to overcome this problem is to resort to Monte Carlo sampling or PF methods (Pulkkinen 1991, Doucet et al. 2000, Doucet et al. 2001, Seong et al. 2002). Assuming that a set of random samples (particles) x0:k, i = 1, 2,..., Ns, of the system state at the time k-1 is available as a realization of the posterior probability p(xk-1|z0:k-1), the predicting step at time k is accomplished by sampling from the probability distribution of the system noise wk-1 and simulating the system dynamics (3) to generate a new set of samples xik which are realizations of the predicted probability distribution p(xk|z0:k-1).
In the update step, based on the likelihoods of the observations zk collected at time k, each sampled particle x k-1 is assigned a weight:
P№) (8)
Zp(z k |x k)
j=i
wk =
An approximation of the posterior distribution p(xk|z0:k) can be obtained in terms of the weighted samples (xk, wik), i = 1, 2,., Ns (Doucet et al. 2001).
One difficulty that arises in the implementation of PF is the degeneracy problem: as the algorithm evolves in time, the weight variance increases (Doucet 1998) and the importance weight distribution becomes progressively skewed, until (after a few iterations) all but one particle have negligible weights (Arulampalam et al. 2000, Doucet et al. 2000, Andrieu et al. 2001). As a result, the approximation of the target distribution p(xk|z0:k) becomes very poor and significant computational resources are spent trying to update particles with minimum relevance. A possible
solution to this problem is to proceed to a resampling of a new swarm of realizations xk from the approximate posterior distribution, constructed on the weighted samples previously drawn; all particles thereby generated are assigned equal weights, w'k = 1/Ns (Doucet et al. 2001).
As the final step, one has to resample from the posterior distribution a new swarm of points xk. The prediction, update and resample steps form a single iteration, recursively applied at each time k.
2.2 Hybrid system model
Let us consider a hybrid system whose dynamic evolution can be described by:
íx^ = fp (x^-i> «^-1 ) [Pk = g k (fik-i,x k-i)
k 'PkVk-i,^ k-i, (9)
pk = {1,2,3,...M} is the discrete state which indicates the mode in which the system is evolving at time k, fgk is the non-linear function describing the (discretized) continuous evolution of system state x when the system is in mode at time k, gk is the discrete mode transition function. In what follows, we shall consider only autonomous transitions between the system modes, i.e. those triggered by the control of the continuous state x which demands transitions among the system modes when reaching specified thresholds.
Let s'k = (P'k,x'k) indicate the 1th sample of the extended hybrid system state, where x'k is the random sample drawn from the importance functionp(xk\xk-i) and f?k is the corresponding discrete mode of system behavior. Then, the posterior probability density of the continuous and discrete states can be represented by the random measure {s1k,wk,i=i...Ns}, where wk is the particle weight of the 1th sample of the hybrid state at time k after resampling.
The estimation of system mode of operation as the most likely one is given by:
Pk = arg max Z wk (i0)
ieGj
where Gj = {pi = j\. Whereas, the posterior estimate mean of the continuous state xk and its variance ó\ are given by:
= Izó Wk x k Tizó Wk
Ns ■ Í ■ \
I wk (xk - xk )
_ ie<3
' k
(ii)
at = ^—— (i2)
Izó wk
where only the particles belonging to the most likely mode ¡3lk are considered Gj = { p'k = ¡3lk j.
3 APPLICATION TO A TANK CONTROL SYSTEM
The particle filter estimation algorithm is applied to a hybrid system of literature (Aldemir et al. 1994, Marseguerra et al. 1996, Labeau et al. 1998, Wang et al. 2002). The system consists of a tank containing a fluid whose level is controlled by three control units which open or close
depending on the fluid level crossing of predefined thresholds (HLVand HLP) (Figure 1). The fluid in the tank is uniformly heated, under adiabatic conditions, by a thermal power source W.
Figure 1. Tank control system (Aldemir et al. 1994, Marseguerra et al. 1996, Labeau et al. 1998, Wang et al.
2002).
The control aims at maintaining the fluid level X1 in the range (x1,min — HLV, X1,max — HLP), while also monitoring the fluid temperature x2 which may become relevant from a safety point of view.
The operational states of the control units at time k are described by the Boolean indicator a/;k, l = 1,2,3, where a/k assumes the value 1 or 0 according to whether the unit is on (a/,k = 1) or off (a/,k = 0). The autonomous control actions modify the states a,k of the units according to the following rules:
1 if x < HLV 0 if x > HLP
0 or 1 dependingon previousswitching
a2,k = '
[1 if x > HLV 0 otherwise
(13)
«1,k = <
1 if x < HLV 0 if x > HLP
0 or 1 dependingon previousswitching
Thus, the following four modes of system dynamic evolution may be identified:
Pk =
1 if xu < HLV
2 if HLV < xu < HLP and = axk_x = 1
3 if HLV < xxk < HLP and ^ ^ = a3k_x = 0
4 if xlk > HLP
(14)
«1,k = <
<
With the additional simplifying physical assumptions that the fluid input in the tank by units 1 and 3 mixes instantaneously and the flow rate through the outlet unit 2 is independent of the fluid level, and the discretisation of the system dynamics, the time evolution of the state x1kk and x2,k can be described by two first-order, decoupled, non-linear difference equations determined by the mass and energy conservation laws (Wang et al. 2002).
The aim of the analysis is that of estimating the discrete mode of the system, i.e. the operational states of the three control units on the basis of Ns trajectories drawn from the system model and a sequence of noisy measurements of the level x1,k and the temperature x2,k:
Z1,k = X1,k + V\,k
Z2,k X2,k + V2,k
(15)
where v1,k and u2,k are the measurement noises. Knowledge of the system mode of operation allows the proper control and maintenance of its components.
Let us suppose that the control system starts from x1;o =6m and x2,0 =10m. The time horizon considered for the evolution of the system dynamics is Nt = 40h, with level and temperature observations at discrete time steps of At = 30min (Nk = 80). As in the application of reference (Wang et al. 2002), the inlet fluid temperature is &m = 15°C, the level thresholds are set at HLV = 4m and HLP = 10m and the fluid flow rates are Q1 = lm/h, Q2 = 4m/h and Q3 = 4.5m/h. A zero -
.... 9
mean Gaussian noise with variance o q = 0.0025 is added to the flow rates, for closer adherence to
reality. The process and the measurement noises are assumed Gaussian with zero mean and 2 2 variances o m = [0.02 0.01] and o v = [0.16 0.05] respectively.
Assuming independence of the level and temperature measurements, the observation
likelihood in (8) can be written as:
\2 f \2 11 z1,k~Pl,k I 1| z2,k~M2,k
p(z k| xk )=n pzul xk )=^ - ' (16)
h=1,2 y]2ft<Jlt yl2K<J,
First, a crude, measurement-based, empirical algorithm is proposed for the estimation of the mode at time k:
Pk =
1 if Z1,k < HLV
2 if HLV < zlk < HLP and = a3k_} = 1
3 if HLV < zlk < HLP and axk_x = aik_x = 0
4 if Z1,k > HLP
(17)
<
where z\kk is the level measurement at time k.
Figure 2 shows the estimated mode ¡3 (dot-dashed line) and the model simulated one fi (solid line). The performance is not satisfactory because the noise v\ generates spurious oscillations in the level measurement z\ with respect to the model-simulated x\ actually driving the mode transitions.
Figure 2. Measurement-based estimated system modes (dotted line) and model-based simulated system modes
(solid line).
To overcome this problem, the particle filter method is implemented with a number of particles Ns — 1000 (Figure 3). Figure 4 shows the particle filter-estimated mode ¡3 (dot-dashed line) and the model simulated one ( (solid line). The agreement is satisfactory, with the only exception in correspondence of the first time when the system enters mode ( — 4, i.e. the fluid level is higher than HLP. This is due to the fact that the first few observations of the fluid level higher than HLP do not provide the filter with enough information for properly performing the mode estimation. This is confirmed in Figure 5, where the estimated level x (dotted) is affected by a larger uncertainty ¿rx k when approaching the threshold HLP for the first time.
'Jl
i H
pf Wil
by4'* j>p
f
* Simulate Measure d level x^ merits z.
--ko band I i
'0 5 ID 15 20 25 30 35 40
Time [h]
Figure 3. Fluid level measurements (dotted line), with measurement noise uncertainty ± 1ct bands (solid line);
model-simulated fluid level (dots).
Figure 4. Particle filter-estimated (dotted line) and model-simulated (solid line) modes.
'0 5 ID 15 20 25 30 35 40 Time [h]
Figure 5. Particle filter-estimated mean fluid level (dotted line), with ± uncertainty bands (solid line) and
model-simulated fluid level (dots).
0 5 10 15 20 25 30 35 40
Time [h]
Figure 6. Particle filter-estimated mean of the fluid temperature (dotted line), with ± uncertainty bands (solid line) and simulated fluid temperature (dots).
Similar satisfactory results (not reported for brevity's sake) have been obtained in the estimation of the temperature state variable, x2.
4 CONCLUSIONS
In this paper, a Monte Carlo-based filter has been devised for estimating both the continuous states and the discrete modes of a controlled system, whose transitions between the discrete modes are autonomously triggered by the continuous states. Comparison with a crude algorithm which bases its estimates directly on the observed measurements, shows the higher performance of the particle filter on a wider range of measurements noises, thus counterbalancing the larger computational effort required.
5 PREFERENCES
Aldemir, T., Siu, N., Mosleh, A., Cacciabue, P.C. & Goktepe, B.G. (1994). Eds.: Reliability and Safety Assessment of
Dynamic Process Systems. NATO-ASI Series F. 120, Berlin: Springer-Verlag. Andrieu, C., Doucet, A., & Punskaya, E. (2001) .Sequential Monte Carlo Methods for Optimal Filtering. Sequential
Monte Carlo Methods in Practice, Doucet, A., de Freitas, N., and Gordon, N., Eds. New York: Springer-Verlag. Arulampalam, M.S., Maskell, S., Gordon, N. & Clapp, T. (2002). A Tutorial on Particle Filters for Online
Nonlinear/Non-Gaussian Bayesian Tracking. IEEE Trans. On Signal Processing. 50, 2, 174-188. Cadini, F., Zio, E. & Avram, D. (2009a). Model-based Monte Carlo state estimation for condition-based component
replacement. Reliab. Eng. and Sys. Saf. 94, 752-758. Cadini, F., Zio, E. & Avram, D. (2009b). Monte Carlo-based filtering for fatigue crack growth estimation. Probabilistic
Engineering Mechanics. 24, 367-373. Djuric, P.M., Kotecha, J.H., Zhang, J., Huang, Y., Ghirmai, T., Bugallo, M. F. & Miguez, J. (2003). Particle Filtering,
IEEE Signal Processing Magazine. 19-37. Doucet, A. (1998). On Sequential Simulation-Based Methods for Bayesian Filtering. Technical Report, University of
Cambridge, Dept. of Engineering, CUED-F-ENG-TR310. Doucet, A., Godsill, S. & Andreu, C. (2000). On Sequential Monte Carlo Sampling Methods for Bayesian Filtering.
Statistics and Computing. 10, 197-208. Doucet, A., de Freitas, J.F.G. & Gordon, N.J. (2001). An Introduction to Sequential Monte Carlo Methods, in Sequential Monte Carlo in Practice. A. Doucet, J.F.G. de Freitas and N.J. Gordon, Eds., New York: SpringerVerlag,.
Kitagawa, G. (1987). Non-Gaussian State-Space Modeling of Nonstationary Time Series. Journal of the American
Statistical Association. 82, 1032-1063. Koutsoukos, X., Kurien, J. & Zhao, F. (2002). Monitoring and diagnosis of hybrid systems using particle filtering methods. Proceedings of the 15th International Symposium on the Mathematical Theory of Networks and Systems (MTNS).
Labeau, P.E. & Zio, E. (1998). The Cell-to-Boundary Method in the Frame of Memorization-Based Monte Carlo Algorithms. A New Computational Improvement in Dynamic Reliability. Mathematics and Computers in Simulation. 47, 2-5, 329-347.
Marseguerra, M. & Zio, E. (1996). Monte Carlo approach to PSA for dynamic process systems. Reliab. Eng. and Sys. Saf. 52, 227-241.
Pulkkinen, U. (1991). A Stochastic Model for Wear Prediction through Condition Monitoring. K. Holmberg & A.
Folkeson Eds.. Operational Reliability and Systematic Maintenance. London/New York: Elsevier. 223-243. Seong, S.-H., Park, H.-Y., Kim, D.-H., Suh, Y.-S., Hur, S., Koo, I.-S., Lee, U.-C.,. Jang, J.-W & Shin, Y.-C. (2002). Development of Fast Running Simulation Methodology Using Neural Networks for Load Follow Operation. Nuclear Science and Engineering. 141, 66-77. Wang, P., Chen X.M. and Aldemir T. (2002). DSD: a Generic Software Package for Model-Based Fault Diagnosis in Dynamic Systems. Reliab. Eng. and Sys. Saf. 75, 31-39.
