Georges Habib
email: georges.habib@mis.mpg.de

Ken Richardson
email: k.richardson@tcu.edu

A brief note on the spectrum of the basic Dirac operator

In this paper, we prove the invariance of the spectrum of the basic Dirac operator defined on a Riemannian foliation (M, F) with respect to a change of bundle-like metric. We then establish new estimates for its eigenvalues on spin flows in terms of the O'Neill tensor and the first eigenvalue of the Dirac operator on M . We discuss examples and also define a new version of the basic Laplacian whose spectrum does not depend on the choice of bundle-like metric.

Introduction and Results

Let (M, F ) be a smooth, closed manifold of dimension n + 1 endowed with a foliation F given by an integrable subbundle L ⊂ T M of rank p. We assume throughout the paper that the foliation is Riemannian and is endowed with a bundle-like metric [START_REF] Reinhart | Foliated manifolds with bundle-like metrics[END_REF], i.e. the metric g on M induces a holonomy invariant transverse metric on the normal bundle Q = T M/L ≃ L ⊥ , of rank q = (n+1)-p; this means that the transverse Lie derivative L X g| Q is zero for all leafwise vector fields X ∈ Γ(L). This condition is characterized by the existence of a unique metric and torsion-free connection ∇ on Q [START_REF] Molino | Riemannian foliations[END_REF], [START_REF] Reinhart | Foliated manifolds with bundle-like metrics[END_REF], [START_REF] Ph | Geometry of foliations[END_REF]. We can then associate to ∇ the transversal curvature data, in particular the transversal Ricci curvature Ric ∇ and transversal scalar curvature Scal ∇ .

Many researchers have studied basic forms and the basic Laplacian on Riemannian foliations [START_REF] Álvarez-López | The basic component of the mean curvature of Riemannian foliations[END_REF], [START_REF] Kamber | De Rham-Hodge theory for Riemannian foliations[END_REF], [START_REF] Ph | Geometry of foliations[END_REF]. Basic forms are locally forms on the space of leaves; that is, forms α satisfying X α = X dα = 0 for all X ∈ Γ(L). They are preserved by the exterior derivative and are used to define basic cohomology groups, which are always finite-dimensional in the case of Riemannian foliations. The basic Laplacian for a given bundle-like metric is a version of the Laplace operator that preserves the basic forms and that is essentially self-adjoint on the L 2 -closure of the space of basic forms. Its spectrum depends on the choice 1 of the bundle-like metric and provides invariants of that metric. See [START_REF] Jung | Transverse conformal Killing forms and a Gallot-Meyer Theorem for foliations[END_REF], [START_REF] Lee | Lichnerowicz and Obata theorems for forliations[END_REF], [START_REF] Lee | Riemannian foliations and eigenvalue comparison[END_REF], [START_REF] Park | The basic Laplacian of a Riemannian foliation[END_REF], [START_REF] Richardson | Asymptotics of heat kernels on Riemannian foliations[END_REF], [START_REF] Richardson | Traces of heat operators on Riemannian foliations[END_REF] for further results.

We now discuss the construction of the basic Dirac operator (see [START_REF] Douglas | Index formulas for geometric Dirac operators in Riemannian foliations[END_REF], [START_REF] Glazebrook | Transversal Dirac families in Riemannian foliations[END_REF], [START_REF] Prokhorenkov | Natural Equivariant Dirac Operators[END_REF]). Let (M, F ) be a Riemannian manifold endowed with a Riemannian foliation. Let E → M be a foliated vector bundle (see [START_REF] Kamber | Foliated bundles and characteristic classes[END_REF]) that is a bundle of Cl(Q) Clifford modules with compatible connection ∇ E . The transversal Dirac operator D tr is the composition of the maps

Γ (E) (∇ E ) tr -→ Γ (Q * ⊗ E) ∼ = -→ Γ (Q ⊗ E) Cliff -→ Γ (E) ,
where the last map denotes Clifford multiplication, denoted by " • ", and the operator ∇ E tr is the projection of ∇ E . The transversal Dirac operator fixes the basic sections Γ b (E) ⊂ Γ(E) (i.e. Γ b (E) = {s ∈ Γ(E) : ∇ E X s = 0 for all X ∈ Γ(L)}) but is not symmetric on this subspace. By modifying D tr by a bundle map, we obtain a symmetric and essentially self-adjoint operator D b on Γ b (E). To define D b , first let H = p i=1 π(∇ M fi f i ) be the mean curvature of the foliation, where π : T M → Q denotes the projection and {f i } i=1,••• ,p is a local orthonormal frame of L. Let κ = H ♭ be the corresponding 1-form, so that H = κ ♯ . Let κ b := P b κ be the L 2 -orthogonal projection of κ onto the space of basic forms (see [START_REF] Álvarez-López | The basic component of the mean curvature of Riemannian foliations[END_REF], [START_REF] Park | The basic Laplacian of a Riemannian foliation[END_REF]), and let κ ♯ b be the corresponding vector field. We now define

D b s := 1 2 (D tr + D * tr )s = q i=1 e i • ∇ E ei s - 1 2 κ ♯ b • s,
where

{e i } i=1,••• ,q is a local orthonormal frame of Q. A direct computation
shows that D b preserves the basic sections, is transversally elliptic, and thus has discrete spectrum [START_REF] Kacimi-Alaoui | Opérateurs transversalement elliptiques sur un feuilletage riemannien et applications[END_REF].

In this paper, we study the invariance of the spectrum of the basic Dirac operator with respect to a change of bundle-like metric; that means when one modifies the metric on M in any way that leaves the transverse metric on the normal bundle intact (this includes modifying the subbundle L ⊥ ⊂ T M , as one must do in order to make the mean curvature basic, for example). We prove Theorem 1.1 Let (M, F ) be a compact Riemannian manifold endowed with a Riemannian foliation and basic Clifford bundle E → M . The spectrum of the basic Dirac operator is the same for every possible choice of bundle-like metric that is associated to the transverse metric on the quotient bundle Q.

We emphasize that the basic Dirac operator D b depends on the choice of bundle-like metric, not merely on the Clifford structure and Riemannian foliation structure, since both projections T * M → Q * and P b as well as κ b depend on the leafwise metric. We clarify this result on the 2-dimensional torus (see Example 2 at the end of Section 3). In Section 4, we show that the square of the standard basic Dirac operator on basic forms is related to but not the same as the basic Laplacian, whose spectrum does depend on the choice of compatible bundle-like metric.

In [START_REF] Habib | Tenseur d'impulsion-énergie et feuilletages[END_REF], the first author studied the spectrum of the basic Dirac operator on spin foliations. That is, a compact Riemannian manifold (M, F ) endowed with a Riemannian foliation F of codimension q such that the normal bundle carries a spin structure (as a foliated vector bundle). In that paper, the author obtained a Friedrich-type estimate [START_REF] Th | Der erste Eigenwert des Dirac-Operators einer kompakten Riemannshen Mannigfaltigkeit nichnegativer Skalar-krümmung[END_REF] for the eigenvalues of the basic Dirac operator for bundle-like metrics with basic-harmonic mean curvature. Since by Theorem 1.1 the spectrum does not change, we deduce that for any bundle-like metric the eigenvalues of the basic Dirac operator satisfy

λ 2 ≥ q 4(q -1) inf M (Scal ∇ ). (1.1)
This estimate is of interest only for positive transversal scalar curvature. Moreover, we show Proposition 1.2 Let (M, F ) be a compact Riemannian manifold endowed with a spin foliation with basic mean curvature κ. Then, we have the estimate

λ 2 ≥ q 4(q -1) inf M (Scal M -Scal L + |A| 2 Q + |T | 2 L ).
If F is a Riemannian flow (i.e. p = 1), then

λ 2 ≥ q 4(q -1) inf M (Scal M + |A| 2 Q + |κ| 2 ). (1.2)
If the limiting case is attained, the foliation is minimal and we have a transversal Killing spinor.

Here A and T denote the O'Neill tensors [START_REF] Besse | Einstein manifolds, Ergeb. Math. Grenzgeb[END_REF][START_REF] O'neill | The fundamental equations of a submersion[END_REF] of the foliation (see Section 3 for the definitions). We point out that Estimate (1.2) may improve Inequality (1.1) (see Example 1). On the other hand, using the Min-Max principle we give a new lower bound for the eigenvalues of the basic Dirac operator on a Riemannian flow (M, F ) in terms of the first eigenvalue of the Dirac operator of M and the O'Neill tensor of the flow. We prove Proposition 1.3 Let (M n+1 , F ) be a compact spin manifold endowed with a Riemannian flow F given by a unit vector field ξ. Then we have the estimate

λ 2 ≥ 1 2 λ 2 (D M ) - n 16 sup M (|A| 2 Q ), (1.3)
where λ(D M ) is the first eigenvalue of the Dirac operator of M .

2 Proof of Theorem 1.1

Suppose that D b is any basic Dirac operator associated to a basic Clifford bundle over a Riemannian foliation with bundle-like metric. Consider a change of the bundle-like metric as mentioned in the introduction. When one does this, there is no effect on the compatibility of the new transversal connection, since the bundles and bundle metrics are unchanged. Furthermore, the set of all basic sections of the Clifford bundle is unchanged. Let D ′ b denotes the new basic Dirac operator, dvol ′ the new volume form, and κ ′ b the basic component of mean curvature with respect to the new bundle-like metric. We see that dvol ′ = h dvol, for some positive smooth function h. Then the new L 2 -metric on basic sections s, t is

< s, t > ′ = M (s, t)h dvol,
where ( , ) is the pointwise inner product on Γ(E). As in [START_REF] Park | The basic Laplacian of a Riemannian foliation[END_REF], we have

< s, t > ′ = M (s, t)P b h dvol. Setting α = P b h, we see that < α -1/2 D b (α 1/2 s), t > ′ = M (α -1/2 D b (α 1/2 s), t)α dvol = M (D b (α 1/2 s), α 1/2 t) dvol = M (α 1/2 s, D b (α 1/2 t)) dvol,
where the last equality is a consequence of the fact that D b is self-adjoint with respect to the old L 2 -metric. Hence,

< α -1/2 D b (α 1/2 s), t > ′ = M (s, α -1/2 D b (α 1/2 t))α dvol = < s, α -1/2 D b (α 1/2 t) > ′ , so that α -1/2 D b α 1/2 = D b + 1 2α (dα)
• is self-adjoint with respect to the new L 2 -metric. Now, we compute

D ′ b -α -1/2 D b α 1/2 = D ′ b -D b - 1 2α (dα) • = D tr - 1 2 κ ′♯ b • -D tr + 1 2 κ ♯ b • - 1 2α (dα) • = - 1 2 κ ′♯ b • + 1 2 κ ♯ b • - 1 2α (dα) • .
By the symmetry of

