Zsófia Lendek
email: zsofia.lendek@aut.utcluj.ro

Paula Raica
email: paula.raica@aut.utcluj.ro

Jimmy Lauber
email: jlauber@univ-valenciennes.fr

Thierry Marie Guerra
email: guerra@univ-valenciennes.fr

Observer design for switching nonlinear systems

Keywords: switching systems, Takagi-Sugeno fuzzy models, observer design, periodic systems, nonquadratic Lyapunov function

In this paper we consider observer design for discrete-time switching systems represented by Takagi-Sugeno fuzzy models. For the design we use a switching Lyapunov function defined on the switches. In order to develop stability conditions for the estimation error dynamics, we consider the variation of this function along possible switches. The conservativeness of the approach is reduced by considering the αsample variation of the Lyapunov function. This approach can bring solutions to observer design for some switching systems with unobservable and unstable local models. The developed conditions are illustrated on a numerical example.

I. INTRODUCTION

Switching systems are usually described by continuous and discrete dynamics, and their interactions. For such systems, one can analyze the stability [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF] of the whole system assuming that the switching sequence is known or not known and design switching observers and controllers. In this paper we consider observer design for switching nonlinear systems represented by Takagi-Sugeno (TS) fuzzy models.

Takagi-Sugeno (TS) fuzzy systems [START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF] are convex combinations of local linear models. For the analysis and synthesis of TS systems the direct Lyapunov approach has been used. The design conditions are derived in the form of linear matrix inequalities (LMIs), which can be solved using available algorithms. To derive the conditions, next to classical quadratic Lyapunov functions [START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF]- [START_REF] Sala | Perspectives of fuzzy systems and control[END_REF], recently nonquadratic Lyapunov functions [START_REF] Guerra | LMI-based relaxed nonquadratic stabilization conditions for nonlinear systems in the Takagi-Sugeno's form[END_REF]- [START_REF] Mozelli | Reducing conservativeness in recent stability conditions of TS fuzzy systems[END_REF] have been used, which have shown a real improvement in the design conditions [START_REF] Guerra | LMI-based relaxed nonquadratic stabilization conditions for nonlinear systems in the Takagi-Sugeno's form[END_REF], [START_REF] Ding | Further studies on LMI-based relaxed stabilization conditions for nonlinear systems in Takagi-Sugeno's form[END_REF]- [START_REF] Lee | Approaches to extended nonquadratic stability and stabilization conditions for discrete-time Takagi-Sugeno fuzzy systems[END_REF].

Switching TS systems have been investigated in the last decades mainly in the continuous case where the stability is based on the use of a quadratic Lyapunov function [START_REF] Tanaka | Switching control of an R/C hovercraft: stabilization and smooth switching[END_REF]- [START_REF] Ohtake | Switching fuzzy controller design based on switching Lyapunov function for a class of nonlinear systems[END_REF] or a piecewise one [START_REF]Controller synthesis of fuzzy dynamic systems based on piecewise Lyapunov functions[END_REF], [START_REF]H∞ controller design of fuzzy dynamic systems based on piecewise Lyapunov functions[END_REF]. For discrete-time TS models, the number of results is limited [START_REF] Doo | Output-feedback control of discrete-time switching fuzzy system[END_REF], [START_REF] Dong | H∞ controller synthesis via switched PDC scheme for discrete-time T-S fuzzy systems[END_REF].

Most of the existing results on switching systems concern linear subsystems, such as the one in [START_REF] Jungers | Finite L 2 -induced gain and λ-contractivity of discrete-time switching systems including modal nonlinearities and actuator saturations[END_REF], where stabilization in the presence of input saturation and uncertainties is considered, or [START_REF] Dehghan | Computations of mode-dependent dwell times for discrete-time switching system[END_REF] which considers the computation of the modedependent dwell-time. Observability for switching discretetime linear systems has been investigated in [START_REF] Santis | On location observability notions for switching systems[END_REF].

In this paper, we derive relaxed conditions for the observer design of switching TS systems. For this, we use a switching nonquadratic Lyapunov function and make use of its variation over possible switches. We assume that although the exact switching sequence is not known, the set of all the admissible switches is known. Furthermore, once a subsystem is activated, it will remain active for a number of samples, for which minimum and maximum bounds are known.

The structure of the paper is as follows. Section II presents the notations used in this paper. The development of observer design conditions is presented in Section III. Section IV discusses the developed conditions and illustrates their use on a numerical example. Section V concludes the paper.

II. PRELIMINARIES

In this paper we consider observer design for discrete-time switching TS systems. We consider subsystems of the form

x k+1 = r i=1 h i (z k )(A j,i x k + B j,i u k ) y k = r i=1 h i (z k )C j,i x k
or, in a shorthand notation

x k+1 = A j,z x k + B j,z u k y k = C j,z x k ( 1
)
where j is the number of the subsystem, j = 1, 2, . . . , n s , n s being the number of the subsystems, x denotes the state vector, r is the number of rules, z denotes the scheduling vector, h i , i = 1, 2, . . . , r are normalized membership functions, and A j,i , B j,i , C j,i , i = 1, 2, . . . , r, j = 1, 2, . . . , n s , are the local models. The scheduling vector is assumed to be available online, but not in advance.

For the easier notation, we use a directed graph representation of the switching system [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF]. The graph associated to (1) is G = {V, E}, where V denotes the set of vertices or subsystems and E denotes the set of admissible switches. As such, (v i , v j ) ∈ E if a switch from subsystem i to subsystem j is possible. Note that we assume that self-transitions are also possible: these correspond to the subsystem being active for more than one sample. A path P(v i , v j ) between two vertices v i and v j in the graph G is a sequence of vertices P(v i , v j ) = [v p1 , v p2 , . . . , v pn p ] so that v i = v p1 , v j = v pn p , and (v p k , v p k+1 ) ∈ E, p k = 1, 2, . . . , n p -1. A path in a graph associated to a switching system corresponds to a switching law. The length of a path is given by the number of edges it contains.

The notations above are illustrated on the following example.

Example 1. Consider a switching system composed of 4 subsystem:

x k+1 = A i,z x k for i = 1, 2, 3, 4
, and with admissible switches (2, 1), (2, 3), [START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF][START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF], [START_REF] Tanaka | Fuzzy Control System Design and Analysis: A Linear Matrix Inequality Approach[END_REF][START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF], [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF][START_REF] Tanaka | Fuzzy Control System Design and Analysis: A Linear Matrix Inequality Approach[END_REF]. Each subsystem can be active for more than one sample. The corresponding graph representation is illustrated in Figure 1. The graph is G = {V, E}, with V = {1, 2, 3, 4} and E = {(1, 1), [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF][START_REF] Tanaka | Fuzzy Control System Design and Analysis: A Linear Matrix Inequality Approach[END_REF], (2, 1), (2, 2), (2, 3), (3, 1), [START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF][START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF], (4, 2), (4, 4)} A path P(1, 3) is given by P(1, 3) = [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched Lyapunov function approach[END_REF][START_REF] Tanaka | Fuzzy Control System Design and Analysis: A Linear Matrix Inequality Approach[END_REF][START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF][START_REF] Tanaka | Fuzzy regulators and fuzzy observers: relaxed stability conditions and LMI-based designs[END_REF]. The length of the path

P(1, 3) = [1, 4, 2, 3] is 3.
Once activated, a subsystem may be active continuously for at least p m i ∈ N + and at most p M i ∈ N + samples, that are assumed known. Our goal is to design a switching observer such that the estimation error dynamics converge to zero, with any admissible switching law.

0 and I denote the zero and identity matrices of appropriate dimensions, and a ( * ) denotes the term induced by symmetry. The subscript z + m (as in A 1,z+m ) stands for the scheduling vector being evaluated at the current sample plus mth instant, i.e., z(k + m).

In what follows, we will make use of the following results:

Lemma 1. [23] Consider a vector x ∈ R nx and two matrices Q = Q T ∈ R nx×nx and R ∈ R m×nx such that rank(R) < n x .
The two following expressions are equivalent:

1) x T Qx < 0, x ∈ {x ∈ R nx , x = 0, Rx = 0} 2) ∃M ∈ R nx×m such that Q + M R + R T M T < 0
Analysis and design for TS models often lead to double-sum negativity problems of the form

x T r i=1 r j=1 h i (z k )h j (z k )Γ ij x < 0 (2)
where Γ ij , i, j = 1, 2, . . . , r are matrices of appropriate dimensions.

Lemma 2. [START_REF] Wang | An approach to fuzzy control of nonlinear systems: stability and design issues[END_REF] The double-sum (2) is negative, if

Γ ii < 0 Γ ij + Γ ji < 0, i, j = 1, 2, . . . , r, i < j Lemma 3. [25] The double-sum (2) is negative, if Γ ii < 0 2 r -1 Γ ii + Γ ij + Γ ji < 0, i, j = 1, 2, . . . , r, i = j

III. OBSERVER DESIGN

Lets us now consider observer design for the switching system (1), repeated here for convenience:

x k+1 = A j,z x k + B j,z u k y k = C j,z x k
We propose to use a switching observer of the form

x k+1 = A j,z x k + B j,z u k + H -1 j,z L j,z (y k -y k ) y k = C j,z x k (3)
for the j-th subsystem, the observer switching together with the observed subsystem. The matrices H j,i and L j,i , j = 1, 2 . . . , n s , i = 1, 2, . . . , r are to be determined.

