Abstraction of continuous system trajectories into timed automata

Threshold

To represent the continuous system dynamics in a discrete formalism, the first step is to discretise the state space. For this, each c;, 's domain is partitioned

each location a timed constraint cP' , named

according to thresholds into a finite number of intervals that can be considered as qualitative states. Thresholds are defined from expert knowledge; e.g., a stock level can be qualified as " insufficient", "Iow", "medium", "high" and "critical". L; denotes

invariant. The system can remain in the same location as long as the invariant is true. Continuous system

the ordered set of thresholds in the i'h dimension, that is, with iE {l, ... ,n} and Q ; EN:

We consider a non-linear piecewise continuous system described by an ordinary differential equation (ODE) system such as:

( L,<)={ lI. o" ",l.~

(I)

where space,

~

ES

~

The continuous state space S is di vided according to the Q , thresholds by the mapping :

n

!Rn , with San-dimensional state

S=g(p ,t),

(EZ~Rm ,

D r :S

-""7 {

I, ... , Q ) }x .. .x{ I, ... ,Q, }

(8)

with Z a m-

dimensional input space, and pEP a parameter. It is

with the assumption

assumed that we have only a partial knowledge of both the initial state ~ (t o) and the input variables

rectangular partition is thus obtained.

values

(7)

I. ~&,

I

Sthat can be estimated by: 'il'i,~;-(to )$~,(to )$~,-(to)

(2)

'il'j,'il't , s )-(t)$(; (t)$(;(t)

(3)

((;,C;;)EZ) (S , cS and Z j cZ).

«

C;, E [l;.o ,l;n.

and ( >

according to the influence of the inter,al bounds on the derivative 's values (taking k ={l, ... ,m} \ (j}):

uncertainty by propagating the interval

f(; ,r;,,:;;-) ~ /,(; ,:;) ~ f(;,t;"r;;- ) = (~=:; ; ,:;: = C

J, 'il'i.

A

3. DISCRETE APPROXIMATION The idea is (i) not to consider as discrete states the cells partitioning the state space themselves but, rather, the cells ' faces and (ii) determine the state transitions from the simulation output of the continuous system within each cell as was done by Kowalewski et af. , 1999 . These authors defined a grid on each cell 's face and simulated all the trajectories from or bounded to each grid-point. In comparison, our approach depends upon each considered initial state and focuses on handling

with iE{I , ... ,n}, jE{I , ... ,m}, (~,-,~,-)ES/ and

Uncertainty representation. We define

that

(4)

within the cells by simulating Eq . (5).

/, (!;'(",C)? f (; ,S) ? J, (; ,;j;; )=;7= :;; ,':= i;;-

310

Lc;- ,~- ]

Trajectories within a cell

obtain T - , the ordered set of time instants when

Thresholds crossing. Let T be the time-point when the system (I) reaches v, a threshold value defined in the /h dimension, i.e., ~, (r) = v, vE lR . Because ~; is estimated by the interval [(- , Si-

crosses the thresholds and L-, the ordered set of threshold crossings:

J, the time window

within which crossing the v threshold may occur can be described by r- and r' such that {: (r- ) =~.- (r-) = v. From Eq. (5), the following implications can then be set:

(~:(r-»O) I\ (~: (r'»O)=H-
(9)

(~,- (r-)r + >r-

(10)

~­

(13)

( L- , ~)={a;

la; = £-(r;)}

(14)

J=card(r)=card(L-)

and

j =I , ... ,J ,

with

cv E

(y-,~)={r;I~-(r; )=I," }

{

I, ... ,.Q }. In the following, the mapping that

associates

'mm

or

'mar

to each 'J~ according to the

sign of the derivative is denoted A r :

Timed automaton representation. Let rmm=min(r-,r-) , rmax=max(r-,r ' ) and F the mapping that associates a label to each threshold crossing:

(15)

A similar approach can be used on the upper bound to obtain T-, L - and associated labels by the

(12)

mapping A r .

with the symbols b. and 'V denoting threshold crossing with increasing and decreasing trends respectively.

Building the timed automaton associated to the system 's trajectory in th e i'" dimension

r min and r max being approximated on Q (cf the

From Eqs . (13-15) , the following sets can be defined: - B, the ordered set of time instants when occur threshold crossings by ~ - and ~ - with label r m in

clock constraint definition) , the system trajectory between two thresholds wand v is represented within the time window [r mln ,rmax J. This window can be

(k8 = I, ... , K 8 , K8 = card( B) = card ( D)

modelled in the timed automata formalism as follows (cf Fig. I) : - S I and S2 are two locations with associated labels "w~" and "v"" respecti ve ly, - x is a clock, - x ~'max is the invariant of SI ,

):

(B,~) = {b,. E Y- [Ar (b,) =T mm }

U{b,.

(16) E

r [AT- (b,. ) =rmm}

- C, the ordered set of time instants when occur

- (s ,, 0,x ~Tm," , 0,sJ is the edge.

threshold crossings by

£- and

~-

with label

' max

(kc =l ,.. ,Kc , K c =card(C )): -, ~,-

' 0 ·'w !::."

w .j.-,l-----''--'--+_-

r-

x:::;

{C k, Er [Ar(C,) = 'm", }

~>T s: mm

(17)

"v::'"

U{CKcEr [Ar(Ck) =r

!"ma..r

mar }

- D , the ordered set of threshold associated to the elements of B:

r

Fig. \. Timed automaton for describing a system ' s trajectory from thresholds w to v with an increasing trend.

1£- (bk)

(D,$) = {dk• E r

U{ d,. E r

Generalisation

- Lo'ii ,

For sake of clarity, the index i is omitted in notations for all the variables in the following sections. Applying the above described approach to the lower bound of Eq. (5) in the iln dimension allows one to

the

decreasing

set

at

= dk.}

I{- (b,. )

of threshold

trends

(18) = d k.

}

crossings

with

initial

time

the

(k,;" =I ,... ,K,;" ' K,;< =card ( Xn):

31 1

crossings

I o'1

={a~k)(a~k,. E L) A (~-
EL)

A (

~

-

=

a~k" )

A ((

::;

o)} (19)

U{a~k)(a~\, EL) A (~- =a~k" )A( ~ + ~o)}

- Lt

'::' j

- i

v+-----+-~~~--_

w~~r_~--------

10

the set of threshold crossings with an

increasing

trend

(k~, =1 ,... , K~"

at

the

initial

a)

time

b)

-,

- i

;:+

K~" =card (I t))

I: ={a~k) (a~k" E L)

A

(~- < a~k"

'm in

~ ,-

'> i

,

< ~-)}

c.-.,

v ~~----~~----

w

U{a~k.l( a~k,E L) A( C= a~kJ A (~-~ a)}

(20) 10

U{a~k~,1 ( a~k"' EL )A (~- = a~k~') A (~-::; o)}

'fm in

c.i)

c.ii)

Fig . 2. Different cases of threshold crossing with increasing trend.

Clock set. The clock set is a singleton since the automaton is obtained from only one simulation of the continuous model.

Based on thi s, a set of locations can be obtained which satisfy one of these conditions . The preceding location is the element to which is associated the

Location set. Let us denote : - S , a mapping that associates a location S to each

nearest

element of D , I o'1 and I :. From S, one gets

S={S k}k;) . "K U SO wi th

~,-

V

T mm

of the

Tm m

of

sn'

Fonnally, al. the

forward edge of sn is defined by the 5-tuple

K=Ks+ K~,+K~ .

(sm ,0, rp ,0, sn) wi th rp:=x~bn

and So, the location corresponding to the initial state of the continuous system . - P , a mapping that associates a proposition to each element of S denoting which threshold is crossed and the trend (with WE Q and with P (so) ~ 10 ):

1

Let the ordered

set S' cS , the union of the different cases described by the Fig. 2:

_ j(P (sn) =I~) = {~'I(P (S') =I~_' )}

'<

(S;> _) -

(23)

(P(Sn) = I: ) = {s,l( p(s,) = I,: .))}

(S( d,,) ~ s, ) (d" =1",) /\ (d" =r ) :=) I~ (s( d,,)~ s, ) /\ (d" =1,,) /\ (d" =;[ -) :=) I~ A

P(s

,

):~

(s( a;,", )~ S,) /\ (a;,", =/,") :=) I~ (s( a~,", )~ s, ) /\ (a~,", = Iw ):=) I~

'<

_ j( P ( sn) =I~) :=) {s'I( P (S') =I:)}

(S,._ ) -

.

(21)

(24)

( P ( sn ) = () :=){s, I(P(S,) =I~ )}

S; = (sn= S (a;,))) v (sn=S (ag,. )) = {so}

(25)

(sn=S (d" ))/\( P (S")=I~)/\

Forward edges. Edges that can reach a new location at the earliest time point T rn in are named F on"ard

S:=

edges. Let Sn.I1E Ks , a location corresponding to a threshold crossing w ith increasing trend, and its associated proposition I~ . The preceding location can

(1,,-, <~-(lo )::;,;'(lo) ::; lw) = {so}

(26)

(s"=S (d" ))/\( P(sn)=C) /\

(1," <~ - (lo )::;~'(lo )::; lw _ ' ) = {so} (S',::;) = S; u S; u S; u S;

reflect the fo llowing cases (cf. Fig . 2): a) crossing a lower threshold with increasing trend, b) crossing the same threshold with decreasing trend, c) the initial state whenever no threshold has already been crossed and (i) the threshold is within the initial state interval , or (ii) the initial state interval lies below the thresho ld. The same reasoning can be made the other way around about threshold crossing wi th a decreasing trend (i.e . about a location \\iith a proposition I: ).

(27)

with k = 1, .. ., n -I (i.e., only reachable locati ons before Sn are taken). The preceding location s rn <;;;; S ' is the latest that could be reached. The set of locations of S forward edges is named A). A similar approach can be used about the latest time points (i. e. , thresh old crossing at T max ) to determine I If the location is an ini tial state threshold crossing. th en the guard is of the form x = I, .

312

the location's invariants. Moreover, when a location has no next location, a Backward edge is defined. Figure 3 shows this case: the system crosses a threshold successively twice, and the edge with guard x= is a backward edge. Doing this, the Inv

r:,

mapping, that associates the invariants to the locations, and A z , the set of backward edges , are defined. A complete description of this approach can be found in Helias, (2003). -=i

I

growth rates of the bacterial biomass X l according to 51 (in the acidogenic phase) and bacterial biomass X 2 according to 52 (in the methanogenic phase) respectively. These growth rates are assumed to

5 follow a Monod-type kinetics, i.e, 11 = I1max 5 ... K s

v+--ri-,_-_~l--~'-,_----,-'~'i~ " 1----~-

/"1

P

co, . 51> 52, Z, and C are the PT -Pea, organic matter concentration, the volatile fatty acids concentration, the total alkalinity and the total inorganic carbon, respectively. k1 to k6 are the yield coefficients for the reactions and 111and 112 are the with kp = ks -k6

with I1max the maximum bacterial growth rate and Ks the half-saturation constant. Finally, Pr is the

r

total pressure in the reactor,

T

Pea, is the partial

\ '!

Fig. 3. Backward edge example. From (16-27) , the continuous system dynamics in the i'h dimension IS approximated by the timed automaton G i =(S" Pi, X" E i,A"Inv i ) where: - S, =

{Sk,} _ Uso., k, - I .... . K,

is a set of locations,

- Pi is a mapping that determines, for each location , the threshold crossing and its trend, - Xi is a singleton composed of the clock Xi, - E i, the set of labels, is an empty set, - A, = AI., U A 2 ., is the set of edges, -lnvi is a mapping that associates the invariants to the locations. 4. EXAMPLE: APPROXIMATION OF A WASTEW A TER ANAEROBIC DIGESTION MODEL To briefly illustrate this approach , the procedure previously described is used on an anaerobic digestion process. A more complete description of this example can be found in Helias et al., (2004). The aim is to predict possible dysfunctions of a wastewater treatment plant by a discrete approximation of a biological reaction model in front of expert partial knowledge of the system ' s dYnamics, Anaerobic digestion is a biological process u~ed for carbon removal from wastewater. The principle is to transform organic matter into biogas (i.e., that is a methane (CH 4 ) and carbon dioxide (C0 2) mixture) in absence of oxygen. This complex process can be modelled in two steps: the acidogenic and methanogenic phases . A simplified version of the mass balance analytical model of our process (Bemard et af. , 200 I) is used:

=D (5 ],n -

51) - kll1l X I

= D( 52,n - 52)';" k2111X I - kJI12X2 =

D(Z," -Z)

=

D(C,n-C),. k 4111XI + kpl12X2

pressure of CO 2 , D is the dilution rate according to input flow and 51in , 52," , Z in and C in the different influent concentrations . Because of a relation holding between Pr and Peo" the cooperativity property holds in this model. Thus, intervals bounding influent and biomass concentrations can be used to handle imprecision, following the approach described in Section 2.. This model is run with parameters measured on an anaerobic digestion pilot plant (Steyer et af. , 2002). The inputs of the process (i.e. , influent concentrations) and the initial state are thus represented by intervals. Qualitative states resulting from expert knowledge (e.g., 52 can be assigned labels "normal", "high" or "critical") also available are used in this approach. The plant operation is apprehended on a short time period (a few days) by combining information on different variables. For example, an important organic overload (i.e. dysfunction due to an inhibition of the reaction by 52) can be associated to a " low" Z value and a "high" 52 value according to thresholds specifically set in these variables' domains. To check the possible occurrence of this situation, we use the model-checking tool on the obtained timed automaton . This approach has been used to check about fifteen operating states during a four days period (Helias et aI., 2004) . The complete abstraction procedure (system simulation, threshold crossing detection, timed automata creation), with 11 thresholds defined on 4 variable domains, takes less than ten seconds on a personal computer with a CPU running at 1.7 GHz. This procedure has been implemented by combining the Matlab ': environment and the Kronos model-checker (Yovine , 1997). The resulting automaton IS composed of 24 locations and 66 edges. 5. DISCUSSION

(28)

The final size of the automaton is set according to the svstem's dimension, the number of thresholds and the crossing number. Consequently, this approach is

3 13

inappropriate when: - the system oscillates around a threshold value, - a large number of thresholds are defined. The idea is not to create an automaton to check all the possible properties of the continuous system but, rather, to create the automaton specifically suited to check a given property according to purposely defined thresholds.

REFERENCES R. Alur and Dill D (1994). The theory of timed automata. Theor. Comp . Sc .. 126: pp. 183-235. O. Bemard, Hadj-Sadok Z., Dochain D., Genovesi A. and Steyer J.-P (2001). Dynamical model development and parameter identification for anaerobic wastewater treatment process. Biotech. and Bioengineering. 75(4): pp. 424-438. J.L. Gouze, Rapaport A. and Hadj-Sadok M.Z (2000). Interval observers for uncertain biological systems. Ecological Modelling. 133(1-2): pp. 45-56. A. Helias (2003). Agregation/abstraction de modetes pour l'analyse et {'organisation de reseaux de flux application a la gestion des effluents d'etevage a la Reunion. PhD thesis. ENSAM. Montpellier (France). 224 p. A. Helias, Guerrin F. and Steyer J.-P (2004). Abstracting continuous system behaviours into timed automata: application to diagnosis of an anaerobic digestion process. In Proc. DX-2004 15th International Workshop on Principles of Diagnosis. Carcassonne (France). TA Henzinger, Ho P.-H. and Wong-Toi H (1998). Algorithmic Analysis of Nonlinear Hybrid Systems. IEEE Trans. on Automatic Control. 43(4): pp. 540-554. S. Kowalewski, Engell S., Prui3ig J. and Stursberg 0 (1999). Verification of logic controller for continuous plants using timed condition/eventsystem models. Automatica. 35(3): pp. 505-518. T. , Moor, and Raisch, J (2002) . Abstraction based supervisory controller synthesis for high order monotone continuous systems. LNCIS 279: pp. 247-265. J. Lunze (1999). A timed discrete-event abstraction of the continuous variable systems. International Journal of Control, 72(13), pp . 1147-1164. H.L. Smith (1995). Monotone dynamical systems: An introduction to the theory of competitive and cooperative systems. Mathematical Surveys and Monographs. 41. 174 p. l-P . Steyer, Bouvier J.-c. , Conte T. , Gras P. and Sousbie P (2002). Evaluation of a four year experience with a fully instrumented anaerobic Water Science and digestion process. Technology. 45(4-5): pp. 495-502 . P. Struss (2002). Automated Abstraction of Numerical Simulation Models - Theory and Practical Experience. in Model Based Systems and Qualitative Reasoning for Intelligent Tutoring Systems. San Sebastian (Spain). pp. 161-168. P. Supavatanakul, Falkenberg C. and Lunze J (2003) . Identification of timed discrete event models for diagnosis. in DX-2003 14th International Workshop on Principles of DiagnosiS. Washington DC (USA). S. Yovine (1997). Kronos: A verification tool for real-Time Systems. Journal of SofMare Tools for Technology Transfer . 1(1-2): pp. 123-133

Whereas in the classical definition of timed automata the clock constraints compare the clock values using rational numbers, integers must be used in the Kronos model-checker. Consequently, m our approach, the threshold crossing time values must be rounded. Thus, if the time unit used in the model is not carefully chosen according to the system dynamics, some additional imprecision may be introduced.

6. CONCLUSION This paper has described a procedure aimed at representing a continuous system ' s dynamics as a timed automaton, where: - some elements (inputs, initial state) are estimated by intervals, which comes down to reformulate the system as extremal ODE systems (upper and lower bounds ~- and ~- approximating the real dynamics), - dynamical envelopes stemming from the simulation of this extremal ODE systems are translated into the timed automata formalism (discrete model with continuous time) according to expert defined thresholds on the state variables ' domain. The principal characteristics of this approach are: - taking threshold crossings themselves as qualitative discrete states (disregarding the interval between two consecutive thresholds); - defining the discrete state transitions from the piecewise simulation of the extremal ODE systems within consecutive thresholds ; - starting the procedure from (i) an initial state, (ii) imprecise values of inputs and (iii) threshold definitions, all defined in accordance to the property to be checked in order to reduce the state explosion. This approach was briefly illustrated on an anaerobic digestion model. The timed discrete-event representation makes model-checking possible to check properties of the continuous system ' s dynamics. An interesting extension would be to couple this representation to the discrete part of a wastewater treatment process in order to check properties on mixed continuous and discrete elements. ACKNOWLEDGEMENTS This work was partially supported by the Region Reunion and the European Social Fund.

314