D ′ b -α -1/2 D b α 1/2 with respect to the new L 2 -metric, the operator κ ′♯ b • -κ ♯ b • -1 2α (dα)
• must be a Hermitian symmetric endomorphism of Γ(E). Since Clifford multiplication by cotangent vectors is skew-Hermitian, this bundle endomorphism must be zero. Since D ′ b is just a conjugate of D b , all the spectrum is unchanged.

Remark. Dom ínguez showed in [START_REF] Dom Ínguez | A tenseness theorem for Riemannian foliations[END_REF] that every Riemannian foliation admits a bundle-like metric with basic mean curvature (so that κ = κ b ). Further, Mason showed in [START_REF] Mason | An application of stochastic flows to Riemannian foliations[END_REF] that the bundle-like metric may be chosen so that the mean curvature is basic-harmonic (∆ b κ = 0). Thus, we may assume that the mean curvature form has this property when calculating the eigenvalues of the basic Dirac operator.

3 Proof of Propositions 1.2 and 1.3

We recall that the geometry of the normal bundle of a Riemannian foliation is characterized by the existence of the O'Neill tensors, given for all X, Y ∈ Γ(T M ) by

A X Y = π ⊥ (∇ M π(X) π(Y )) + π(∇ M π(X) π ⊥ (Y )), and
T X Y = π ⊥ (∇ M π ⊥ (X) π(Y )) + π(∇ M π ⊥ (X) π ⊥ (Y ))
, where π ⊥ : T M → L denotes the projection. Indeed, for all Z ∈ Γ(Q) the tensor A Z is a skew-symmetric map interchanging L and Q and we have that

A Z W = 1 2 π ⊥ [Z, W ] for all W ∈ Γ(Q).
Moreover, the scalar curvatures of M , Q and L are related by the following relation [2, p. 244]

Scal M = Scal ∇ + Scal L -|A| 2 Q -|T | 2 L -|κ| 2 + 2div Q κ (3.4)
where

| • | Q (resp. | • | L ) is the norm evaluated on a local frame of Γ(Q) (resp. Γ(L)).
For the proof of Proposition 1.2, we assume that the mean curvature is basic and use the transversal Schrödinger-Lichnerowicz formula for D 2 b [7, Eq. 1.6] to write for any basic spinor ψ

D 2 b ψ = ∇ * ∇ψ + 1 4 (Scal ∇ + |κ| 2 )ψ - 1 2 δ b κ • ψ, (3.5)
where δ b is the adjoint of d on the space of basic forms [START_REF] Park | The basic Laplacian of a Riemannian foliation[END_REF]. Plugging Equation (3.4) in (3.5) finishes the proof of the proposition after integrating over M the scalar product of (3.5) by ψ and using the fact that div

Q κ = -δ b κ + |κ| 2 if κ is basic and that |∇ψ| 2 ≥ 1 q λ 2 |ψ| 2 .
Remark. In [START_REF] Habib | Tenseur d'impulsion-énergie et feuilletages[END_REF], we consider a Riemannian flow F and modify the metric of M in the direction of the flow, i.e. we consider the family of metrics g t = t 2 g ξ ⊕g Q . The spectrum of the Dirac operator of (M, g t ) tends to the spectrum of the basic Dirac operator D b when t tends to zero. Hence by Friedrich's estimate, we have for the Dirac operator of (M, g t ) that λ 2 t ≥ q+1 4q inf M Scal t M . A direct computation for the scalar curvature leads to Scal t M = Scal M + (1 -t 2 )|A| 2 Q . Thus when t tends to zero, we obtain the estimate for the eigenvalues of the basic Dirac operator λ 2 ≥ q+1 4q inf M (Scal M + |A| 2 Q ) which is not as strong as (1.2) since q (q-1) ≥ q+1 q .

The proof of Proposition 1.3 is based on the use of the Min-Max principle by computing the Rayleigh quotient ||DM ψ|| 2 ||ψ|| 2 associated with a basic spinor ψ, where || • || is the L 2 -norm on the spinor bundle ΣM induced by the Hermitian metric. For this, we recall that the spinorial Levi-Civita connections ∇ M and ∇ are related for any spinor field ψ ∈ Γ(ΣM ) by [START_REF] Habib | Tenseur d'impulsion-énergie et feuilletages[END_REF]Eq. 2.4.7]   

∇ M ξ ψ = ∇ ξ ψ + 1 4 n i=1 e i • A ei ξ • ψ + 1 2 ξ • κ • ψ ∇ M Z ψ = ∇ Z ψ + 1 2 ξ • A Z ξ • ψ, where Z ∈ Γ(Q) and {e i } i=1,••• ,n is an orthonormal frame of Γ(Q).
We easily deduce the relation between the basic Dirac operator and the Dirac operator of M :

D M ψ = D b ψ - 1 4 n i=1 ξ • e i • A ei ξ • ψ,
for each basic spinor field ψ. Using the pointwise inequality

| n i=1 ξ • e i • A ei ξ • ψ| ≤ n i=1
|A ei ξ||ψ| and the fact that |a + b| 2 ≤ 2a 2 + 2b 2 , we deduce that for any basic eigenspinor ψ with eigenvalue λ, we have

|D M ψ| 2 ≤ 2λ 2 |ψ| 2 + 2( 1 4 ) 2 ( n i=1
|A ei ξ|) 2 |ψ| 2 ≤ 2λ 2 |ψ| 2 + n 8 |A| 2 Q |ψ| 2 ,
where the last inequality is a consequence of the Cauchy-Schwarz inequality. Thus we get that ||DM ψ|| 2

||ψ|| 2 ≤ 2λ 2 + n 8 sup M (|A| 2 Q )
, and by the Min-Max principle the proof of the proposition follows. Now we will consider particular Riemannian flows on the 3-dimensional sphere and the 2-dimensional torus and will estimate the eigenvalues of the basic Dirac operator.

Example 1. Let S 3 ⊂ C 2 be the round 3-sphere equipped with its standard metric of sectional curvature equal to 1 and its standard spin structure. For r > 0, we define the 1-parameter group γ r t of isometries by [START_REF] Hebda | An example relevant to curvature pinching theorems for Riemannian foliations[END_REF] γ r t (z, w) = (e irt z, e it w) for t ∈ R and (z, w) ∈ S 3 . This group generates the Killing vector field Z r ∈ T S 3 defined by Z r (z,w) = (irz, iw) and thus the vector field ξ r (z,w) = Z r (z,w) /|Z r (z,w) | defines a Riemannian flow F r on S 3 . Let us consider on T (z,w) S 3 the two vector fields X r On the one hand, Inequality (1.1) allows to deduce an estimate for the eigenvalues of the basic Dirac operator

λ 2 > 1 2 inf S 3 (2 + 6 r 2 r 2 |z| 2 + |w| 2 ) =    1 + 3r 2 if r < 1 4 if r > 1.
On the other hand, by Inequality (1.2) we get that

λ 2 > 1 2 inf S 3 (6 + 2r 2 + (1 -r 2 ) 2 |z| 2 |w| 2 (r 2 |z| 2 + |w| 2 ) 2 ) =    r 2 + 3 if r < 1 1 r 2 + 3 if r > 1.
We see that for r < 1 the last estimate is better than the first one. The lower bound in Inequality (1.3) is weaker than those above since we get that

λ 2 > 9 8 -1 4r 2 if r < 1 9 8 -r 2 4 if r > 1. Example 2. Consider the manifold (S 1 × S 1 , f 2 dθ 2 ⊕ dt 2 ) where f (θ, t) is a function on S 1 × S 1 = R 2 /(2πZ) 2 . The frame ξ = 1 f ∂ θ , e 1 = ∂ t is
∂ ∂t ψ + i ġ(t) 2g(t) ψ = λψ.
The solution of this equation is ψ(t) = cg -1 2 (t)e -iλt , where c is a complex constant. For the trivial spin structure, the spinor ψ is 2π-periodic if λ = k ∈ Z. This verifies that the spectrum does not depend on the choice of the function f , as Theorem 1.1 shows. This example may be easily modified to provide a similar example on any suspension of an isometry ϕ : M → M , yielding a one-dimensional foliation of an M -bundle over a circle. That is, the manifold is R × ϕ M , the quotient of R × M by the equivalence relation ( θ, x ) ∼ ( θ -2π, ϕ(x) ).

Remarks about the basic Laplacian

As mentioned in the introduction, it is well-known that the eigenvalues of the basic Laplacian depend on the choice of bundle-like metric; for example, in [START_REF] Richardson | Traces of heat operators on Riemannian foliations[END_REF]Corollary 3.8], it is shown that the spectrum of the basic Laplacian on functions determines the L 2 -norm of the mean curvature on a transversally oriented foliation of codimension one. If the foliation were taut, then a bundlelike metric could be chosen so that the mean curvature is identically zero, and other metrics could be chosen where the mean curvature is nonzero. This is one reason why the invariance of the spectrum of the basic Dirac operator is such a surprise.

If one considers the bundle * Q * as a Cl(Q) bundle with Clifford multiplication α ♯ • = α ∧ -α for (local) basic one-forms α, then the corresponding basic Dirac operator is

D b = d + δ b -κ b - 1 2 κ ♯ b • .
The square of this operator is certainly not the basic Laplacian ∆ b , but ∆ b and D 2 b do have the same principal symbol. Observe that since dκ b = 0 (see [START_REF] Álvarez-López | The basic component of the mean curvature of Riemannian foliations[END_REF], [START_REF] Park | The basic Laplacian of a Riemannian foliation[END_REF]), the following operator is a differential on the space of basic forms:

d := d - 1 2 κ b ∧,
and its L 2 -adjoint restricted to the space of basic forms is

δ b = δ b - 1 2 κ b . Then d + δ b = d + δ b - 1 2 κ b - 1 2 κ b ∧ = D b .
Thus, with this new twisted differential on basic forms, the analogy with ordinary manifolds carries over, and the resulting basic D 2 b -harmonic forms would play the same role in Hodge theory with the twisted basic d-cohomology. Further, unlike the ordinary and well-studied basic Laplacian, the eigenvalues of D 2 b are invariants of the Riemannian foliation structure alone and independent of the choice of compatible bundle-like metric.