The error dynamics e k = x kx k using this observer, under the assumption that the scheduling variables are available online at sample k, can be written as

e k+1 = A j,z e k -H -1 j,z L j,z (y k -y k ) = (A j,z -H -1 j,z L j,z C j,z )e k (4)
which in itself is a switching system.

To derive the observer design conditions, meaning that the error dynamics (4) should be asymptotically stable, consider the switching Lyapunov function

V (e k ) = e T k P m,j,z e k (5)
defined during the switches, i.e., on the edges of the associated graph G = {V, E}, with (v m , v j ) ∈ E.

Remark: If a subsystem j may be active for several number of samples, the edge (v j , v j ) is also considered.

With this, the following result can be formulated:

Theorem 1. The error dynamics (4) is asymptotically stable, if there exist P m,j,k = P T m,j,k > 0, H j,k , (v m , v j ) ∈ E, (v j , v l ) ∈ E, k = 1, 2, . . . , r, such that -P m,j,z ( * ) H j,z A j,z -L j,z C j,z -H j,z -H T j,z + P j,l,z+1 < 0 (6)
for all admissible paths where [v m , v j , v l ] is an admissible path.

P(v m , v l ) = [v m , v j , v l ], v m ∈ V.
During the transition for j to l, the dynamics of the error system are described by

A j,z -H -1 j,z L j,z C j,z -I e k e k+1 = 0
Using Lemma 1, the difference in the Lyapunov function is negative, if there exists M such that -P m,j,z 0 0 P j,l,z+1

+ M A j,z -H -1 j,z L j,z C j,z -I + ( * ) < 0 By choosing M = 0 H j,z
we have directly [START_REF] Guerra | LMI-based relaxed nonquadratic stabilization conditions for nonlinear systems in the Takagi-Sugeno's form[END_REF].

Using Lemma 3 LMI conditions can be formulated, as follows:

Corollary 1. The error dynamics (4) is asymptotically stable, if there exist P m,j,k = P T m,j,k > 0, H j,k , (v m , v j ) ∈ E, k, l = 1, 2, . . . , r, such that

Γ m,j,l,γ kk < 0 2 r -1 Γ m,j,l,γ kk + Γ i,j,l,γ kβ + Γ i,j,l,γ βk < 0 k, β, γ = 1, 2, . . . , r with Γ m,j,l,γ kβ = -P m,j,k ( * ) H j,k A j,β -L j,k C j,β -H j,k + ( * ) + P j,l,γ for all admissible paths P(v m , v l ) = [v m , v j , v l ], v j ∈ V.

Remark:

In order to reduce the conservativeness by exploiting the knowledge available of the switching sequence, one can also use the observer

x k+1 = A j,z x k + B j,z u k + H -1 m,j,z L m,j,z (y k -y k ) y k = C j,z x k (7)
for the j-th subsystem, if the last switch has been from vertex m to vertex j.

The error dynamics e k = x kx k using this observer, can be written as

e k+1 = A j,z e k -H -1 m,j,z L m,j,z (y k -y k ) = (A j,z -H -1 m,j,z L m,j,z C j,z )e k (8)
Following the same steps as above, and in Finsler's lemma choosing

M = 0 H m,j,z
we have Corollary 2. The error dynamics (8) is asymptotically stable, if there exist P m,j,k