2 r 2 r 1 e r 2 = r r 2

 222 (z,w) = (|w| 2 z, -|z| 2 w) and Y r (z,w) = (i|w| 2 z, -ir|z| 2 w) for (z, w) ∈ S 3 . The frame {ξ r = Z r /|Z r |, e r 1 = X r /|X r |, e r 2 = Y r /|Y r|}, is orthonormal of the tangent space of S 3 . Hence by a direct computation with respect to this frame, on gets that the transversal scalar curvature is equal to Scal ∇ (z,w) = 2 + 6r |z| 2 +|w| 2 and the mean curvature κ r of the flow is a basic 1-form and equal to(1-r 2 )|z||w| r 2 |z| 2 +|w| 2 e r 1 . The O'Neill tensor of the flow is A e |z| 2 +|w| 2 ξ r .

  orthonormal with respect to the metric and we have [ξ, e 1 ] = 1 f ∂f ∂t ξ. The vector field ξ defines a Riemannian flow on M with vanishing O'Neill tensor and mean curvature is κ ♯ = -1 f ∂f ∂t e 1 . The basic component of the mean curvature satisfies κ ♯ b = -ġ(t) g(t) e 1 , where g(t) = 1 2π 2π 0 f (θ, t) dθ. In order to compute the spectrum of the basic Dirac operator, we will find solutions of the differential equation D b ψ = λψ, with λ is an eigenvalue of D b . Clifford multiplication by e 1 is multiplication by i. Hence the equation D b ψ = λψ can be written as ∂ ∂θ ψ = 0 and i


Rony Ghostine
email: rony.ghostine@ensem.inpl-nancy.fr

Jean-Marc Thiriet
email: jean-marc.thiriet@ujf-grenoble.fr

Jean-François Aubry
email: jean-francois.aubry@ensem.inpl-nancy.fr

Michel Robert
email: michel.robert@esstin.uhp-nancy.fr

A Framework for the Reliability Evaluation of Networked Control Systems

Keywords: Networked control systems, Reliability analysis, Petri-net, transient faults

this paper presents a framework to enable the analysis of the influence of the transmission faults on the reliability of a networked control system (NCS). The approach is composed of two parts: a modelling part in which all the basic components of a networked control system are modelled and a simulation part in which simulation is done on the models to evaluate the reliability. Due to external perturbations transmission faults may occur on the medium decreasing network quality of service and system performance. These aspects are difficult to assess with traditional dependability method like fault trees and reliability blocks. Our approach is applied to a case study example. The results show that our framework is an effective way for the reliability evaluation of networked control systems (

INTRODUCTION

Control systems with spatially distributed components have existed for several decades. Examples include chemical processes, power plants, airplanes, etc…. in such systems the components were connected via point-to-point connections and the systems were designed to bring all information from the sensors to a central location to take a decision on how to act (Halevi and Ray, 1988). Physical setups and expanding functionality are pushing the limits of the point-to-point architecture. Hence, such centralized point-to-point control systems are no longer suitable to meet new requirements such as modularity, decentralization of control, integrated diagnostics, quick and easy maintenance, and low cost. Technology advances and the availability of network connectivity have prompted the idea of introducing network facilities to control systems. Such systems are called networked control systems (NCSs): Their sensors, actuators, estimator units, and control units are connected through communication networks. This type of system provides several advantages such as modular and flexible system design, simple and fast implementation, and powerful system diagnosis and maintenance utilities. The disadvantage is that the analysis and design of an NCS becomes complex. Conventional control theory with many ideal assumptions, such as synchronized control and non-delayed sensing and actuation, has to be re-evaluated before it can be applied to NCSs. Specially, the problems of time-varying transmission periods, network schedulability, network induced delay, and packet loss, are of concern [START_REF] Soglo | Networked Control System Simulation Design and Its Application[END_REF]. [START_REF] Branicky | Networked controlsystem co-simulation for co-design[END_REF] developed a simulation tool that combines dynamic-system simulation for the control agents and environment with packet-level network simulation for the communications by extending network simulator-2 (NS2).

Due to the complexity of networked control system architectures, it is not trivial to evaluate their dependability level. The present paper aims at bringing a contribution relative to this aspect.

RELIABILIY OF CONTROL SYSTEM

Fault-trees and reliability blocks diagrams are the easiest and most often used techniques in complex systems dependability assessment. Many people have refined these techniques which have been applied to various industries, including aerospace, medical, and nuclear… These techniques are also called Boolean models. Their aim is to show how a binary system (with two states) state depends on the binary states of the system's components. These methods are not at all suited to modelling systems in which there are strong dependencies between components.

The assumption of components independence is precisely what makes Boolean models so powerful, but this assumption is extremely restrictive, and may prove to be totally unrealistic and lead to grossly erroneous results for some kinds of systems. These methods can not be directly applied to control systems study for two main reasons: -These methods do not take into account the transmission delay between nodes. This is a problem since control systems are known to be real time systems and delays are known to be able to degrade or destabilize the system [START_REF] Zhang | Stability of networked control systems[END_REF]). -The failure definition for control systems depends on the difference between the desired value of a physical property and its actual value, hence the need to calculate this actual value at each step and to compare it with the reference input. One of the major problems inherent to any reliability study of an industrial system is to take into account, in an effective and realistic way, the dynamic interactions existing between the physical parameters (pressure, temperature, flow rate, level,...) and the nomina1 or dysfunctional behaviour of the components of the system itself. In [START_REF] Dutuit | Dependability modelling and evaluation by using stochastic Petri nets: application to two test cases[END_REF] the authors present a Petri nets approach for a dependability evaluation of a control process, only permanent faults are considered.

In [START_REF] Askerdal | A Control Theory Approach for analysing the effects of data errors in safety-critical control system[END_REF] the authors developed a control theory methodology for analyzing the effects that data errors may have on the control system dependability. The effect is measured as the resulting control error (defined as the difference between the desired value of a physical property and its actual value).

In [START_REF] Moncelet | Dependability evaluation of a simple mechatronic system using coloured Petri nets In: Workshop on Practical Use of Coloured Petri Nets and Design CPN[END_REF] authors proposed an approach to find feared scenarios in control system. Like the works listed below only permanent faults are considered.

Unlike regular control systems, in networked control systems the synchronization between different sensors, actuators and control units is not guaranteed. Furthermore, there is no guarantee for zero delay or even constant delay in sending information from sensors to the control units and control signals from the control units to the actuators; moreover when there is congestion in the communication network, some packets are dropped. In real time systems, particularly control systems, delays or dropped packets may be catastrophic and may cause instability in the control system.

In [START_REF] Barger | Safety analysis and reliability estimation of a networkedcontrol system In: SAFEPROCESS 2003[END_REF], the authors take into consideration the message loss by assigning a fixed probability.

In [START_REF] Jumel | Safety Evaluation of controlled system distributed on TTA Architecture[END_REF], the authors present a methodology to study the safety of a mechatronic function distributed on a network. They have illustrated this methodology on a brake function mapped on the Time Triggered Architecture (TTA) in presence of transient faults. Thanks to Markov's Chains, they compute some important properties of safety. They take into account different types of faults (byzantine, and non_byzantine) they consider only one type of degradation and they assume that the vehicle's speed is decreased from Mmax in absence of error and Mmin in presence of errors. This assumption tends to simplify the problem and doesn't reflect the real behaviour of a dynamic system, where the faults effect depends on the system dynamics.

Industrial environment is specified by the existence of electromagnetic interferences (EMI). These interferences generate faults in electronic circuits that affect the normal operations. These faults are transient faults which means that the component affected is temporary unavailable, for example if a communication network is affected by a transient fault it may be unable to transmit data for a certain interval of time. These faults lunch the error detection, fault location and recovery mechanism, during this time the component is unavailable for the system mission, this time degrades the system performance and may even lead to a dynamic failure if the delay exceeds a certain limit. In communication systems, transient faults usually affect the medium leading to transmission error. In this study we try to propose an approach to include the influence of transmission error in the reliability analysis of networked control system. Note that external interference occur stochastically in time, this will lead to variable delays on affected messages.

Our work can be seen as an extension of Barger et al work by including variable delays in the dependability study, and Jumel et al work by taking into account the real behaviour of dynamic systems. While this study can be extended to other networks, only one network is considered which is the Control Area Network (CAN).

The failure definition for control systems depends on the actual output of the system. It is no more a simple Boolean function on the components failure. Hence the need to calculate the actual output at each step to decide if the system is in a failure situation or not. This remark makes the analysis very difficult for networked control system, especially when delays are variable.

General approaches to reliability evaluation and fault tolerant design may resort to analytical methods or to experimental evaluations. In the latter case, fault injection techniques are commonly adopted, and can basically be grouped into simulation-based or hardware-based techniques.

The purpose adopted in this paper is to use a simulation model for the entire system. First we model all the components of the traditional control system (sensor, controller, actuator, process…), second we model the behaviour of the network and the whole system, and since we are interested on the transmission faults, we will add a third model to inject perturbations on the network. The goal is to evaluate the impact of external perturbations and additional traffic on the reliability of networked control systems.

At the beginning of this paragraph we have shown why traditional technique like fault trees and block diagrams can't be used for our approach. A solution is to recur to dynamic models. The most popular are Markov processes, because of their numerous nice mathematical properties. In practice, the direct use of Markov processes has virtually been given up, to be replaced by some higher level formalisms (i.e. Stochastic Petri Net) that enable the automatic generation of a (potentially huge) Markov process. In our work we decided to work with Petri Net extensions which are well known in reliability evaluation, and well adapted to model complex systems [START_REF] Malhotra | Dependability Modelling Using Petri-Nets[END_REF].After the modelling phase variables to evaluate are declared by means of reward function [START_REF] Malhotra | Dependability Modelling Using Petri-Nets[END_REF].

STOCHASTIC ACTIVITY NETWORK

Stochastic activity networks (SANs) are a stochastic generalization of Petri nets. These models permit the representation of concurrency timeliness, fault-tolerance and degradable performance in a single model. They conserve all the modelling power of Petri nets, and in the same time give the possibility for a compact representation. Structurally, they consist of activities, places, input gates, and output gates. Activities which are similar to transitions in normal Petri nets are of two types: timed and instantaneous. Timed activities represent activities of the modelled system whose duration impact the performance of the modelled system. Instantaneous activities, on the other hand, represent system activities which occur immediately. Input gates and output gates control the enabling of activities and define the marking changes that will occur when an activity completes. SAN models have been used to evaluate a wide range of systems and are supported by several powerful modelling tools such as UltraSAN and Möbius. SAN is defined with the express purpose of facilitating unified Performance/dependability evaluation as well as more traditional performance and dependability evaluation. Dependability evaluation is performed by defining a set of measures in the model.

Table 1. Graphical notations of the elements

NETWORKS

A model is always a compromise between faithfulness and simplicity. A model can be very faithful and represents all the details, but on the other hand it can be intractable. In general assumptions are taken to simplify the model, while leaving it relatively faithful and easily tractable.

Many modelling levels of abstraction can be chosen, with different characteristics and precision. In particular, for network communications we may have: • Bit-accurate models that reproduce the actual clock cycles in the network interfaces.

• Message-accurate models that represent network messages as atomic entities and abstract about detailed timings, while still explicitly representing the protocol implementation. For our model we have chosen the message-accurate models; it is sufficient to our need and can express all the features of the networks (delays, collision, priorities, and errors).

Controller Area Network

CAN is a broadcast bus, with a priority-based access to the medium and non destructive collision resolution. Data to be transferred is encapsulated within communication objects called frame. Each frame contains an identifier (Id), unique to the whole system, which serves two purposes: assigning a priority for the transmission and allowing messages filtering upon reception. The following section presents the CAN behaviour in Fault Scenarios, and some related works.

In [START_REF] Unruh | Error detection analysis of automotive communication networks[END_REF] the authors estimate the expected number of undetected transmission errors during lifetime of a vehicle is lower than 10 -12 . This performance is the result of a very efficient error detection mechanism being used in CAN. This mechanism can be divided into: a message level like CRC code, and a bit level. At the bit level, the transmitter monitors the bus signals and detects errors. Each transmitting station observes the signal on the bus and thus, detects the difference between the bit sent and the bit received. If one error is discovered by at least one station using the above mechanisms, the current transmission is aborted by sending an error flag. This prevents other messages to accept this message. After sending the error flag, the sender automatically re-attempts transmission, and the message reenters a scheduling list, and the one with the highest priority is selected and transmitted on the bus. The number of retransmissions is an important parameter. In the specification of CAN, this number is not defined. To calculate the worst response time, the authors in [START_REF] Cheong | Delay Modelling And Controller design for networked Control systems[END_REF] include the number of retransmission as a variable that characterizes the frame periodic. As cited below, the EMI can affect the medium and produce a transmission error. A detection mechanism is added to cope with these errors. Error recovery mechanisms take some time in detecting and retransmitting the affected message. This lost time is defined as an inaccessibility period where the network isn't ready to provide its service. This behaviour and its consequences are studied in [START_REF] Rufino | A study on the inaccessibility characteristics of the controller area network[END_REF]. A scheduling analysis of CAN is done by (Tindell, et al., 1995). Tindell extended his study to integrate the presence of faults, by adding an additional term error recovery function. The main disadvantage of the analysis is the use of a deterministic model to represent the fault occurrence, which is not realistic. A stochastic fault model which is closer to EMI behaviour and more realistic was proposed by [START_REF] Navet | Worst-Case Deadline Probability in Real-Time Applications Distributed over Controller Area Network[END_REF]. Generalized Poisson Process is used to model the frequency of interference, as well as their duration (single errors and error bursts). In this work the authors introduced for the first time the WCDFP (worst case deadline failure probability) which provides a valuable knowledge on the system's reliability. This information is very important when dealing with a hard real time system where losing a message can lead to a global failure of the whole system.

As it was mentioned early, in this study we assume that transmission errors due to interferences on the medium are always detected and directly proclaimed by the error recovery mechanisms.

We assume also that messages can be lost if one of these conditions is satisfied: 1) the maximum number of retransmission [START_REF] Cheong | Delay Modelling And Controller design for networked Control systems[END_REF] for this message is reached, 2) a new message of the same type is ready to be sent (arrival of a new message always deletes the old one.)

Network model

The Network model is the central model, all other models are linked to it by sharing places (Figure 2). It represents a model of a CAN network. This model takes account of the CSMA/AMP strategy used in CAN that is mainly the management of the priorities. Lines at the left represent station linked to the network. The model can be easily modified to support new components by adding one line to each component. TM place represents the transmission medium, the medium is idle if there is a token inside, and otherwise the medium is busy. The CSMA/AMP inputcase gives the access to the medium to one station; it is done by assigning priority to the places. The activity on_ the_medium is fired to indicate the beginning of a message transmission. If an interference occurred (presence of a token in IN), an error is produced, transmission is aborted, and a token is removed to the error place which will activate the model error_mechanism, this model will send an error frame on the network.

As we mentioned before this model is linked to other models, like sensors, controllers, actuators, and other components. These components will send and receive messages over the network. For a running time, some performance parameters on the network can be evaluated from this model like: -the sum and the average time delays for each node on the network, -the efficiency of the network, defined as the ratio of the total transmitting time to the time used to send messages, including queuing time, blocking time, etc., -network utilization defined by the ratio of the total time used to transmit data and the total running time, -The ratio of lost messages for each node on the network… These parameters were defined by [START_REF] Lian | Performance evalaution of control networks, Ethernet, ControlNet, and DeviceNet[END_REF] to analyze and compare three control networks: Ethernet, ControlNet and DeviceNet. Note that the only considered faults are the perturbations on the network. The fault_injection model injects faults representing an EMI via the medium. We consider that faults have a fixed duration, time between faults is assumed to be exponentially distributed. Since little is known about the frequency of transient faults on local area networks, we will vary their rate occurrence and evaluate their influence. All above reliability for networked control systems can't be evaluated in an analytical way. We propose to evaluate it by Monte-Carlo simulation. The method proposed is described with an example of an NCS with three independent loops that are closed over a control network.

The process controllers, which execute PID algorithms, are well-designed respectively based on different timing parameters. Before the evaluation, we must define the failure definition of a control loop. We based our failure definition on thresholds of the magnitude of the control error, i.e., the difference between the reference output 1 and the actual value. 3 and for two different values of the perturbation fault rate. To understand the results, we evaluate some performance on the network. Table 4 summarizes network parameter in each case. Perturbation on network decrease the network performance, leading to more delay, and high loss probability. For the same value of fault rate results are different, in fact due to the priority assignment in case 2, loop-1 messages are delayed by the higher priority messages of loop_2. The In this paper, we have presented an environment to assess the impact of the transmission errors on the reliability of a NCS. The approach consists in modelling the functional behaviour of classical control system components, and both functional and dysfunctional behaviour of the network. We base our definition of system failure on thresholds of the magnitude of the control error. Thanks to our environment we were able to predict how the transmission errors on the network may affect the system reliability. Special focus was given to the network; all the other components are considered as free fail. One extension will be the integration of faulty components, with the possibility to analyse the consequences on the failure of the system of various sequences of related or unrelated faults. The possibility to introduce variable in the modelling phase and the ease of joining atomic model in Möbius give the possibility to make a library of models ready to be use in other studies. Implementing this library and adding other models is one of our perspectives.

Fig. 1 .

 1 Fig. 1. Networked control system

Fig

  Fig. 2 Network model

  Fig. 3. NCS with three control loops We evaluate the reliability of the first loop by changing two parameters: the priority given to the first loop components and the perturbation fault rate.All the results are obtained by simulation under Möbius tools(Deavours, et al,. 2002) with a confidence interval of 0.1 and a confidence level of 0.95.Case 1: higher priority attributed to loop_1 components.

Figure. 4

 4 Figure. 4 case 1 reliability

Figure

  Figure. 5 case 2 reliability Figure.5 shows the reliability curve for loop_1, under priority assignment given in Table3and for two different values of the perturbation fault rate. To understand the results, we evaluate some performance on the network. Table4summarizes network parameter in each case. Perturbation on network decrease the network performance, leading to more delay, and high loss probability. For the same value of fault rate results are different, in fact due to the priority assignment in case 2, loop-1 messages are delayed by the higher priority messages of loop_2.

Table 2 . Priority assignment

 2

	component	period	Priority
	Sensor1	0.01	1
	Controller1 Event-driven	2
	Sensor2	0.011	3
	Controller2 Event-driven	4
	Sensor3	0.012	5
	Controller3 Event-driven	6
	Figure.4 shows the reliability curve for loop_1, under priority
	assignment given in Table 2, and for two different values for
	perturbation fault rate.

1 Outputs without perturbations on the network are called reference outputs.

Table 3 . priority assignment

 3

	component	period	priority
	Sensor1	0.01	3
	Controller1 Event-driven	4
	Sensor2	0.011	1
	Controller2 Event-driven	2
	Sensor3	0.012	5
	Controller3 Event-driven	6

Table 4 . Network performance

 4

		Average	efficiency utilization	Lost
		delay			probability
		Loop_1
	Case1	λ tf =75 0.0013 0.6670	0.873	0.0041
		λ tf =100 0.0014 0.6181	0.9706	0.0071
	Case2 λ tf =75 0.00132 0.6611	0.8738	0.0044
		6. CONCLUSIONS


Abdelhak Mkhida
email: mkhida@ensam.ac.ma

Jean-Marc Thiriet
email: jean-marc.thiriet@ujf-grenoble.fr

Jean-François Aubry
email: jean-françois.aubry@ensem-inpl-fr

Toward an Intelligent Distributed Safety Instrumented Systems: Dependability Evaluation

In this paper, the modelling and thus the performance evaluation relating to the dependability are studied for structures which have intelligence in the instruments constituting the Safety Instrumented Systems (SIS) in order to determine the contribution of the intelligent instruments in the safety applications. Dynamic approach using Stochastic Petri Nets (SPN) is proposed and the metrics used for the evaluation of the dependability of the Intelligent Distributed Safety Instrumented Systems (IDSIS) refer to two modes of failures mentioned by the safety standards: mode of dangerous failure and mode of safe failure.

INTRODUCTION

Safety Instrumented Systems (SIS) have been used for many years to perform safety instrumented functions in the process industries (IEC61511 2003). The SIS are defined as the collection of all safety related sensing elements to determine an emergency situation, all safety related logic solvers to determine what action to take and all safety related final elements to implement the action [START_REF] Knegtering | Safety lifecycle management in the process industries: the development of a qualitative safety-related information analysis technique[END_REF].

Adding a network communication and intelligent instruments such as intelligent sensors and actuators to the SIS makes this system intelligent. These instruments with a network communication provide several advantages such as integration of diagnosis, cabling reduction, increasing of reconfigurability etc [START_REF] Lian | Network design consideration for distributed control systems[END_REF].

The problem is to quantify the contribution of the use of the intelligent instruments in the safety loops in compliance with related standards. In safety systems many studies were developed for dependability evaluation. Theses approaches ignore the use of intelligent instruments.

The dependability evaluation of such systems is a difficult task and can concern two approaches: a static approach, and a dynamic approach which takes into account the progressive evolution of system states and functioning modes [START_REF] Mkhida | Influence of the control strategy choice on the safety level of the distributed control system[END_REF] [START_REF] Dutuit | Dependability modelling and evaluation by using stochastic Petri nets : application to two test cases[END_REF]. The potential application of such systems is very widespread from electronic trippers to some critical real-time systems such as X-by wired cars, drones [START_REF] Berbra | Diagnosis of a safe networked quadrotor[END_REF]

or aeronautics…

In this paper, we are interested by estimate the metrics that are used in the safety context according to safety standards such as IEC 61508 and IEC 61511. These metrics are the probability of failure on demand (PFD) and the probability of fail safe (PFS).

The methodology used for the dependability evaluation for these systems consists on the structuring and the modelling of these systems; the purpose is to make the verification and analysis by mean of stochastic Petri nets in order to exploit the models.