= P T m,j,k > 0, H m,j,k , (v m , v j ) ∈ E, (v j , v l ) ∈ E, k = 1, 2, . . . , r, such that -P m,j,z ( * ) H m,j,z A j,z -L m,j,z C j,z -H m,j,z + ( * ) + P j,l,z+1 < 0 (9) for all admissible paths P(v m , v l ) = [v m , v j , v l ], v m ∈ V.
The result above can be further extended to take into account more previous switches. However, this comes with the added computational cost and will eventually lead to considering all possible switching trajectories. To avoid this, but still reduce the conservativeness of the approach, let us now consider the α-sample variation of the Lyapunov function. As proven by [START_REF] Kruszewski | Stabilization of a class of nonlinear model with periodic parameters in the Takagi-Sugeno form[END_REF], a system is asymptotically stable, if the associated Lyapunov function decreases every α samples, α ≥ 1, instead of every sample. Thus, for observer design, let us consider the error dynamics (4) and the Lyapunov function [START_REF] Sala | Perspectives of fuzzy systems and control[END_REF]. Then, the following result can be formulated: Theorem 2. The error dynamics (4) is asymptotically stable, if there exist α ∈ N + , P i,j,k = P T i,j,k > 0, H j,k , (v i , v j ) ∈ E, k = 1, 2, . . . , r, such that

      -P v0,v1,z (
      < 0 (10) for all admissible paths P(v 0 , v α+1 ) = [v 0 , v 1 , . . . , v α+1 ],
where

Ω i = H vi,z+i-1 A vi,z+i-1 -L vi,z+i-1 C vi,z+i-1 .
Proof. Consider the switching Lyapunov function [START_REF] Sala | Perspectives of fuzzy systems and control[END_REF], defined on the edges of the associated graph, with P vi,vj ,z being active during the transition from vertex i to vertex j.

e k+α      = 0 with G i = A vi,z+i-1 -H -1 vi,z+i-1 L vi,z+i-1 C vi,
     + M      G 1 -I 0 . . . 0 0 G 2 -I . . . 0 . . . . . . . . . . . . 0 0 0 . . . -I      + ( * ) < 0
with G i defined as above.

By choosing

M =        0 0 . . . 0 H v1,z 0 . . . 0 0 H v2,z+1 . . . 0 . . . . . . . . . . . . 0 0 . . . H vα,z+α-1       
we have directly [START_REF] Dong | Dynamic output feedback H∞ control synthesis for discrete-time T-S fuzzy systems via switching fuzzy controllers[END_REF].

IV. EXAMPLE AND DISCUSSION

Let us first illustrate the conditions of Theorem 1 on an example.

Example 2. Consider the switching system -actually a periodic switching system -illustrated in Figure 2. Assuming that none of the subsystems may be active for more than one sample, i.e., p m = p M = 1, the graph is G = {V, E}, with V = {1, 2, 3} and

E = {(1, 2), (2, 3), (3, 1)}
. Consider the following local models of the TS system above:

A 1,1 = 0.61 0.10 1.82 0.50 A 1,2 = 0.21 1.59 0.34 0.77 Note that there are no observer gains L 2,1 and L 2,2 . This is because the second subsystem on its own is not observable. However, this observer is able to estimate the states of the switching system above. A trajectory of the error dynamics for the case when x(0) = [-1, 1] T and x(0) = 0 is illustrated in Figure 3. For the simulation, the membership functions used were h 1 = 1 2 (1 -sin x 1 ), h 2 = 1h 1 , and for L 2,1 and L 2,2 zero matrices were used. The initial subsystem was the first one, from which the system switched according to the periodic law.

C 1,1 = 1 0 C 1,2 = 1 0 A 2,1 = A 2,2 = 1.1 0 0.2 0.8 C 2,1 = 0 0 C 2,2 = 0 0 A 3,1 = 0.
Let us now discuss the conditions developed for the αsample variation of the Lyapunov function. In the example above, the conditions required that the Lyapunov function decreases with every switch/every sample, i.e., for all (v i , v j ) ∈ 1 The values have been truncated to two decimal places. E. That is, we had the conditions:

-P 1,2,z ( * ) H 2,z A 2,z -L 2,z C 2,z -H 2,z + ( * ) + P 2,3,z+1 < 0 -P 2,3,z ( * ) H 3,z A 3,z -L 3,z C 3,z -H 3,z + ( * ) + P 3,1,z+1 < 0 -P 3,1,z ( * ) H 1,z A 1,z -L 1,z C 1,z -H 1,z + ( * ) + P 1,2,z+1
< 0

On the other hand, a 2-sample variation means that the Lyapunov function has to decrease along paths of length 2, i.e., we have the conditions:

0 B B B @ -P1,2,z ( * ) ( * ) " H2,zA2,z -L2,zC2,z « -H2,z + ( * ) ( * ) 0 " H3,z+1A3,z+1 -L3,z+1C3,z+1 « " -H3,z+1 + ( * ) +P3,1,z+2 « 1 C C C A < 0 0 B B B @ -P2,3,z ( * ) ( * ) " H3,zA3,z -L3,zC3,z « -H3,z + ( * ) ( * ) 0 " H1,z+1A1,z+1 -L1,z+1C1,z+1 « " -H1,z+1 + ( * ) +P1,2,z+2 « 1 C C C A < 0 0 B B B @ -P3,1,z ( * ) ( * ) " H1,zA1,z -L1,zC1,z « -H1,z + ( * ) ( * ) 0 " H2,z+1A2,z+1 -L2,z+1C2,z+1 « " -H2,z+1 + ( * ) +P2,3,z+2 « 1 C C C A < 0
Furthermore, a 3-sample variation means that the Lyapunov function has to decrease along paths of length 3, which, in this case, is equivalent to the whole period of switching.

In the example above, one of the reasons for which the LMIs are feasible is that the subsystems may be active for a finite number of samples. In the case considered above, the switching occurs at every sample, i.e., each subsystem is active for one sample. Feasible LMIs are also obtained if the subsystems may be active for at most 2 samples. However, assuming that each subsystem, once activated may remain active for an infinite number of samples, the corresponding LMIs become unfeasible.

We assumed that the switching sequence is not known in advance and it cannot be directly influenced. Due to this assumption the conditions require that the estimation error dynamics is stable if the Lyapunov function decreases along every path of length α. This is the worst-case assumption, i.e., all possible combinations on switches between the estimation error subsystems have to be taken into account. If the switching sequence can be chosen, or it is known in advance, the conditions can be relaxed.

A shortcoming of the proposed conditions is the computational complexity of generating all the switching paths of length α, in particular for large α and large-scale switching systems and in consequence, a large number of LMIs that has to be solved. However, reducing the conservativeness of the conditions by introducing additional sums in the Lyapunov function increases the number of LMIs.

V. CONCLUSIONS

In this paper we developed observer design conditions for discrete-time switching systems represented by Takagi-Sugeno fuzzy models, by using a switching Lyapunov function defined on the switches. We assumed that the switching sequence is not known in advance and it cannot be directly influenced. The developed conditions have been formulated as LMIs.

Fig. 1 .

 1 Fig. 1. Graph representation of the switching system in Example 1.

Proof.

  Consider the switching Lyapunov function[START_REF] Sala | Perspectives of fuzzy systems and control[END_REF], defined on the edges of the associated graph, with e T k P m,j,z e k being active during the transition from vertex m to vertex j. The difference in the Lyapunov function for two consecutive samples is ∆V = e T k+1 P j,l,z+1 e k+1e T k P m,j,z e k

  The difference in the Lyapunov function for α consecutive samples is ∆V = e T k+α P vα,vα+1,z+α e k+αe T k P v0,v1,z e k = e k e k+α T -P v0,v1,z 0 0 P vα,vα+1,z+α • e k e k+α where [v 0 , v 1 , . . . , v α+1 ] is an admissible path.Along the switching sequence [v 0 , v 1 , . . . , v α+1 ], the error dynamics are described by 

1 2 3 Fig. 2 .

 32 Fig. 2. Periodic switching system for Example 2.

Fig. 3 .

 3 Fig. 3. Estimation error for Example 2.

ACKNOWLEDGMENT

  This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS -UEFISCDI, project number PN-II-RU-TE-2011-3-0043, contract number 74/05.10.2011, by International Campus on Safety and Intermodality in Transportation the European Community, the Délegation Régionale à la Recherche et à la Technologie, the Ministére de L'Enseignement supérieur et de la Recherche the Region Nord Pas de Calais and the Centre Nationale de la Recherche Scientifique: the authors gratefully acknowledge the support of these institutions.

  The second subsystem is linear, but it is unstable and unobservable. The second local model of the third subsystem is again unstable and unobservable. Due to this, methods available in the literature yield unfeasible LMIs. However, using Lemma 2 to formulate LMI conditions for Theorem 1, we obtain1

		11 0.09 0.09 0.07	A 3,2 =	0.32 0.30 0.09 0.38
	C 3,1 = 0 1	C 3,2 = 0 0
	P 1,2,1 =	4.69 -3.97 -3.97 7.79	P 1,2,2 =	4.69 -3.97 -3.97 7.79
	P 2,3,1 =	3.45 -3.35 -3.35 5.05	P 2,3,2 =	4.09 -3.91 -3.91 6.41
	P 3,1,1 =	3.48 -3.97 -3.97 5.18	P 3,1,2 =	4.74 -5.67 -5.67 7.39
	H 1,1 =	4.83 -2.45 -2.50 6.26	H 1,2 =	3.77 -3.21 -1.84 6.46
	H 2,1 =	3.57 -3.21 -2.54 5.35	H 2,2 =	3.57 -3.21 -2.54 5.35
	H 3,1 =	4.55 -5.05 0.45 5.84	H 3,2 =	4.14 -4.61 -4.85 6.27
	L 1,1 =	-2.14 12.08	L 1,2 =	2.40 3.45
	L 3,1 =	2.42 4.27	L 3,2 =	3.89 7.93


Christoph Lüders

Eléonore Bellot

François Fages

Ovidiu Radulescu

Sylvain Soliman

Symbolic Methods for Biological Networks D2.1 Report on Scalable Methods for Tropical Solutions (T1.2)

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Tropical geometry is a natural approach to find the scalings necessary for slow/fast decompositions and model order reduction of biochemical reaction networks [SGF + 15]. These reductions extending the quasi-steady state and quasi-equilibrium approximations, well known in biochemistry, are based on automated algorithms. The first step in these algorithms consists in computing tropical equilibration solutions.

The basic scaling idea is to consider that both variables x i ∈ R + and parameters k i ∈ R + of the ODE system describing the chemical kinetics are functions of a small, positive parameter. The dominant terms of these functions are considered to be powers of ε. Thus, x i = xi ε a i and k i = ki ε γ i . For biochemical reaction networks, the reaction rates are rather generally rational functions of the variables x i , therefore the chemical kinetic equations are rational or polynomial ODEs. Using the variable and parameter scaling, all the monomials terms in these ODEs can be also expressed as powers of ε. Furthermore, the order of the dominant terms in each ODE provide the timescale order of each variable. A scaling satisfies the tropical equilibration conditions when the dominant monomial terms can be compensated, i.e. when not all the dominant monomials have the same sign. A tropical equilibration is complete when the above condition is satisfied for all the ODEs and partial when the above condition applies at least to the ODEs describing the fast variables dynamics. In the latter case, the timescale condition on the variables that are equilibrated has to be solved consistently with the monomial order condition.

In this report, we present three methods and study their scalability to solve complete equilibration problems in the BioModels model repository. The first two methods, a naive polyhedral method using PtCut [START_REF] Lüders | PtCut: Calculate Tropical Equilibrations and Prevarieties[END_REF], and a Satisfiability Modulo Theories (SMT) method recently introduced in [Lüd20a] using the SMT solver CVC4 [BCD + 11], compute the set of solutions over the real numbers. Satisfiability Modulo Theories (SMT) is usually build on top of SAT (Boolean satisfiability), which was the first problem that was proved, in the form of 3SAT, to be NP-complete. SMT allows to test a logical formula with unknowns and relations for satisfiability and, if it is so, for an assignment of the unknowns that leads to the formula's satisfiability [START_REF] Monniaux | A survey of satisfiability modulo theory[END_REF]. SMT checking is used today in verification of computer hardware and software and has advanced much in recent years due to advances in technology and industrial applications [START_REF] De | Satisfiability modulo theories: introduction and applications[END_REF]. The SMT approach introduced in [Lüd20a] and used here is faster than the polyhedral approach for models that would otherwise take more than one minute to compute, in many cases by a factor of 25 or more. Furthermore, this method is an anytime algorithm, thus offering a way to compute parts of the solution when the polyhedral approach is infeasible.

The third method evaluated on the same benchmark here, is the Constraint Programming (CP) method presented in [START_REF] Soliman | A constraint solving approach to model reduction by tropical equilibration[END_REF] and implemented in Biocham-4. This method searches for integer solutions. The CP approach is similarly efficient as the SMT method, generally under two minutes for the polynomial and rational fractional ODE systems of this benchmark. Interestingly, this method reveals that a large part of the models that can be equilibrated over the reals have in fact no integer solution.

Both SMT and CP methods thus appear to be scalable on BioModels.

The rest of the report is organized as follows. The next section describes the idea of tropical geometry and tropical equilibrations in ODE systems. Section 3 describes the SMT approach with some details. Section 4 describes the CP approach for searching for integer solutions. Section 5 presents the benchmark of polynomial and fractional ODE systems obtained from the curated part of the BioModels model repository. Then in Section 6, we present the results and compare the performance of the three methods on this benchmark. We conclude in the last section on the scalability of SMT and CP solvers on real and integer tropical equilibration problems for real-size problems and on some possible improvements for future work.

Tropical Equilibration Problems

Let x ∈ R d + and α, β ∈ N m 0 be multi-indices with x α = i x α i i .
We express a polynomial function as

P (x) = α∈P k α x α - β∈N k β x β ,
where P and N are index sets of the monomials with positive resp. negative sign.

The order of a monomial k α x α is given by the power of ε in that monomial and can be straightforwardly calculated as γ α + i α i a i , where γ α is the order of k α and a i are the orders of x i .

When → 0 dominance relation between terms is given by their orders, typically ε a << ε b iff a > b. As a consequence, the dominant term is the term of minimal order, and orders are computed according to the rules of min-plus algebra, where multiplication becomes plus and addition becomes min.

The critical observation is that during dynamics on a slow invariant manifold, dominant terms can not remain uncompensated. For the polynomial P (x) this means that

min α∈P γ α + i α i a i = min β∈N γ β + i β i a i , (1)
called tropical equilibration condition.

A total tropical equilibration solution satisfies (1) for all the polynomial ODEs or for all numerators for rational ODEs in a common denominator representation.

Observe that (1) comprises now only of minima of linear functions and the set of solutions can thus readily be expressed as a set of polyhedra in the space of variable orders (a 1 , a 2 , . . . , a n ).

One polyhedron is defined by each combination of α ∈ P and β ∈ N that yields a hyperplane via

γ α + i α i a i = γ β + i β i x i , (2)
while for all η ∈ P ∪ N half-spaces are defined by

γ α + i α i a i ≤ γ η + i η i a i . ( 3
)
The whole set of polyhedra is defined by cycling over all possible choices for α and β.

Since the ODE system often consists of multiple equations, we get multiple sets of polyhedra. Since all l.h.s. have to be zero, the constraints that define the set of polyhedra all have to hold at the same time, i.e., the sets of polyhedra have to be intersected. The resulting set of polyhedra is called the set of tropical equilibration solutions.

Solving Tropical Equilibrations Over the Real Numbers

Satisfyability Modulo Theories (SMT) Method

We are interested in the intersection of several sets of polyhedra. That is, given polyhedra P ij and sets of polyhedra B i = {P i1 , P i2 , . . . , P in i }, we are interested in their intersection

i B i = {R 1 , R 2 , . . . , R },
where the R i are again polyhedra. In this article we first show how to use SMT checking to solve this question. Satisfiability Modulo Theories (SMT) allows us to decide if a logical formula, with atoms that are themselves equations or inequalities, is satisfiable or not. For example, x > 1 ∧ x < 2 is an SMT formula. One has to specify a theory of numbers that unknowns in the formula can assume. In the above example, the problem is satisfiable in the theory of real numbers, but not in the theory of integers. If an SMT problem is satisfiable, SMT can return a model, which is an assignment for all unknowns in the formula.

SMT solvers may accept logical formulas in general, but often formulas are written in conjunctive normal form (CNF), i.e., as a number of AND clauses called assertions. If SMT solvers are used in incremental mode, one can, after they have found a solution, add additional assertions and "continue" to look for further solutions, but keep using their internal search heuristics. We make use of that later.

Polyhedral Representation of the Set of Solutions

A polyhedron is defined as the intersection of finitely many hyperplanes and half-spaces. Furthermore, each hyperplane can be expressed as two (closed) half-spaces, thus a polyhedron can be described as a finite number of half-spaces [START_REF] Günter | Lectures on Polytopes[END_REF].

Since we are using computers to do the work, we represent numbers ∈ R as elements of Q for reasons of numerical stability. All numbers we are dealing with are either provided with finite precision or can be computed to arbitrary precision.

Hence, given ambient space Q d , a (closed) half-space H is the set

H = {x ∈ Q d | λ 0 + λ 1 x 1 + λ 2 x 2 + . . . + λ d x d ≤ 0, λ i ∈ Q}. (4)
Given half-spaces H k , we define a polyhedron as

P = k H k . ( 5
)
A bag is what we call a set of polyhedra P j and describe it as a union, i.e.,

B = j P j .
Finally, we are looking for the intersection of said bags B i , that is

V = i B i = i j P ij = k R k . ( 6
)
The naive solution to the problem of computing the intersection is to cycle successively through all combinations. To do that, pick two bags B j and B j , j = j , and intersect all polyhedra from one with all polyhedra from the other to form a new bag B . Then, remove B j and B j from the set of bags and insert B instead. Continue this procedure until there is only one bag left, which will then consist of the sought polyhedra R 1 , R 2 , . . . , R . That is the solution that was used in [SGF + 15] and, with some refinements, in PtCut [START_REF] Lüders | PtCut: Calculate Tropical Equilibrations and Prevarieties[END_REF].

The problem with this solution is that the complexity is exponential in the number of bags.

In practice, it often happens that the number of intermediate results, i.e., the number of polyhedra in some B , can be very high, even if in the end there are only a few solution polyhedra. This intermediate expression swell makes computing the intersection for some models infeasible.

Table 4 shows BioModels from our survey, their number of resulting polyhedra and their maximum number of intermediate polyhedra. More details on the computation can be found in Section 5.

The SMT Procedure

First, we have to formulate our problem as an SMT problem. Fortunately, it is easy to convert a polyhedron as defined in (5) into a logical formula. Set theory maps easily to logical formulas with union mapping to logical OR and intersection to logical AND. In the following, H denotes the logical formula that defines the set H.

Thus, (6) expands to Ṽ = i j Hijk (7)

and definition (4) of H ijk employs a linear function that can be used as a formula in SMT. Call the resulting SMT formula f . Thus, we can use an SMT solver to decide the satisfiability of formula f and, what's more important, get a point x ∈ Q d that satisfies the constraints.

Next, we are looking for a matching polyhedron that includes point x and is included in the solution V . Since point x is contained in the intersection of the B i , it must be contained in at least one polyhedron P ij per bag B i . Thus, we cycle through all B i to find a containing P ij , call it P i . (There may be more than one P ij x, but any will do.)

Obviously, the intersection R = i P i includes point x, but most likely R has higher dimension than that. Furthermore, since R is the intersection of exactly one polyhedron per bag it is included in V as well. Hence, we have found a whole polyhedron that includes point x.

In the next step, we modify our initial formula f to exclude the polyhedron R , like this:

f = f ∧ ¬R . ( 8
)
Note that we are only adding another assertion to the formula, so we can utilize the incremental mode of SMT solvers to save (a lot of!) time for its next computation.

The important observation here is that we are expanding the original formula f -which describes all solution points-to exclude what we already know to be a solution and continue the search. We can iterate this process until formula f is unsatisfiable. This is the algorithm in Python-style pseudocode:

# input : a list ' ll ' of sets of p o l y h e d r a . # output : a list ' rr ' of p o l y h e d r a . def c o m p u t e _ p r e v a r i e t y ( ll ): # set the solver to re -use its h e u r i s t i c s solver = Solver ( incremental = True ) f = c o n v e r t _ t o _ S M T _ f o r m u l a ( ll ) rr = [] # results list while True :

# add the formula to the ( e x i s t i n g ) a s s e r t i o n s . solver . add_assertion ( f ) # get the model ( a v a r i a b l e a s s i g n m e n t ) that fits the # constraints , or None if 'f ' is u n s a t i s f i a b l e .

x = solver . get_model ()

if x is None : break R = [] # r

e s u l t i n g p o l y h e d r o n # cycle through all bags 'B ' and # collect c o n s t r a i n t s of p o l y h e d r o n c o n t a i n i n g 'x '. for B in ll :

# cycle through all p o l y h e d r a 'P ' in bag 'B '. for P in B : if P . contains ( x ): R . append ( P . constraints ())

break # now 'R ' defines a p o l y h e d r o n surely in the i n t e r s e c t i o n . # exclude 'R ' from further s e a r c h e s . f = Not ( R ) # new a s s e r t i o n for next round rr . append ( R ) return rr # list of p o l y h e d r a .

The result of this function is a list of polyhedra. Mathematically, the set of polyhedra describes the equilibrium (resp. prevariety) V . Yet, there are some problems that we address in the next section.

The logic used for SMT formulas is QF_LRA, that is, quantifier-free linear real arithmetic ("real" here means rational). This allows Boolean propositional logic of equations/inequalities consisting of linear polynomials over elements of Q [BST + 10].

Improvements to the SMT Procedure

Non-maximal Polyhedra

The main issue we experience with the procedure compute_prevariety() is that the polyhedron R computed from point x is often not maximal. That is, R is only a lower dimensional face of a higher-dimensional polyhedron. The full high-dimensional polyhedron will eventually be found by the procedure, but earlier-found lower-dimensional faces would still remain in the result list rr, albeit superfluous.

To avoid this, we test if each newly found polyhedron R is included in one of the already computed polyhedra of result list rr. Unfortunately, this causes quadratic run-time in the number of resulting polyhedra. But there is an observation that can reduce the time needed.

If the newly found polyhedron R is included in some polyhedron R , then obviously, point x ∈ R is included in R as well. Testing if a point is included in a polyhedron is simple and fast, so one can test this first and only if this test succeeds one must perform the full polyhedron inclusion test. Measurements show that with this heuristic, almost all polyhedron inclusion tests can be avoided. See Section 5 for details.

In our procedure, we would have to modify function contains() in line 22 and function append() in line 23 according to these observations.

Redundant Constraints

Another issue is the redundancy of the constraints that are collected in line 23. Efficiency can be increased by minimizing the set of constraints: the larger the number of constraints, the larger the memory demand and, of course, SMT search times.

One can simply cycle through all constraints, test if each of them is really required and if not, drop it. The remaining set is not necessarily a minimal set, though.

Here's how this can be done: Let c be the constraint in consideration, g the formula before and g = g ∧ c after the addition. If g is more restrictive than g (i.e., c makes a difference), then the following is unsatisfiable:

g ∧ ¬g = g ∧ ¬(g ∧ c) = g ∧ (¬g ∨ ¬c) = (g ∧ ¬g) ∨ (g ∧ ¬c) = g ∧ ¬c.
That formula can be applied by SMT to all constraints to drop superfluous ones.

Preprocessing

Different possibilities were explored to improve the speed of the procedure by preprocessing the input, i.e., the sets of polyhedra.

For one, one can collect all constraints from all bags with only one polyhedron each. Call the resulting polyhedron C. Because of distributivity these constraints hold for all polyhedra of the equilibrium. Hence, we can intersect all polyhedra in their bags with C to test if the intersection is empty, in which case we can drop the polyhedron from its bag and thus reduce problem complexity.

A more potent version of this technique can be used to test polyhedra in all bags on if they are required for the definition of the equilibrium. Let B be the polyhedron in question, A the union of all other polyhedra in B's bag and C the intersection of all other bags. Then the equilibrium is

(A ∪ B) ∩ C. If B is required, then A ∩ C = (A ∪ B) ∩ C and in particular A ∩ C (A ∪ B) ∩ C. Thus, the following set is non-empty: ((A ∪ B) ∩ C) \ (A ∩ C) = (A ∪ B) ∩ C ∩ A ∩ C = (A ∪ B) ∩ C ∩ (A ∪ C) = (A ∪ B) ∩ C ∩ A = B ∩ C ∩ A.
This can easily be tested with SMT for each polyhedron B per bag and the superfluous polyhedra can be dropped, again reducing problem complexity.

Solving Tropical Equilibrations Over the Integers

Constraint-based Modeling

The constraint programming (CP) approach is a numerical approach that encodes the tropical equilibration problem as a Constraint Satisfaction Problem (CSP) [START_REF] Soliman | A constraint solving approach to model reduction by tropical equilibration[END_REF]. This method is implemented in Biocham, which is itself implemented in SWI-Prolog and uses a constraint-solving library over integers.

It doesn't work with a symbolic polyhedral representation of the solutions, but searches for numerical solutions that are points in the intersection of the sets of polyhedra explained earlier. Both complete equilibrations (i.e., for all variables) and partial equilibrations (i.e., for some fast variables) can be computed. In this report and in the following evaluation benchmark, only complete equlibration problems are considered.

Constraint programming, as SMT, deals with combinatorial optimization problems, but it uses constraints in an active manner during search to prune the search space, instead of a symbolic manner to perform clause resolution.

A constraint-based model is defined by a finite set of variables given with their domain, and a finite set of constraints.

• Variables. Here the variables are the orders of magnitude a i ∈ Z used to rescale the system by posting x i = xi ε a i , as seen in Section 2. Taking integer values ensure that these orders of magnitude are well separated.

• Constraints.

-Tropical equilibration constraints. The degree in ε of the monomials in P and in N are integer linear expressions in a i , the variables of the problem. Such linear expressions are used to filter the domain of the variables by computing their minimum and maximum values according to the bounds of the domain of the variables. Tropical equilibration needs to enforce that the minimum degrees in P and in N are equal (1).

-The constraint of minimum, A = min(B, C), is implemented with simple inequality constraints A ≤ B, A ≤ C, plus reified constraints. A reified constraint is a constraint that links bidirectionnaly the value of one boolean variable X to the satisfaction of another constraint, say on two variables Y and Z, for instance

X ⇔ Y = Z to state that X is true if and only if Y = Z. Two reified constraints, (X ⇔ A = B) ∧ (Y ⇔ A = C) ∧ X + Y ≤ 1 ∧ A ≤ B ∧ A ≤ C,
are used to define the minimum constraint to enforce that A is equal to B or C, in addition that A is greater than B and C.

-Linear invariants. A conservation law is a linear invariant, i.e., an equality beetwen a linear combination of the x i and a constant. Having the initial values of x i , and thus the value of the constant, one can add a constraint on the a i so that the equality of the degrees in the rescaled invariant equation still holds. This allows us to get tropical equilibrations that are compatible with the initial system.

Constraint-based Search

Constraints come with domain filtering algorithms which dynamically prune the domain of variables when the domain of other variables change in a constraint. However, this strategy is not sufficient to decide the satisfiability of a system of constraints and must be combined with a search procedure for virtually enumerating all possible values of the variables, or more generally, for enumerating all the alternatives of disjunctive constraints.

For solving tropical equilibration problems over finite domains, we obtained good results with dichotomic search by bisecting the domain of the variables (bisect option in SWI-Prolog) without any particular heuristics for choosing the variables. Finite bounds must be given on the values of the variables, but this is not a restriction in practice, since biochemical species' concentrations usually do not vary by more than a hundred of magnitude orders. Furthermore, in order to speed-up the computation of all solutions in such large domains, we used an iterative domain expansion strategy: the problem is first tried with a domain of [ε -2 , ε 2 ] for all variables , i.e., equilibrations are searched by rescaling in the [10 -2 , 10 2 ] interval. If that fails, the domain is doubled and the problem tried again until the limit interval [ε -256 , ε 256 ] is reached.

Benchmark of Tropical Equilibration Problems from BioModels

This benchmark uses models from the BioModels repository [LNBB + 06]. Since BioModels are written in SBML, and to avoid inconsistency in parsing, ODEparse [START_REF] Lüders | ODEparse: generate ODEs from SBML[END_REF] was used to convert SBML into ODEs. The resulting ODE systems can be downloaded from ODEbase [LEN + 19] at http://odebase.cs.uni-bonn.de. A parser to read files from ODEbase was implemented in Biocham. Features of SBML models like events, constraints, time dependency, multiple compartments, or function definitions, were ignored.

The ODE systems selected contain only polynomials or rational functions. For the latter, each equation is multiplied with its principal denominator to get a polynomial. The polynomials are then tropicalized as sketched in Section 2.

The parameter ε was set to 1/11 for all methods.

Although the three methods can handle partial equilibration problems, we consider here only the problem of complete equilibration, i.e., when all variables equilibrate.

SMT Method Parameters

The computations work for both equilibria and prevarietyies, since the difference is only in the input. Here, run-times only for computation of equilibria will be shown.

Furthermore, for computation, the logarithms log ε (k α ) in (2) & (3) were rounded to integers or rationals for reasons of numerical stability.

The sets of polyhedra created by tropicalization were saved with polyhedra expressed as sets of equalities and inequalities.

The software used for polyhedral computation was PtCut 3.5.1 [START_REF] Lüders | PtCut: Calculate Tropical Equilibrations and Prevarieties[END_REF] as available from the web site of the first author. Polyhedral computation were done with the help of the Parma Polyhedra Library (PPL) [START_REF] Bagnara | The Parma Polyhedra Library: Toward a Complete Set of Numerical Abstractions for the Analysis and Verification of Hardware and Software Systems[END_REF], version 1.2.

For SMT computation, SMTcut 4.6.0 and the PySMT framework [START_REF] Gario | PySMT: a solver-agnostic library for fast prototyping of SMT-based algorithms[END_REF], version 0.8.0 with SMT engines MathSAT 5.5.1 [CGSS13] and Z3 4.8.4 [START_REF] De | Z3: An efficient SMT solver[END_REF] were used.

Tests were run with Python 3.7.7 on an Intel Core i7-5820K CPU with 48 GB of memory under Windows 10 64-bit. Neither PtCut nor SMTcut make active use of multithreading.

CP Method Parameters

Concerning the CP solver over the integers, the domain of the degree in epsilon was set to the set {0, 1, . . . , 256}. Enlarging this domain had no impact on the results obtained.

In this benchmark we did not include initial values of the variables.

The CP tropicalization was run on Biocham 4.4.14.

Availability of Code and Data

The used software program, SMTcut, is available under a free software licence from https://gitlab.com/cxxl/smtcut. The ODE systems that were used can be downloaded from ODEbase, http://odebase.cs.uni-bonn.de and the data that was used as input as well as output resides in one large (≈ 20 MiB) repository available under https: //gitlab.com/cxxl/smtcut_data_1.

The software Biocham is available at http://lifeware.inria.fr/biocham4/, and the data used at http://lifeware.inria.fr/wiki/Main/Software#DeliverableSYMBIONT

Performance Results

Table 1 shows the run-times for all methods on a subset of models that have complete tropical equilibration solutions over the real numbers, except for model 397 which has no solution. The run-times value our attention for certain models:

• BioModels 14, 151, 410, 560 & 730: each of them could not be computed with PtCut and it is likely infeasible to that with this approach. SMTcut was at least able to compute many polyhedra, even though it is unknown how large a part of the full equilibrium this constitutes. On the other hand, the CP method found integer solutions for model 14, and proved that models 151 and 560 are unsatisfiable over the integers in less than 3 seconds.

• BioModels 183 and 491, 492: Here SMTcut was able to play out its full potential: with only one polyhedron in the equilibrium, it took only some rounds until the maximal polyhedron was found. PtCut, on the other hand, terminated after a day's work with an intermediate number of 15000 polyhedra and still only 5 of 65 iterations done. The CP method failed to solve model 183 in less than 15 minutes, but found integer solutions for the other two models in 6 seconds.

• BioModel 397 doesn't have a solution at all. SMTcut realizes this in 0.2 s, whereas CP takes 1.2 s and PtCut almost 15 s.

• BioModels 48, 73, 74, 93, 105, 407, and especially 501: computation time of SMTcut was at least 20 times lower than with PtCut. The CP method was even faster by one or two orders of magnitude on models 74, 93, 105 and 407, but was in timeout on model 73.

• BioModels 102, 328, 430 & 498: here PtCut was at least 3 times faster than SMTcut. The reason might be that for models of medium dimension (< 20) and with many polyhedra, PtCut can be faster if the intermediate expression swell is not too large, see Table 4. All these models were solved by CP in less than a second or 7.3 s Over the reals, SMTcut was the better choice in all cases where PtCut needed more than 58 s of computation time. Over the integers, the CP method could be even faster to conclude by up to two orders of magnitude, but could also fail with a timeout of 2 minutes whereas SMTcut took a few seconds. The CP method reveals that nearly 30% of the tropical equilibrations satisfiable over R have no solution over Z.

Further Results with the CP Method

The complete equilibrations were computed and the existence of a tropical equilibrium or not was reported in the Tables 1, 2 and 3.

Results for all the polynomial models are presented in Table 2. On the over 150 models, only one had a timeout over 15 minutes of computation, less than 10% ran over 2 seconds, with one of them running over 20 seconds. The extreme majority of the models took then less than 2 seconds to be computed.

Rational functions models results, presented in Table 3, were way slower than original polynomial models (multiplication time not taken in account in the table) to tropicalize. The multiplication by the common denominator, even after simplification, leaves large equations to work with, and was involved in the slow-down of the computation. A timeout over 15 minutes happens for more than half the models, and often exceed more than an hour of computation. The computation is taking more than 20 seconds to compute in around one eight of the cases not timed-out.

Further Results for the SMT Method

Benchmarking of Different SMT Solvers

Under Windows, we could only use MathSAT [START_REF] Cimatti | The MathSAT5 SMT Solver[END_REF] and Z3 [START_REF] De | Z3: An efficient SMT solver[END_REF] and were unable to install the CVC4 [BCD + 11] and Yices [START_REF] Dutertre | Yices 2.2[END_REF] solvers. The other solvers that are in principle supported by PySMT are CUDD, PicoSAT and Boolector, but they don't support the QF_LRA logic, so they were no option.

Hence, we switched to Ubuntu Linux 18.04 64-bit and benchmarked some models again on the same machine to get a comparison. Under Ubuntu we could work with MathSAT 5.5.1, Z3 4.8.4 and Yices 2.6.0. Unfortunately, CVC4 1.7-prerelease could be installed, but didn't work properly.

Table 5 shows the run-times for selected models.

A side-by-side comparison of operating systems suggests that Windows and Ubuntu implementations are similar. MathSAT on Windows was some 3-9% slower than on Ubuntu, whereas Z3 was not always faster on any one operating system.

A comparison of MathSAT and Z3 solvers suggests that they yield comparable run-times, with Z3 in one case almost twice as fast as MathSAT (BioModel 73). Yices could not be evaluated properly, since it crashed with all larger models.

Preprocessing

Table 6 shows a comparison of times with and without preprocessing the input. Several observations can be made:

• BioModel 183: the process was killed after over 11 hours of preprocessing. This is due to the extraordinary number of polyhedra in some input bags.

• With the exception of BioModel 74 (and, of course, 183), all models with a run-time over 24 s ran faster with preprocessing (including the time for preprocessing itself) than without.

• All models with run-times less than 17 s ran slower (or as fast in the case of BioModels 22 & 498) with preprocessing than without.

So it seems that this kind of preprocessing is advisable only in special cases.

Profiling

The relative time needed in different parts of the procedure to compute the whole equilibrium varies with the number of resulting polyhedra. Table 7 shows the relative time spent in three different parts of the procedure:

• S: Searching for another point outside the already known polyhedra,

• M: Minimizing the found polyhedron,

• I: Inserting that polyhedron into the list of already known ones and testing for inclusion.

From the numbers it is obvious that the inclusion check time takes relatively more time as the number of polyhedra grows. Since that routine requires quadratic time, this is as expected.

Furthermore, the SMT checking for every next point needs more time as the number of polyhedra grows as well. This is due to the fact that the formula grows over time (to exclude all already found polyhedra) and thus the search gets more complicated, and we should not forget that we would expect an exponential run-time from theory.

The time needed for minimization, in contrast, is getting less prominent, which is no surprise, since the time required to check for superfluous constraints in a polyhedron does not depend on the number of already found polyhedra, nor do we expect to find polyhedra with more constraints as the search takes more time.

Polyhedral Inclusion Testing

As was described in Section 3.3, the test for inclusion of a newly found polyhedron in the set of already found ones can be sped up by testing containment of the included point x, which was found by the SMT solver.

Some cursory investigation shows that the number of full checks that still have to be done is about 0.26-0.68 times the number of polyhedra in the result list.

Yet, the main message is that even though the constant is low, the run-time of the inclusion test is still quadratic and it becomes the most dominant part of the computation for large result lists.

Conclusion

We have shown that tropical equilibration problems of real size, coming from the curated part of the BioModels model repository, can be solved quite efficiently both over the integers and over the real numbers.

For solving over the reals, a new method, based on SMT, was presented and compared to a known algorithm using purely polyhedral methods. The run-times were always smaller for problems that would otherwise take more than 58 seconds to compute, sometimes by a factor of 25 or more. The novel method has the further advantage that, even if computation of the entirety of the solution is infeasible, parts of it can be computed and more computation power would lead to more parts being computed. Furthermore, extensive benchmarks with different SMT solvers under different operating systems were run on tropicalizations of 46 different BioModels. The choice of operating systems doesn't make a significant difference. The choice of SMT solver sometimes makes a difference, but there is no clear winner there. Only MathSAT and Z3 worked under all circumstances and had comparable run-times.

Tropical solutions over the integers ensure that the orders of magnitude of the variables are well separated. We have shown that the Constraint Programming method introduced in [SFR14] exhibits good performances similarly to the SMT method on this benchmark, with run-times mostly inferior to 2 minutes per model. However, this evaluation in BioModels revealed that 30% of the tropical equilibration problems solvable over the real numbers have no integer solution.

Both the SMT and CP methods thus appear to be scalable on BioModels. More work is now needed to compare the real solutions with the integer solutions found by the two methods, and assess the scope of their respective use for model reduction applications.

A. Run-times

Table 1 :

 1 Model number in BioModels, number of variables, and comparison of run-times between PtCut and SMTcut on Windows with MathSAT solver on one side, and constraint programming (CP) in Biocham on the other side. Models marked with a star (*) have a rational vector field. Models with time in italic for CP have no integer tropical equilibration while they have a real-valued equilibration. Timeouts are indicated by a lower bound on the computation time.

	BM R	Dim Combos Polyhedra	PtCut [s]	SMTcut [s] Speed-up CP Biocham [s]
	14		86	10 67	> 6996 > 95594.9	> 79964.1	-	14.1
	21	*	10	10 8	46	3.0	1.4	2.1	> 900
	22	*	10	10 9	147	2.0	3.7	0.5	> 900
	26		11	10 6	6	0.1	0.3	0.5	0.06
	28		16	10 9	17	0.4	1.2	0.4	0.2
	30		18	10 10	6	0.3	0.5	0.8	0.2
	32		36	10 21	244	110.0	37.2	3.0	0.6
	41		10	10 14	4	1.4	0.3	4.8	> 900
	48	*	23	10 22	5	29.2	1.4	21.2	12.9

Table 2 :

 2 Satisfiability over the integers and run-times in seconds with CP in Biocham for the polynomial models of the curated part of BioModels.

	BM	output Biocham [t]	261	unsat	0.2
	2 5 9 11 14 26 28 30 35 38 40 50 52 57 60 69 72 80 82 85 86 91 92 99 102 103 105 108 123 147 150 156 159 163 173 178 197 198 199 200 205 220 226 229 230 233 243 257 259 260	sat sat sat unsat sat unsat sat sat unsat sat sat unsat unsat unsat sat unsat sat sat sat sat sat unsat unsat unsat sat sat unsat sat sat sat sat sat sat unsat unsat unsat unsat sat sat sat unsat unsat unsat unsat unsat unsat unsat sat unsat unsat	0.2 0.1 0.3 0.3 14.1 0.1 0.2 0.2 0.1 0.1 0.01 0.1 0.04 0.1 0.01 0.1 0.03 0.1 0.05 0.3 0.5 0.2 0.01 0.1 0.3 0.6 1.5 0.2 0.3 14.4 0.01 0.01 0.01 0.2 0.4 0.02 0.1 0.05 0.05 190.3 18.0 1.3 0.3 0.1 0.7 0.01 0.1 0.1 0.2 0.2	262 263 264 267 270 271 272 282 283 286 289 292 306 307 309 310 311 312 314 315 321 332 333 334 335 357 359 360 361 362 363 364 365 389 405 407 413 416 430 431 459 460 475 478 483 484 485 486 487 491 492	unsat unsat sat unsat sat unsat unsat unsat unsat unsat unsat unsat unsat unsat sat unsat unsat unsat unsat sat unsat unsat unsat unsat unsat unsat unsat unsat sat unsat unsat unsat sat unsat sat unsat unsat unsat sat sat unsat unsat unsat unsat sat unsat sat unsat unsat unsat unsat	0.1 0.1 0.1 0.01 0.4 0.04 0.04 0.01 0.01 2.2 0.04 0.01 0.02 0.01 0.01 0.01 0.01 0.01 0.05 0.4 0.02 2.7 1.2 2.4 0.4 0.04 0.1 0.05 0.04 0.5 0.01 0.1 0.5 0.2 0.1 2.0 0.02 0.7 0.4 0.4 0.01 0.01 0.2 0.5 0.1 0.002 0.01 0.003 0.02 6.3 6.2

Table 3 :

 3 Satisfiability over the integers and run-times in seconds with CP in Biocham for the rational fractional models of the curated part of BioModels.

	BM	output Biocham [t]
	10	timeout	> 900
	16	sat	175.4
	17	timeout	> 900
	21	timeout	> 900
	22	timeout	> 900
	23	timeout	> 900
	27	sat	0.2
	29	sat	1.1
	31	sat	482.4
	32	unsat	0.6
	33	unsat	1.1
	37	unsat	0.1
	39	sat	0.5
	41	timeout	> 900
	42	timeout	> 900
	43	sat	0.4
	44	timeout	> 900
	45	unsat	6.2

Table 4 :

 4 Intermediate expression swell for some BioModels computed with the SMT method.

	BM	Combos PtCut [s] Polyhedra Max. int. ph's
	21	10 8	2.978	46	3408
	22	10 9	2.036	147	1170
	32	10 21	109.952	244	1092
	41	10 14	1.416	4	924
	48	10 22	29.175	5	1160
	93	10 26	6113.411	596	47772
	102	10 10	7.199	322	4784
	103	10 18	763.480	1938	111402
	105	10 23	644.176	130	21088
	147	10 16	18.738	54	5069
	200	10 20	63.969	20	4704
	221	10 9	2.150	50	2573
	222	10 9	10.968	192	12516
	230	10 16	24.554	68	3330
	315	10 16	2.529	13	432
	328	10 9	1.351	86	140

Table 5 :

 5 Run-times of SMT method for some models under Ubuntu (U) compared to Windows (W). A star (*) signifies a program crash.

	BioModel Polyhedra	MathSAT (U) [s] Z3 (U) [s]	Yices (U) [s]	MathSAT (W) Z3 (W) [s]
	22	147	3.583	4.656	2.097	3.703	4.000
	32	244	34.598	42.938	35.319	37.234	37.156
	73	13449	2089.749	1314.049	*	2232.641	1288.156
	93	596	154.564	198.951	*	169.375	199.156
	183	1	90.184	84.033	131.141	94.922	102.938
	397	0	0.179	0.160	0.132	0.188	0.156
	430	1676	260.927	248.023	*	282.484	260.625
	501	916	216.917	297.374	*	237.719	292.672

Table 6 :

 6 Run-times of SMT method in seconds for some models without and with preprocessing using the MathSAT solver under Windows. Column "PP time" contains the time for preprocessing and column "after PP" the total time minus the preprocessing time.

	BioModel Dim Combos Polyhedra	Without PP	With PP	PP time	After PP	Speed-up
	21	10	10 8	46	1.3	1.5	0.6	1.0	0.87
	22	10	10 9	147	3.7	3.7	0.5	3.3	1.00
	26	11	10 6	6	0.2	0.4	0.3	0.1	0.50
	28	16	10 9	17	1.2	1.8	0.8	1.0	0.67
	30	18	10 10	6	0.4	1.3	1.0	0.3	0.31
	32	36	10 21	244	37.2	29.8	1.8	28.0	1.25
	41	10	10 14	4	0.3	2.0	1.9	0.1	0.15
	48	23	10 22	5	1.2	12.5	12.3	0.2	0.10
	61	22	10 25	10084	3191.3	2651.9	40.6	2611.3	1.20
	73	16	10 14	13449	2232.6	2203.2	1.4	2201.8	1.01
	74	19	10 15	9685	1873.7	2179.1	1.2	2177.8	0.86
	93	33	10 26	596	169.4	143.8	7.3	136.5	1.18
	102	13	10 10	322	24.6	20.9	2.4	18.5	1.18
	103	17	10 18	1938	382.8	343.1	38.2	304.8	1.12
	105	27	10 23	130	16.6	73.8	60.2	13.6	0.22
	147	24	10 16	54	11.7	12.3	6.0	6.3	0.95
	183	67	10 84	1	94.9	> 40000.0	> 40000.0	-	-
	200	22	10 20	20	3.9	10.1	8.0	2.2	0.39
	221	8	10 9	50	1.1	2.1	1.2	0.9	0.52
	222	8	10 9	192	6.1	6.9	1.5	5.4	0.88
	230	24	10 16	68	9.7	17.0	7.1	9.9	0.57
	315	19	10 16	13	1.6	4.0	2.9	1.1	0.40
	328	18	10 9	86	3.8	4.0	0.6	3.4	0.95
	365	30	10 23	70	14.2	18.1	10.6	7.5	0.78
	397	50	10 38	0	0.2	1.8	1.8	0.0	0.11
	407	47	10 16	212	37.1	30.2	2.0	28.1	1.23
	430	27	10 13	1676	282.5	255.8	1.6	254.2	1.10
	431	27	10 16	155	17.7	18.3	3.2	15.1	0.97
	477	43	10 17	467	59.0	55.0	1.7	53.3	1.07
	482	23	10 11	17	1.0	2.0	1.2	0.8	0.50
	489	35	10 21	42	7.5	22.3	14.8	7.5	0.34

Table 7 :

 7 Relative run-time of SMT method used for Searching another point, Minimizing a polyhedron and Inclusion checking. The second part of the models didn't finish and the numbers signify only the times until the process was terminated.

DFG SYMBIONT Project ANR-17-CE40-0036


Mariane El Kassis
email: mariane.el-kassis@mines-ales.fr

⇤ Franc ¸ois Trousset

Nicolas Daclin

Gregory Zacharewicz

Business Process Web-based Platform for Multi modeling and Simulation

Keywords: Model-Driven Engineering (MDE), Business Process Modeling and Notation (BPMN), discrete event simulation (DES), Discrete Event System Specification (DEVS), Model Transformation, Meta-Modeling, Business Process Simulation Interchange standard (BPSIM)

Enterprises and organizations need more streamlined business processes (BPs) to increase their performance. One of the main factors in building successful BPs is the stakeholders' engagement in BP management and improvement projects via modeling and simulation methodologies, which aspires to reduce conflicts, encourage innovation, and increase ownership and inclusive decision making. This paper attempts to highlight the existing development of concepts, methods, and tools bridging the gap between the business process modeling and simulation worlds. This paper will summarize current works and propose perspectives about a transformation approach of Business Process Modeling and Notation (BPMN) models that are extended with Business Process Simulation Interchange Standard (BPSIM) into Discrete Event System Specification (DEVS) simulation models based on a metamodel approach.

INTRODUCTION

Business process management is widely deployed and used in many organizations, and business processes are a significant part of the organization's running. It is used on three levels -first at a descriptive level to share and capitalize knowledge ; second at an executable level to orchestrate and automate business processes ; and third at an analysis level to improve performance. To this purpose, BPMN language is becoming an in-demand skill in the business industry [START_REF] Anwar-Ur-Rehman Pasha | bloom's taxonomy for standardizing bpm education of it undergraduates students[END_REF]. BPMN workflows are usually designed by business users or developers using standalone applications allowing stakeholders to translate their knowledge into conceptual models away from development, implementation, and environments. On the other hand, business process simulation is an essential element for the performance optimization of organizations [START_REF] Mourtzis | Simulation in the design and operation of manufacturing systems : state of the art and new trends[END_REF]. Although very important in business process management, simulation is still not systematically used in most business process improvement projects. One reason could be that simulation usually requires additional acquired techniques performed by technical users rather than the business stakeholders. Models can bring added-value support in the design of complex systems. Many large industries adopt the Model-Driven Engineering (MDE) approach in developing complex systems. It has proven its efficiency in various business domains such as aerospace, automobile, telecommunications, etc. [START_REF] Whittle | The State of Practice in Model-Driven Engineering[END_REF]. MDE concentrates on the source model conformity to a defined metamodel. The source model is then transformed into a target model following mapping rules and transformation techniques. In modeling extensive processes, business stakeholders collaborate to share their knowledge in defining the business process. Modeling and Simulation (M&S) methodologies recommend using conceptual models during the simulation model development process, which requires engaging business stakeholders in discrete event simulation improvement projects [START_REF] Fernandes | Methodology for the management of discrete event simulation projects based on pm-bok® : action research in a high-tech company[END_REF]. The purpose of this paper is to allow us to extend BPMN workflows using discrete-event simulation elements specified in BPSIM ; the extended model will be transformed into a DEVS model to operate the simulation. The paper is structured as follows. After a brief introduction in section 1, problems and expected outcomes are presented in section 2. Community contribution is detailed in sections 3,4,5. Our proposition is presented in section 6 (extension and transformation approach). Section 7 concludes the paper and opens outlines.

PROBLEMS AND EXPECTED OUTCOME

Enterprise business processes need to be formalized to be better understood and monitored. For this purpose, several standards have been developed. This paper will specifically study BPMN, a standard set by the Object Management Group (OMG) that provides a conventional notation allowing collaboration between business users, analysts, and implementers. BPMN is currently mature with over 20 years of use and development and widely expresses business processes. Several business process tools integrate graphical notation and software and are able to interpret process models in order to orchestrate them. However, these tools are limited in operational simulation capabilities [START_REF] Omg | Business Process Model and Notation[END_REF]. Selected free and open-source BPMN applications exist, which business stakeholders may use to design business processes. Such tools may allow model simulation but integrating the modeling part with the simulation part is still not comprehensive. For example, external simulation elements (scenario parameters, resources, duration, etc.) need to be added to the model to be simulated. Therefore, these two worlds remain separate, limiting the use of simulation to users accustomed to simulation principles. On the other hand, BPSIM is a standardized specification developed to simulate business models by enriching the conceptual model with scenario parameters that affect the process execution [START_REF] Wfmc | Business Process Simulation Specification[END_REF]. BPSIM aims to enrich selected BPMN elements with information on related parameters and logic significant for the simulation model. This paper discusses the transformation and simulation of a BPMN mo-del extended with BPSIM into a simulation model that conforms to the DEVS formalism to permit verification and simulation of the model properties. Consequently, we intend to develop a web-based platform to empower business stakeholders and engage them in discrete event simulations and business process improvement projects.

BPMN & ITS EXTENSION

BPMN is a business process modeling method used to describe an organization's value chains and business activities using a graphical representation (ISO/IEC 195103 [START_REF] Omg | Business Process Model and Notation[END_REF] ). The model extension is a technique that supports the enrichment of the model with external concepts originating from different domains. Extensibility is one of the main features of the BPMN metamodel as defined by OMG. BPMN extensions are utilized to visualize elements not available in BPMN like resources [START_REF] Jesús | Extending BPMN 2.0 : Method and tool support[END_REF] or key performance indicators [START_REF] Onggo | A BPMN extension to support discrete-event simulation for healthcare applications : an explicit representation of queues, attributes and data-driven decision points[END_REF]. On the other hand, resource allocation mechanisms and failure of resources are represented as BPMN extensions in the context of discrete event simulation. BPMN extensions are divided into conceptual enrichment and operational enrichment. The conceptual enrichment includes the extensions of resources by adding attributes like resource description, nature, and capacity. The operational enrichment consists of describing the behavior of the enriched element. For example, tasks may be extended in order to describe the interoperability behavior [START_REF] Mallek | The application of interoperability requirement specification and verification to collaborative processes in industry[END_REF]. BPMN users may enrich a specific metamodel by adding extension elements and attaching additional attributes to existing BPMN elements, and the model would remain BPMN-Compliant [START_REF] Omg | Business Process Model and Notation[END_REF].

BPMN TRANSFORMATION & SIMULATION

MDE is a software engineering methodology that promotes the development of models at different levels of abstraction, raising the level of development to models rather than pro-grams. The core principles of MDE are the use of models and model transformation. The principle of the models concentrates on the notion of meta-model and conformance, while the principle of model transformation includes the concept of mappings and transformation rules [START_REF] Mens | A Taxonomy of Model Transformation[END_REF]. The transformation of high-level workflows into simulation models was a main area of study by many researchers in the last two decades.

In the context of MDE and metamodel based transformation, both source and target models must conform to their respective metamodel.

While BPMN conforms to the BPMN metamodel standardized by OMG [START_REF] Omg | Business Process Model and Notation[END_REF], simulation models like DEVS do not have a specific standardized metamodel yet researchers worked to define a model-driven development framework for M&S. Such framework presented BPMN and DEVS metamodels and a defined set of rules to transform a conceptual source model BPMN into a simulation model DEVS [START_REF] Cetinkaya | MDD4MS : A model driven development framework for modeling and simulation[END_REF]. BPMN elements differ in their internal behavior. These elements may be considered black boxes compared to white boxes DEVS models. To simulate a BPMN flow in DEVS, elements of this flow need to be mapped to their corresponding DEVS representation. Over the past years, researchers worked on providing a mapping for different BPMN elements. However, only a few elements were studied and taken into consideration. Since no standard DEVS metamodel is currently adopted, every group of researchers developed their own DEVS metamodel [START_REF] Cetinkaya | MDD4MS : A model driven development framework for modeling and simulation[END_REF] [7] [START_REF] Bazoun | Business Process Simulation : Transformation of BPMN 2.0 to DEVS Models[END_REF]. Authors in [START_REF] Cetinkaya | MDD4MS : A model driven development framework for modeling and simulation[END_REF] introduced the first mapping concepts which were then extended covering additional areas which includes additional types of BPMN tasks [START_REF] Bazoun | Business Process Simulation : Transformation of BPMN 2.0 to DEVS Models[END_REF] and the notion of resource allocation and failure [START_REF] Ambrogio | Resource-based modeling and simulation of business processes[END_REF]. Nevertheless, only a few BPMN concepts are integrated into these works. They still miss some essential features such as intermediate events, interrupting events, message flows, etc. In addition to the transformation of BPMN to DEVS, many researchers analyzed the transformation to other discrete event simulation models such as Petri nets [START_REF] Mutarraf | Transformation of Business Process Model and Notation models onto Petri nets and their analysis[END_REF]. However, it is more convenient to simulate models using DEVS rather than Petri nets as the notion of time is not natively defined in Petri nets [START_REF] Bernard | Theory of Modeling and Simulation[END_REF]. Moreover, DEVS provides a generalized framework of the modeling and simulation systems and offers a formal definition of the simulator [START_REF] Zacharewicz | Flattening G-DEVS / HLA structure for Distributed Simulation of Workflows[END_REF].

SIMULATION APPLICATIONS

To have simulation results as realistic as possible, it is necessary to add simulation parameters that best describe the real world. The BP simulation tools currently available in the market differ in terms of simulation capabilities which affects the accuracy of the results. Simulation tools may be segmented into three categories [START_REF] Pufahl | Design of an Extensible BPMN Process Simulator[END_REF] :

-Business process modeling or management tools with simulation support (e.g., ADONIS, ARIS Toolset, Bizagi Modeler, L-SIM, Simul8, PragmaDEV) -General-purpose simulation tools (e.g., Arena, AnyLogic) -Stand-alone business process simulators (e.g., Bimp) Academia and the industry developed several DEVS tools. These tools differ by design goals and specific programming language implementations. An evaluation of available DEVS simulators was performed comparing them in terms of formalism, compliance, features, and performance [START_REF] Van Tendeloo | An evaluation of DEVS simulation tools[END_REF].

BPMN is designed to represent business processes without any concerns for simulation. Our goal is to be able to simulate BPMN models by adding the necessary simulation parameters. BPSIM is a specification by the Workflow Management Coalition (WFMC) that defines the parameterization of business process models from different perspectives allowing process analysis, simulation, and optimization [START_REF] Wfmc | Business Process Simulation Specification[END_REF]. We have then chosen to extend BPMN2.0 standard with BPSIM2.0 standard for the specification level. The extended BPMN (BPMN+BPSIM) model may then be simulated with several existing BPSIM-able Simulation engines, which understands and process the semantics existing in the extended model. At the operational level, we have chosen the transformation into the DEVS model as it already considers the temporal dimension and because a lot of work has already been conducted in this area. Moreover, the existing BPSIM simulators are proprietary and thus difficult to adapt to specific contexts. In this context, we aim to transform BPMN models extended with BPSIM into DEVS models that can be simulated/operated to detect errors and improve the business process as described in Figure 1. The extended BPMN model (BPMN+BPSIM) is transformed to its respective DEVS model based on a defined set of mapping and transformation rules. The final goal is to obtain a simulation tool that is open, extensible, and based on standards.

Extension with BPSIM

BPSIM 2.0 open-source specification document offers a straightforward method to incorporate BPSIM attributes and semantics onto BPMN. The purpose of this work is to enrich the conceptual model by extending selected BPMN elements with information on related parameters and logic that are significant for the simulation model issued from BPSIM. The BP-SIM perspective consists of five parameters described in the BPSIM specification [START_REF] Wfmc | Business Process Simulation Specification[END_REF] (time, control, cost, property and priority).

Every parameter of this perspective is a class group containing sub-elements defined based on their role in the simulation. BPSIM standard checks the validity of assigning the parameter to a specific BPMN element using a set of rules listed in the BPSIM specification document. Not all parameters are valid to be assigned to different BPMN elements. For instance, time parameters may be assigned to BPMN tasks and activities which have their details defined, unlike subprocesses. Control Parameters may not be assigned to connecting objects, swim lanes, artifacts, and attributes. Functionality is also represented at the level of the BPMN element type, such that within property parameters, message flows may be assigned property but not queue length. BPSIM provides several additional features to BPMN, such as : probabilistic distribution, assigning resources to activities, probabilistic duration, probability of Branch selection, duration parameter, priority of resources, etc.

Transformation from (BPMN + BP-SIM) to DEVS

As our goal is to transform BPMN+BPSIM models into DEVS models, we need a metamodel for both source and target and a language to express rules of transformation. The BPMN metamodel is defined clearly as UML in the BPMN 2.0 OMG specification [START_REF] Omg | Business Process Model and Notation[END_REF], while BPSIM elements are also described in the BPSIM specifications [START_REF] Wfmc | Business Process Simulation Specification[END_REF]. Nevertheless, we still need to acquire a unique BPMN+BPSIM metamodel that we intend to build from existing BPMN and BPSIM ones. For the target model, currently, there is no unique metamodel for DEVS. The existing DEVS metamodels are dependent on the environment where the model is being simulated. The need for a universal metamodel for DEVS has been addressed in the research community who proposed partial metamodels of DEVS. Although these propositions have been used in conceptual and technical works, no consensus and standardization were reached.

In the most recent work, a conceptualization of DEVS formalism using a universal representation of DEVS metamodel was suggested to combine both theory and practice points of view. It incorporates all DEVS views (abstraction, formalization, and implementation) into a traceable network of metamodels to be used in a compliant fashion [START_REF] Maria | Towards a Universal Representation of DEVS : A Metamodel-Based Definition of DEVS Formal Specification[END_REF]. To set rules for the transformation, we need a specialized language : several ones have been proposed to achieve model transformations, out of which are the extensible stylesheet language transformation (XSLT) and Atlas Transformation Language (ATL). XSLT is a standard used to process XML documents by transforming the XML to another schema or format after identifying mapping patterns while ATL is a model transformation language and toolkit developed on top of Eclipse. ATL provides ways to produce target models from a set of source models. ATL tools support the major tasks involved in using a language : editing, compiling, executing, and debugging [6]. ATL method was chosen to transform BPMN to DEVS in [START_REF] Bazoun | Business Process Simulation : Transformation of BPMN 2.0 to DEVS Models[END_REF].

FIGURE 1 -

 1 FIGURE 1 -Extension and Transformation architecture