Modelling is achieved by a stochastic approach using the Stochastic Activity Network (SAN). The SAN is a powerful formalism of the stochastic Petri nets (Moghavar and Meyer, 1984) [START_REF] Ghostine | Influence of transmission faults on the dependability level of networked control systems: application to a control loop[END_REF].

The following sections describe the intelligent distributed safety instrumented systems, how they are modelled and the procedures proposed for their evaluation. The two metrics related to safety performance (PFD & PFS) are calculated. The results are compared, presented and conclusions are drawn.

INTELLIGENT INSTRUMENT

Intelligent Distributed Control Systems (IDCS)

Intelligent or smart instruments have been known for more than two decades. These instruments are more sophisticated than traditional instruments [START_REF] Mekid | Further structural intelligence for sensors cluster technology in manufacturing[END_REF]).

An intelligent instrument is a component part of Intelligent Distributed Control Systems. These systems are thus composed of smart sensors and actuators with calculations capabilities and communication interfaces. Traditional control systems are connected by analogue current (4-20 mA) and voltage loops whereas IDCSs are connect ed by a communication network (generally a fieldbus).

The particularity of a distributed system relates to data exchange between the devices via a communication medium supposed to be a network or a fieldbus. Thereby, the intelligent distributed control system becomes more and more complex and sophisticated, and consequently, makes the design step more difficult [START_REF] Cauffriez | Design of intelligent distributed control systems: a dependability point of view[END_REF].

Structure and functions of an intelligent instrument

The intelligent instruments offer the possibility of a local processing of the information which is distributed on the various entities thus allowing a distribution of the execution of the tasks and appearing a distributed control.

The intelligent instruments are basic instruments that contain microprocessors [START_REF] Nobes | Smart instruments in protective measures, Is your product safe[END_REF]), these instruments meet the following criteria:

-the main purpose of the instrument is to measure or directly control a single process variable, -these instruments include some flexibility in their use due to parameters that are set by the manufacturer or operator, -intelligent instruments are not restricted to measurements but also include actuators (valves, motors) and other control equipment.

Intelligent instruments offer considerable lifetime cost reductions. These cover the whole life cycle, including specification, installation, operation, maintenance, decommissioning…

Let us now consider an architecture and functions of intelligent sensor. The architecture comprises the following [START_REF] Mekid | Further structural intelligence for sensors cluster technology in manufacturing[END_REF]):

-A sensing element that links the external world to a sensor system by generating electrical signal (e.g. voltage, current) with response to physical properties of the environment, -An interface element for signal conditioning and data conversion. The signal obtained from the sensing element is modified and converted to a discrete time digital data, -A processing element that includes a microcontroller with an associated memory and software, -A communication element, which provides a two-way communication between the processing and users.

Various functionalities have been suggested for intelligent sensors. [START_REF] Robert | Smart sensors in flexible manufacturing systems[END_REF] proposed that the intelligent sensors should contain configuration, communication, measurement and validation functionalities. [START_REF] Meijer | Concepts and focus point for intelligent sensor systems[END_REF] includes three functionalities: compensation, computing and communication. [START_REF] Tian | A fieldbusbased intelligent sensor[END_REF] suggested the functions of compensation, validation, data-fusion, and communication.

While [START_REF] Mekid | Further structural intelligence for sensors cluster technology in manufacturing[END_REF] propose that what is called an intelligent sensor should have the functions of compensation, processing, communication, validation, integration and datafusion.

We propose for the generic intelligent instrument [START_REF] Mkhida | Modélisation formelle d'un instrument intelligent dans le cadre d'analyse de sûreté de fonctionnement[END_REF] the functions of measurement, configuration, test, validation, control and communication. Figure 1 illustrates these functionalities.

SAFETY INSTRUMENTED SYSTEMS

Safety instrumented systems are used in many industrial processes to reduce the consequences of process demands on humans, the environment and material assets. Different standards can be used to design safety instrumented systems for process industry like ISA S84, IEC 61508 and IEC 61511 (IEC61508 2000, IEC61511 2003, ISA84 1996). These standards have been developed to ensure that the SIS is designed, implemented and operated according to the specified needs.

Objectives of a safety instrumented system

The IEC 61511 defines a safety function as : "function to be instrumented by a SIS, other technology safety-related system or external risk, which is intended to achieve or maintain a safe state for the process, with respect to a specific hazardous event". The Safety Instrumented Function (SIF) is used to describe the safety functions implemented by instrumented technology. The SIS is the physical system implementing one or more SIFs.

The objectives of SIS is to reduce the frequency at which hazard may occur to an acceptable level (Wiegerinck 2002). The safety function only reduces the risk (multiplication: probability x consequences) and never completely eliminates the risk. Some safety functions do not reduce the probability of the consequences.

All combined instrumentation, devices, and equipment that are fulfilling an intended safety function are considered to be part of the safety instrumented systems. The SIS could be composed by a set of safety-related sensing elements, all safety-related logic solver and all safety-related final elements.

Safety Integrity levels (SIL)

Once the required level of risk reduction to be achieved by the SIS is established, expressed as risk reduction factor(RRF), this level can be translated into the required Safety Integrity Level (SIL) [START_REF] Knegtering | Safety lifecycle management in the process industries: the development of a qualitative safety-related information analysis technique[END_REF].

IEC 61508 and IEC 61511 split the safety integrity into hardware safety integrity and systematic safety integrity. The safety integrity is split into four discrete safety integrity levels (SILs), SIL1 to SIL4. Each SIL represents a maximum allowed probability of failure on demand of the SIS. SIL4 is the level with the most stringent requirements. To fulfil a specific SIL, the SIL requirements related to hardware safety integrity as well as systematic safety integrity must be met [START_REF] Lundteigen | Assessment of hardware safety[END_REF].

To calculate the SIL of a safety function it is required that the complete safety loop from sensor to actuator is considered. Therefore, it is not sufficient to only analyze one subsystem of the safety process, such as the logic solver, and determine the realized SIL.

INTELLIGENT SAFETY

The intelligent safety loops (or intelligent distributed safety instrumented systems) enable us to implement safer plants by using the intelligent devices for safety application in compliance with the new IEC 61511 standard.

The concept of intelligent safety is inherent in the use of intelligence within safety instrumented systems. Several manufacturers claim the certification of intelligent products in the applications of safety while launching "smart SIS" but without real justification to the level of the reliability performances. Researchers employ also the term "intelligent safety" in various applications covering nuclear power [START_REF] Yang | Assessment of the MDNBR enhancement methodologies for the SMART control rods banks withdrawal event[END_REF][START_REF] Chung | Thermal hydraulic analysis of SMART for heat removal transients by a secondary system[END_REF][START_REF] Bae | Safety evaluation of the inherent and passive safety features of the smart design[END_REF], railway or transport systems (Farritor and Goddar 2004).

The use of the intelligent instruments in the process industry was facilitated by the increase of the performances in the microprocessors used in industry, in particular in instrumentation. The justification of the use of these instruments in the applications of safety is not completely proven. These instruments have important benefits useful to this type of applications [START_REF] Nobes | Smart instruments in protective measures, Is your product safe[END_REF]).

There is a trend toward the use of intelligent instruments in safety applications. The inherent ability to diagnose failure is the primary reason. The ability of intelligent instruments to reliably measure complex parameters is another.

The incorporation of the intelligence in the SIS leads towards safety intelligence expressed by IDSIS (Intelligent Distributed Safety Instrumented Systems) whose methodology of evaluation is proposed in the following section.

MODELLING OF THE SYSTEM BEHAVIOUR

Modelling approach

The combinatory methods (fault tree, event trees, reliability diagrams) inherent to reliability study of dynamic systems allow us to identify and evaluate the combinations of events leading to the occurrence of another event [START_REF] Dutuit | Dependability modelling and evaluation by using stochastic Petri nets : application to two test cases[END_REF]. Such combinations do not take into account the order of the events and they exclude any taking into account of time dependency or delays, which can be of high importance in dynamic systems.

The dynamic approach using the Stochastic Petri Nets (SPN or SAN, Stochastic Activity Nets) is proposed to overcome the difficulties mentioned above.

Petri nets, as tools for discrete event simulation are of large use in dependability evaluation. Dynamic changes in the PNs are induced by transition firing. Firing of one or of several transitions changes the net marking. In other terms, it induces a discrete change of states. Thus, dynamic properties of PNs such as parallel firing, successive firing or firing of concurrent transitions may be used to simulate complex events sequences [START_REF] Vernez | Perspectives in the use of coloured Petri nets for risk analysis and accident modelling[END_REF]. The modelling problem consists in achieving the evaluation of dynamic systems which evolve as a function of the time. SPN give the possibility to model a system in which cooperate discrete-time variables and discrete events, which can occur on stochastic basis. The formalism used for the behavioural modelling will naturally be able to represent these characteristics. The use of this formalism will allow the possibility to represent component failure in interaction with the functional behaviour.

In the present approach, it is assumed that the failure distributions of individual components of a system are given, and the dependability measures of the stochastic system are sought. Furthermore, the system is assumed to be dynamic (its properties change with time).

Möbius tool support of SAN

SAN models have been used to evaluate a wide range of systems and are supported by several modelling tools such as Möbius [START_REF] Deavours | The Möbius framework and its implementation[END_REF]. The models are developed using this SAN-based tool. Möbius tool consists on the description of the net structure and of the desired performance variables and solutions methods to be used in the evaluation.

System modelling

In the modelling of the system, the functional and dysfunctional aspects coexist; the failures are divided into safe failures and dangerous failures. A dangerous failure results in an absence of reaction of the safety function. A safe failure results by the setting in a safe position of the system or in an unexpected execution of the safety function. The detection of a safe or dangerous failure, results in a setting into safe position of the system or a forced execution of the safety function. At the beginning, we present SIS without incorporation of intelligent instruments nor communication network. The model of the logic solver is shown in figure 3. The model of the logic solver presented figure 3 shows a disposition of two parts, one functional and the other dysfunctional. In the functional part (on the left), the cycles of the logic solver are carried out by a periodic clock (periode). The self-tests of the various devices are managed locally according to a policy allocating the same duration of test for the various devices and starting with the test of the sensor (Tc), then the actuator (Tv) and finally the logic solver (Ta). This policy is not the only possible and other policies can be possibly established. For the dysfunctional part, it should be made sure that the token is carried of functional part when the system fails safe or dangerous. The logic solver can be also restored in the event of safe failure and it also has a coverage rate of diagnosis which is inherent to him. Now, we will be interested in a model of the system with a communication network (CAN: Control Area Network) and then we will also introduce intelligent instruments instead of traditional instruments. Indeed, CAN protocol is used as a communication architecture for safety critical applications [START_REF] Carvalho | A time-triggered controller area network platform with essentially distributed clock synchronisation[END_REF]. In this article, we focus on CAN while in this study can be extended to other types of communication network. [START_REF] Ghostine | Influence of transmission faults on the dependability level of networked control systems: application to a control loop[END_REF]. The network is made up in this example of three transmitting stations which will be able to send messages on a shared medium (transmission channel). The messages sent are placed in buffer places which are extended places, messages then await the release of the channel to be transmitted towards their destination. The channel access is based on messages priorities. This means that each subscriber will be affected by a priority allowing him or not to send these characteristics via the channel, to avoid the collisions on the channel. The medium (channel) is affected of a delay representing the time of transmission of the messages. The output message of the channel is ready to be sent and it is available in the place ″ received ″.

In the intelligent sensor model (see figure 5), the order of test is not transmitted any more via communication network by the logic solver and the functionality relating to the test as well as the other functionalities of the intelligent sensor is treated locally. The modules of test can be requested either when it is necessary, or in a cyclic way, or permanently according to an active monitoring principle (the intelligent sensor has at one's disposal a self-tests). In our model, the intelligent sensor has self-tests in a cyclic way according to the instants of period. The other functionalities are also synchronized by this period.

The model of the actuator does not differ much from that of the intelligent sensor and the dysfunctional parts are similar. It should be noted that the model of the sensor, the actuator and the logic solver are sampled models and information is sent to the communication network periodically. The coverage rate of diagnosis is taken in this example equal to 60% for the whole of the devices. The period of the selftests operated by the logic solver is selected equal to one hour. I.e., that in compliance with the policy of test chosen, the duration of cycle of the logic solver is 3 hours Also let us note that this duration affects the PFD and PFS when it is changed. For example, for a value of period of test of 1 hour such as it is presented in figure 7, the PFD is 0.113 and the PFS is 0.1604 for a time of 8760 hours mission. While these values apply to one duration of test of the logic solver equal to 2 hours, the PFD passes to 0.107 and the PFS becomes equal to 0.165. It is to show the influence of the frequency of the period of the test of the logic solver on the total performances in safety of the system. Now, we will be interested in the performance evaluation in safety after the introduction of the intelligent instruments. The results are given for a duration of simulation equal to 8760 hours. The evolution of the two metrics (see figure 7) which describe the performances in safety of the studied system show a clear improvement of the probability of dangerous failure according to the addition of intelligence and a weak degradation of the probability of safe failures. Indeed, more there are a high rate of detection by the integrated means of self-tests, more the failure rate weakens and the safe failure rate increases by the means of the dangerous failures which were transformed into some kind of safe failures.

CONCLUSION

This paper presents the dependability evaluation of intelligent distributed safety instrumented systems. The performance of these systems is evaluated and compared to the systems without intelligence. The results express the safety performance of these systems, using the particular modelling assumptions made, and the reliability data chosen. Some parameters have a significant impact on the PFD and PFS values, but the accent was put on the contribution of intelligence on the SIS. The setting of the intelligence in SIS improves the PFD but not in a significant way and degrades the PFS slightly. We can say this is the first level of our studies, we have to emphasise the modelling of intelligence in the future model, in order to have a more realistic representation of it.

Special focus was given to the intelligent instruments introduced in the SIS with their failures without counting the failure of network that can influence the performance of safety.

Fig. 1 .

 1 Fig. 1. Functional architecture of an intelligent instrument

Figure 2 :

 2 Figure 2 : Example of sensor modelling and its failures

Figure 3 :

 3 Figure 3: Structure of Logic solver

Figure 4 :

 4 Figure 4 : Network model This submodel represents the behaviour of the CAN network[START_REF] Ghostine | Influence of transmission faults on the dependability level of networked control systems: application to a control loop[END_REF]. The network is made up in this example of three transmitting stations which will be able to send messages on a shared medium (transmission channel). The messages sent are placed in buffer places which are extended places, messages then await the release of the channel to be transmitted towards their destination. The channel access is based on messages priorities. This means that each subscriber will be affected by a priority allowing him or not to send these characteristics via the channel, to avoid the collisions on the channel. The medium (channel) is affected of a delay representing the time of transmission of the messages. The output message of the channel is ready to be sent and it is available in the place ″ received ″.

  Figure 5: intelligent sensor model 6. RESULTS OF SIMULATIONThe procedure used for the calculus of probability of dangerous failures PFD and the probability of safe failures PFS consists of the presence of tokens in the places which describe the safe failures and the dangerous failures for the whole system.

Figure 6 :

 6 Figure 6: Evolution of PFD and PFS according to timeFigure6shows the evolution of PFD and PFS for a 10000 hours duration which is a little higher than one year (8760 hours).

Figure 7 :

 7 Figure 7: Safety performances for SIS and IDSIS


Hoa Van

Nguyen
email: vhnguyen@irisa.fr

Dominique Lavenier
email: lavenier@irisa.fr

Speeding up Subset Seed Algorithm for Intensive Protein Sequence Comparison

Sequence similarity search is a common and repeated task in molecular biology. The rapid growth of genomic databases leads to the need of speeding up the treatment of this task. In this paper, we present a subset seed algorithm for intensive protein sequence comparison. We have accelerated this algorithm by using indexing technique and fine grained parallelism of GPU and SIMD instructions. We have implemented two programs: iBLASTP, iTBLASTN. The GPU (SIMD) implementation of the two programs achieves a speed up ranging from 5.5 to 10 (4 to 5.6) compared to the BLASTP and TBLASTN of the BLAST program family, with comparable sensitivity.

I. INTRODUCTION

The sequence similarity search finds similar segments, or local similarities, between two DNA sequences or protein sequences, measured by match, mismatch and gap scores. To establish similarity, an alignment score is calculated. Its objective is to locate similar regions in the DNA or protein sequences because high degree of similarity often implies similar function or structure. A typical search is to query a bank with a gene whose function is unknown. The results of the request correspond to similar segments.

Recent biotechnology improvements in the sequencing area lead to an increasing amount of genomic databases. Genbank [1], for example, contains more than 180 billion of nucleotides (November 2007) and its size is multiplied by a factor ranging from 1.4 to 1.5 every year. Thus, the similarity search is a very time consuming task which may need the use of supercomputers.

Similarity search is primarily a data content search. Banks are systematically scanned from the first to the last sequence. It is a challenging task requiring a huge computing power. With the stagnation of microprocessor frequencies, there is a great interest in optimizing this task.

There are several algorithms to search alignments. One of the first is the Smith-Waterman algorithm developed in 1981 [START_REF] Smith | Identification of common molecular subsequences[END_REF]. It uses dynamic programming techniques and has a quadratic complexity. Latter ones, developed in 1990 [START_REF] Altschul | Gapped BLAST and PSI-BLAST : A new generation of protein database search programs[END_REF] [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF], as the BLAST program family, are based on heuristic. The key of the BLAST heuristic is that a statistically significant alignment is likely to contain a common word of W characters. Thus, instead of scanning the whole search space, the search only focus to these specific common words, called anchor points. Consequently, the computation time is drastically reduced. However, the heuristic does not guarantee to detect the best alignments, but it reports most of the significant alignments.

More recently, several works have proposed to use spaced seeds to perform detection of anchor point, including Part-ternHunter [START_REF] Ma | PatternHunter: faster and more sensitive homology search[END_REF], PartternHunter II [START_REF] Li | PatternHunter II: Highly sensitive and fast homology search[END_REF]. Instead of considering a word of W contiguous characters as a seed, a word of no necessarily consecutive W characters may be considered. More precisely, a spaced seed looks like a binary mask string, such as 10101, where matches in all positions 1 are required for an anchor point; 0 denotes allowed mismatches. The next seed generation doesn't consider one seed but a set of several seeds which can be of different length and possibly made of different spaced seeds (multiple seeds) [START_REF] Daniel | Optimizing multiple seeds for protein search[END_REF]. Another seed family (subset seed) groups different characters in the same set [START_REF] Kucherov | A unifying framework for seed sensitivity and its application to subset seeds[END_REF], still providing better expressiveness.

This paper presents an effective implementation of a recently proposed seed-based heuristic, called subset seed, that exploits fine-grained parallelism using data indexing for intensive protein sequence comparison, such as the comparison of two protein banks or the comparison of a protein bank and a genome translated to its six reading frames. We use Intel SSE2 SIMD instructions and the SIMD architecture of GPU (Graphics Processing Unit) device to speed up critical section of the algorithm.

The remainder of this paper is organized as follows. In section 2, a brief background of the subset seed algorithm is given. Section 3 describes how to use Intel SSE2 SIMD instructions and SIMD architecture of GPU for implementing this algorithm. In section 4, we present numerical test results before section 5 concluding and presenting future work.

II. SIMILARITY SEARCH

In this section we present the subset seed similarity search algorihtm, including seed-based heuristic.

A. Seed based similarity search

The idea of heuristic search is based on an observation: In general, significant alignments have a region of high similarity which is characterized by a zone of W identical charactersword -between two sequences that make up this alignment. Therefore, first, we can locate places where these words appear. When the same word is located on two sequences, it provides an anchor point for calculating alignment by 57 978-1-4244-2379-8/08/$25.00 (c)2008 IEEE extending to the right and to the left the number of match between identical symbols. If the extension finds a significant similarity (measured in the number of match / mismatch), a final phase is processed to build a better alignment by considering the inclusion or deletion of symbols (gap).

The scheme below shows a construction example of alignment between two protein sequences. The word ARV, in both sequences, serves as an anchor point (stage 1). Starting from this anchor point, we extend to the right and to the left for getting an interesting match / mismatch ratio (stage 2). Finally, we extend the alignment for taking into account the errors of gap (stage 3).

stage #1 K V I T A R V T G S A Q W C D N T G V K N I H M | | | K L I S A R V K G S Q F C T N P T G M K A N I H stage #2
This is a 3-stage process: (1) anchor point detection; (2) ungapped alignment; (3) gapped alignment. This heuristic is implemented with many variants in a number of programs such as programs of the BLAST family which are the reference programs in bioinformatics. Depending on the type of database inputs, each BLAST [START_REF] Altschul | Gapped BLAST and PSI-BLAST : A new generation of protein database search programs[END_REF] [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF] program has a different name: BLASTP when the query is an amino acid sequence and the database is a protein bank, TBLASTN when the query is an amino acid sequence and the database is a nucleotide bank.

In intensive protein sequence comparison between two banks (bank1 and bank2), BLASTP executes comparison of sequence S (bank1) against the set of sequences in bank2. It first searches for all hits between sequence S and bank2 by using an index structure such as a lookup table. Then, these hits are extended, first without gaps, then allowing them. Hits meeting specified threshold are returned as statistically significant alignments.

The advantage of BLASTP is its speed compared to dynamic programming techniques. The speed depends on the hit size. Actually, longer the hit size, lower the possibility to locate anchor point.

B. Subset seed similarity search

In the context of intensive protein sequence computation of two banks: bank1 (N sequences) and bank2, BLASTP must scan and process bank2 N times. To avoid this drawback, we dynamically index the two banks in main memory. Based on this idea, the two banks are scanned only once just for making the index.

The advantage of a double indexing is that the systematical scan of bank2 disappears. Through an appropriate indexing structure, we can directly point to all identical words in the two banks. If a W-AA word appears |nb1| times in bank1 and |nb2| times in bank2, there are |nb1|×|nb2| hits. It means that there are also |nb1| × |nb2| ungapped alignments to calculate. The subset seed algorithm (iBLASTP) proceeds in 5 successive stages as shown on TABLE I. Stage 0 index the two banks, stage 1 constructs two neighboring blocks, stage 2 performs ungapped alignments, stage 3 computes gapped alignments and stage 4 displays alignments. We describe in detail these stages in the following subsections. Stage 0: Bank indexing. Spaced seeds [START_REF] Li | PatternHunter II: Highly sensitive and fast homology search[END_REF] and their relative, vector seeds [START_REF] Daniel | Optimizing multiple seeds for protein search[END_REF], can augment selectivity and reduce the false positive rate in the seeding step. But spaced seeds can not distinguish between different types of mismatches in the case of protein sequence. Because all amino-acids are not equivalent, we can create groups of characters. Based on the groups, [START_REF] Kucherov | A unifying framework for seed sensitivity and its application to subset seeds[END_REF] [7] created subset seeds, the extension of spaced seeds. Subset seeds have an intermediate expressiveness between spaced and vector seeds but they allow an efficient hashing method to locate hit. The main advantage of subset seeds is that they provide a powerful seed definition and, in the same time, preserves the possibility of direct indexing. A subset seed of W characters (see detail in [START_REF] Peterlongo | Protein similarity search with subset seeds on a dedicated reconfigurable hardware[END_REF]) is defined as a word (s 1 s 2 ...s W ). Figure 1 presents a subset seed example of 4 characters. To build indexing structure for bank1, we store all sequences into memory. Then each sequence is parsed into fixed length overlapping words of W characters -W-AA words. Each W-AA word is matched in subset seed for finding a corresponding word. For example, to apply the subset seed of Figure 1, the words AQAA, APAA, ASAA are translated to a similar word AgAA (or a seed). Furthermore, each occurrence of W-AA word has a position (index). We use a buffer integer (size equal to size of bank1) to store all positions of occurrences of similar W-AA words (translated), for example the words AQAA, APAA, ASAA. The positions of similar W-AA words (ex. AgAA) are considered as a list of W-AA words. We have 20 W (20 is the number of amino acids) lists corresponding to 58 978-1-4244-2379-8/08/$25.00 (c)2008 IEEE 20 W W-AA words. The number of linked lists depends on the subset seed being used. For example, if W = 4, there are 20 4 words: AAAA, ..., YYYY and when applying a subset seed such as in Figure 1, there are 8800 (20 × 2 × 11 × 20) linked lists. We use an other buffer to store the heads of the linked lists.

Bank2 is processed in the same way. At the end of this stage, for each W-AA word (in subset seed), we have two lists which content |nb1| occurrence positions of W-AA word in bank1 and |nb2| occurrence positions of W-AA word in bank2. Stage 1: Construct two neighboring blocks. Hits detected in the first stage are considered as the base of ungapped alignments. The objective is to be able to rapidly decide if one hit has favorable environment to build an alignment. The indexing structure above gives us the lists of positions of W-AA words and their neighboring knowledge for performing ungapped alignments. In this stage, for each W-AA word, we construct two neighboring blocks: one block (nb1) for bank1 and another block (nb2) for bank2. To build neighboring block nb1, index1 provides a list of all W-AA word positions, then for each position a subsequence is copied. Each subsequence contents W-AA word and its neighboring knowledge. Two blocks of subsequences allow us to reuse cache memory when we perform ungapped alignment. Stage 2: Ungapped alignments. For each W-AA word, we perform |nb1|×|nb2| ungapped alignments of two neighboring blocks. More precisely, for each subsequence in block nb1, we calculate the ungapped score with |nb2| subsequences in block nb2. Ungapped extension is performed on subsequences of bounded length, and ignores insertion and deletion events. The aim is to quickly compute a score according to a substitutionscore matrix. Starting from 0, the score is increased or decreased depending on the substitution-score values in this matrix. If ungapped extension score is superior than a specified threshold S 1 , this ungapped alignment is passed to the next stage. Stage 3: Gapped alignments. We use the dynamic programming technique to build gapped alignments. This stage is divided into two sub-stages: The first computes small gapped alignments (stage 3.1), it aims to limit the searching space of dynamic programming. In this procedure, there is a constraint on the number of allowed gap. If score exceeds a specified threshold S 2 , the standard dynamic programming procedure is launched (stage 3-2).

A gapped alignment contents at least one ungapped alignment. If a gapped alignment contents two ungapped alignments, they are not always on the same diagonal. For this reason, sometime, we must recalculate gapped alignments for different ungapped alignments. We can avoid this problem by storing all ungapped alignments contented in the gapped alignments in a list, and performing gapped alignment only for hit which doesn't belong to this list. stage 4: Traceback & display. In this stage, trace-back information optimizes alignments recorded in previous stage and displays them to user.

The structure of our approach has been developed with objective of exploiting fine grained parallelism on specific hardware structure. Thus, the size of the subsequences in stage 2 and 3-1 is fixed.

C. Profiling of the code

We have implemented two sequential programs: iBLASTP and iTBLASTN by reference to the programs BLASTP and TBLASTN programs of the BLAST family. The first compares two protein banks and the second compares a protein bank with a genome translated into 6 reading frames (6 ways to match a DNA to proteins). Both versions run on standard processors. They allow us to measure more precisely the speed-up of fine grained parallelism.

We have examined the average runtime for each stage. Results are shown in TABLE II. We found that the second and third (3-1) stages consumed respectively 60-72% and 23-25% of the total execution time. The processing time of ungapped and small gapped alignment takes most of the time. Obviously, to get a significant speed up, execution time of these stages needs to be reduced. To increase the performance of these stages, we exploit two types of fine-grained parallelism: SIMD instructions of modern microprocessors and SIMD architecture of recent GPU.

A. SIMD instructions

With introduction of MMX microprocessors in 1997, Intel made computing with SIMD technology available in a generalpurpose microprocessor. Multimedia extensions exploit finegrained parallelism, where computations are split into subwords with independent units operating on them simultaneously. SIMD processing was enhanced with the addition of the SSE2 (Streaming SIMD Extension) in the Pentium 4 microprocessor. It permits handing sixteen simultaneous byte operations in 128-bit XMM registers.

1) Ungapped alignments: Instead of performing |nb1| × |nb2| ungapped alignments for a W-AA word by using the Single Instruction Single Data (SISD), we have used SIMD instructions. In this case, a score fits into 1 or 2 bytes instead of 4 bytes.

The SIMD register simultaneously contains k scores related to k subsequences of block nb1. When k subsequences are computed with all subsequences of block nb2, a simple speed improvement is achieved by creating a kind of profile score for the k subsequences. This profile score is a specific substitutionscore matrix of k subsequences. Instead of indexing the original substitution-score matrix by subsequences symbols, the 59 978-1-4244-2379-8/08/$25.00 (c)2008 IEEE new matrix is indexed by the position on the k subsequences and a subsequence symbol (block nb2).

This profile score permits to load substitution-scores in a single read operation. When using 1 byte integer k = 16 and when using 2 byte integer k = 8. Each element of the SIMD registers maps one subsequence. The first element maps s1 1, the second element maps s1 2, still the last element maps s1 k. Thus, score(i) (between k sequences in block nb1 and one subsequence of block nb2) is considered as a score result between the i-th subsequence and the subsequence in block nb2.

When k = 16, the 128-bit wide registers are divided into 16 8-bit elements. Dividing the registers into 8-bit elements limits the range of scores from 0 to 255. In mots cases, the scores fit in the 8-bit range unless the subsequences are totaly similar.

block nb1 block nb2 s1_1: ..K C A G A A S D.. s2_1: ..N M A H A A S Q.. s1_2: ..A G A S A A V N.. s1_3: ..L N A A A A W E.. s1_4: ..Y E A N A A L T.. s1_5: ..M T A S A A K H.. s1_6: ..S H A G A A S Q.. s1_7: ..W Q A D A A T S.. s1_8: ..T S A E A A E M.. The above scheme illustrates an ungapped alignment computation between 8 subsequences in block nb1 with one subsequence in block nb2 when the position on 8 subsequences is at the beginning of seed and subsequence symbol (s1 of black nb2) is "A".

2) Small gapped alignments: In addition, we also use SIMD instructions to compute small gapped alignments. The ungapped extensions passed to the stage 3.1 are stored in a list. When this list contains at least K elements, all elements on this list are considered for small gapped alignment with SIMD instructions.

Unlike ungapped alignment, if there is one ungapped alignment passed equally there is a pair of subsequences. The pair of subsequences in small gapped alignment is usually similar. Furthermore, the length of subsequence is longer than the length of subsequence in ungapped alignment. Consequently, when k =16, scores usually fit out the range of 8-bit. Thus, we decided for a unique value of k = 8.

In addition, we must perform only one small gapped alignment for a pair of subsequences. The computing begins at seed position, extends to the right and then to the left. Thus, we cannot use the subsequence profile because it is just used once. Consequently, in the Smith-Waterman algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF], to compute each cell score, the value of substitution-score must be sequentially calculated. More precisely, symbols from the two subsequences have to be read and look-up into the substitution-score matrix in order to calculate the corresponding substitution-score. This process has to be repeated for the 8 elements in the SIMD register.

The scheme below illustrates small gapped alignment computation for 8 pairs of subsequences. The computing begins at seed position, extends to the right and then to the left. The score of pair(i) is a sum of score(i) of two extentions.

bank1 bank2 s1_1: V..K S A G A A S..V s2_1: N..A D A H A A G..N s1_2: N..R N A S A A F..S s2_2: G..Q Q A G A A F..L s1_3: K..V R A A A A K..M s2_3: T..G Q A N A A K..G s1_4: F..I V A N A A Q..K s2_4: F..Y G A P A A I..T s1_5: A..R I A S G A E..I s2_5: D..L I A K A A E..F s1_6: G..H R A G A A L..F s2_6: C..H L A R A A V..A s1_7: H..D L A D A A E..E s2_7: Y..N F A N A A D..G s1_8: S..R S A E A A Q..A s2_8: I..G F A T A A Q..S

B. SIMD architecture of GPU

During the last decade, GPUs [START_REF]NVIDIA CUDA Compute Unified Device Architecture[END_REF] have been developed as highly specialized processors for the acceleration of raster graphics. The GPU has several advantages over CPU architectures for highly parallel intensive workloads, including higher memory bandwidth, significantly higher floating-point capabilities, and thousands of hardware thread contexts with hundreds of parallel compute pipelines executing programs in a single instruction multiple data (SIMD) mode.

NVIDIA has introduced a new GPU, i.e. Geforce 8800 GTX and a C-language programming called CUDA [START_REF]NVIDIA CUDA Compute Unified Device Architecture[END_REF] (Compute Unified Device Architecture). Geforce 8800 GTX architecture comprises 16 multiprocessors. Each multiprocessor has 8 SPs (streaming processors) for a total of 128 SPs. Each group of 8 SPs shares one L1 data cache. A SP contains a scalar ALU (Arithmetic Logic Unit) and can perform floating point operations. Instructions are executed in a SIMD mode.

The threads on GPU are organized in blocks. The grid of thread blocks is executed on the device. Thread blocks have the same dimensions. A thread block is processed by only one multiprocessor, so that the shared memory space resides in the on-chip shared memory leading to very fast memory accesses. A multiprocessor can process several thread blocks concurrently by partitioning among the set of registers and the shared memory.

For each W-AA word, |nb1| × |nb2| ungapped alignments are performed in parallel by the GPU. The ungapped extensions passed to the stage 3.1 are stored in a list. When this list contains at least K elements, all elements on this list are considered for small gapped alignment on the GPU.

1) Ungapped alignments: We implemented ungapped alignment stage on graphic card by using the matrix multiplication algorithm [START_REF]NVIDIA CUDA Compute Unified Device Architecture[END_REF]. For each W-AA word, there are two subsequence blocks. Suppose that block A[wA, hA] corresponds to block nb1; wA is the length of subsequences, hA is the number of subsequences in block nb1; block B[wB, hB] corresponds to block nb2, wB is the number of subsequences in block nb2, hB is the length of subsequences. On the other hand, block B is the transposition of block nb2. Furthermore, we use an other block C[hA, wB] to store scores of ungapped alignments between block nb1 and block nb2. The value of each cell[i,j] in block C corresponds to the score of subsequence j (row j of block nb1) and subsequence i (column i of block nb2). The task of |nb1| × |nb2| ungapped alignments between block A and block B is split among threads on GPU as followed: each thread block is responsible for computing on square subblock C sub of C. Each thread within the block is responsible for computing one element of C sub (Figure 2). The dimension block size of C sub is chosen equal to 16. Thus, there are hA 16 × wB 16 thread blocks in grid. Threads within thread block share memory each other.

Two subsequence blocks are mapped to the texture memory of the GPU. The texture memory is shared by all the processors, and speed up comes from the texture memory space by being implemented as a read-only region of device memory. At the beginning of the computation, each thread loads a character from the texture memory to shared memory by a texture reference, called texture fetching. 2) Small gapped alignments: The use of high performance computing on GPU is efficient only to perform large tasks. Thus, we use the GPU to execute small gapped alignments when there are at least K elements ready for computation. With K small gapped alignments, there are 2K extensions extending in two directions. To be enable to compute 2K extensions on the GPU, we have to construct two subsequence blocks as already done for ungapped alignment: one block (A) for bank1 and one block (B) for bank2. Compared to ungapped alignment, there is a difference: for one small gapped alignment, we copy two subsequences in bank1 -one at the left and one at the right of seed -to block A, and two subsequences (bank2) to block B. Consequently, there are 2K subsequences in each block. The GPU divides 2K extensions into thread blocks, each thread within a block is responsible for computing one extension.

The form of thread block is 1 × 256. Thus, there are 2K 256 thread blocks. Two subsequence blocks are mapped to the texture memory. At the beginning of the computation, each thread copies its pair of subsequences from the texture memory by the texture reference to its local memory for reducing memory access conflict. The scores of 2K extensions are stored in block C[2K,1] (Figure 3).

IV. RESULTS

We tested the implementations on an Intel 2.6 GHz processor with 2 MB cache L2, and 2 Gb RAM, running Linux (fedora 6). We used the graphic card GeForce 8800 GTX (version GPU). The characteristics of this board are as follows: The execution time is calculated using the Linux command, "time". The BLASTP and TBLASTN of the BLAST family (version 2.2.16, 2007) were launched with default parameters, except for the statistical parameter, E-value, that is chosen at 10 -3 , which is a reasonable value in a context of intensive sequence comparison.

The sensitivity is evaluated in relation to the number of alignments found by the two programs. Specifically, we test whether alignments begin and end at the same places in the two banks with a margin of P % calculated on the average size of the two alignments. For example, to compare two alignments of size 100 at 5 %, we check that the start and stop positions of the alignments are in the range of 5 amino acids.

TABLE III compares execution time of (SIMD/GPU) iBLASTP and BLASTP with parameter set as: PIR-Protein bank and four SWISS-PROT (SWP) banks, With fine grained parallelism of SIMD instruction, the performance of ungapped extension achieve a speed up ranging from 5 to 7 compared to the standard implementation. Acceleration factors for 13 and 16 are obtained for GPU as reported in Figure 4.

Furthermore, the SIMD small gapped alignment achieves a speed up ranging from 2.5 to 2.8 compared to the standard implementation. Acceleration factors for 13 and 19 are obtained for the GPU as reported in Figure 4.

The main disadvantage of SIMD small gapped alignment is that element scores are sequentially accessed from memory. Consequently, speed up is less than speed up of SIMD ungapped alignment. Compared to SIMD ungapped alignment, the GPU can achieve about 2 times faster. The powerful performance of SIMD ungapped alignment is that it uses effectively the subsequence profile cache and element scores can be accessed in parallel. We believe that the GPU lost performance because of several factors, including poor cache memory and lack of high bandwidth access to cached data.

Based on the definition of two equivalent alignments, we compared the sensibility between iBLASTP and BLASTP (iT-BLASTN and TBLASTN). The same data sets as previously are considered. The number of alignments found by the two iBLASTP (iTBLASTN) implementations is the same, thus, we have compared the sensibility between BLASTP and GPU iBLASTP (TBLASTN and GPU iTBLASTN). This sensibility was evaluated by considering three values of P: 2 %, 5 %, and 10 %. TABLES V & VI summarize the results. For each performance, the number of alignments found by BLASTP and GPU iBLASTP (TBLASTN and GPU iTBLASTN) is specified as a percentage of equivalent alignments. The two programs BLASTP and GPU iBLASTP (TBLASTN and GPU iTBLASTN) detect the same number of alignments, and approximately 95 % (96 %) of the alignments are equivalent. In fact, the difference could be explained as followed: BLASTP and GPU iBLASTP (TBLASTN and GPU iTBLASTN) do not use the same seed. Thus, there are some alignments found by BLASTP (TBLASTN) and not by GPU iBLASTP (GPU iTBLASTN) and inversely.

V. CONCLUSION

In this paper we have presented a parallel subset seed algorithm for similarity search between protein sequences by using SIMD instruction and GPU. To our knowledge, this is the first attempt for this type of algorithm on these structures, the other implementations for this search are based on the dynamic programming algorithm.

Two programs have been developed and tested: iBLASTP and iTBLASTN. They refer to BLASTP and TPLASTN of the BLAST program family that are daily used by thousands 62 978-1-4244-2379-8/08/$25.00 (c)2008 IEEE of biologists. For both GPU (SIMD) programs, we have obtained acceleration factors of 5.5 and 10 (4 and 5.6) overall execution time, and gains ranging from 13 to 19 (2.3 to 7) for computation parallelized on GPU (SIMD).

One of the limitations of our current implementation is the sequential treatment of stage 3-2, which limits significantly the overall speed up. The solution may come from works demonstrated in [5] [2]. This stage uses essentially the dynamic programming algorithm which is effective on the GPU and through SIMD instructions. Therefor, we expect to obtain a substantial gain by parallelizing this stage on GPU or using SIMD instruction.

Fig. 1 .

 1 Fig. 1. Subset seed example with length equal to 4

60 978- 1 -

 1 4244-2379-8/08/$25.00 (c)2008 IEEE

Fig. 2 .

 2 Fig. 2. Computing ungapped alignment on GPU.

Fig. 3 .

 3 Fig. 3. Computing small gapped alignment on GPU.

•

  16 multiprocessors SIMD at 675 MHz; each multiprocessor is composed of eight processors running at twice the clock frequency; • maximum number of threads per block: 512; • amount of device memory: 768 MB at 1.8 GHz; • maximum bandwidth observed between the computer memory and the device memory: 2 GB/s. The following banks have been used for testing: (1) a protein bank contenting 141,708 sequences from the PIR-Protein bank (the version 80 01/2005) with an average length of 340; (2) a nucleotide bank contenting 27,360 sequences (gbvrt3 in GenBank, the version 156) with an average length of 5,454; (3) a set of four protein banks (extracted from the SWISS-PROT bank, the version 05/2007) contenting respectively 5,000, 10,000, 20,000, and 40,000 sequences with an average length of 367.

Fig. 4 .

 4 Fig. 4. Speed up of SIMD and GPU for ungapped alignment

Fig. 5 .

 5 Fig. 5. Speed up of SIMD and GPU for small gapped alignment

  TABLE IVcompares execution time of (SIMD/GPU) iTBLASTN and TBLASTN with parameter set as: GenBank and four SWISS-PROT banks. Acceleration factors of 4.2 and 5.6 are respectively obtained for SIMD iBLASTP and SIMD iTBLASTN.

	978-1-4244-2379-8/08/$25.00 (c)2008 IEEE	61

TABLE V

 V ,794 1,062,130 94.8% 95.0% 95.4% 40k 2,237,076 2,237,140 94.5% 95.2% 95.6%

	COMPARISON OF SENSITIVITY BETWEEN IBLASTP AND BLASTP
	nb. seq	nb. of alignments		value of P
	(SWP)	BLASTP	iBLASTP	2%	5%	10%
	5k	305,435	305,663	94.7% 94.9% 95.4%
	10k	611,031	610,993	95.3% 95.3% 95.8%
	20k	1,047

TABLE VI COMPARISON

 VI OF SENSITIVITY BETWEEN ITBLASTN AND TBLASTN

	nb. seq	nb. of alignments		value of P
	(SWP)	TBLASTN	iTBLASTN	2%	5%	10%
	5k	290,016	290,406	95.6% 95.7% 95.9%
	10k	572,608	573,880	96.0% 96.2% 96.4%
	20k	1,172,466	1,173,378	96.4% 96.5% 96.6%
	40k	2,208,330	2,211,935	96.4% 96.6% 96.7%