Maïtine Bergounioux
email: maitine.bergounioux@univ-orleans.fr

Xavier Bonnefond
email: xavier.bonnefond@mip.ups-tlse.fr

Pierre Maréchal
email: pierre.marechal@mip.ups-tlse.fr

Optimal control techniques for thermo-acoustic tomography

Keywords: Tomography, Image Processing, Control Theory

Thermo-acoustic and photo-acoustic tomography are imaging techniques that combine high electromagnetic absorption contrast between two media with ultrasound high resolution. Both techniques lead to an ill-posed inverse problem. One currently has a choice between three main types of reconstruction procedures namely the filtered backprojection formulae, eigenfunction expansion methods and time reversal method. In this paper we propose to investigate this inverse problem with an alternative control formulation : in our model the function to be recovered is the control function while the (acoustic) pressure is the state function which satisfies a wave equation. We stress that our objective is to give modelling hints, and not to provide new results. Otherwise expressed, the originality of this paper lies essentially in the way the problem is stated.

Introduction

Thermo-acoustic and photo-acoustic tomography are imaging techniques that combine high electromagnetic absorption contrast between two media with ultrasound high resolution. These hybrid systems use an electromagnetic pulse as an input and record ultrasound waves as an output. The electromagnetic energy is distributed at a given time as uniformly as possible through the object. The induced increase of temperature depends on the local absorption properties (for example, cancerous tissues absorb more energy than healthy ones). This opens the way to the detection of heterogeneities via measurements of the pressure field. Heterogeneities behave like internal acoustic sources, and the signals recorded by pressure detectors outside the medium under study provide information on their distribution. One speaks of thermo-acoustic tomography (TAT) when the heating is realized by means of microwaves, and of photo-acoustic tomography (PAT) when optical heating is used. While in TAT waves of radio frequency range are used to trigger the ultrasound signal, in the PAT the frequency lies in the visual or near infra-red ranges. In brief, TAT and PAT are two hybrid techniques using electromagnetic waves as an excitation (input) and acoustic waves as an observation (output). For the mathematical purpose there is no distinction between these methods. Both techniques lead to an ill-posed inverse problem of the same form which, under simplifying assumptions, entails inversion (in the wide sense) of the spherical Radon transform.

More precisely, the main problem can be formulated as follows [19,20,21]: given the sound speed c(x) and measured data y obs on S ⊂ R n (n = 2, 3), find the initial value u o (x) of the pressure Figure 1: Realization of a tomograph with integrating transducers, from Patch and Scherzer [29] y(x, t) where y is the solution to the problem

                 ∂ 2 y ∂t 2 (t, x) -c 2 (x)∆y(t, x) = 0, (t, x) ∈ [0, T ] × R n , y(0, x) = u o (x), x ∈ R n ,
∂y ∂t (0, x) = 0, x ∈ R n , y(t, x) = y obs (t, x), x ∈ S, t ∈ [0, T ].

(1)

We shall briefly justify this model in the sequel. The initial value u o is the TAT (or PAT) image. This problem is highly ill-posed. In most reconstruction methods in TAT, additional assumptions are performed such as, for example, conditions on the support of the function to be recovered and/or the observation surface, or a constant sound speed. A nice overview of the state of art has been done in [21]. We summarize the introduction (for more details one can refer to the whole paper). One currently has a choice between three main types of reconstruction procedures for closed observation surfaces, namely the filtered backprojection formulae, eigenfunction expansion methods and time reversal methods.

The filtered backprojection approach is the most popular [13,15,16,19,22]. However, it is not clear that backprojection-type formulae could be written for any closed observation surface S. In [16], inversion formulae are provided assuming odd dimensions and constant sound speed. Indeed, in this case the Huygens' principle holds : for any initial source with a compact support, the wave leaves any bounded domain in a finite time. This is no longer true if the spatial dimension is even and/or the sound speed is not constant. All known formulae of filtered backprojection type assume constant sound speed and thus are not available for acoustically inhomogeneous media. In addition, the only closed bounded surface S for which such formulae are known is a sphere.

Expansion series are useful in the case where the Huyghens principle is valid. This approach was extended to the constant speed and arbitrary closed observation surface and modified by the usage of the eigenfunctions of the Laplacian with Dirichlet conditions on S [8]. It theoretically works for any closed surface and for variable sound speeds [32].

The time reversal method [20,21] can be used for approximating the initial pressure when the sound speed inside the object is variable. It works for arbitrary geometry of the closed observation surface S. The sound speed can be variable. Ammari et al. [3,4,[START_REF] Ammari | Direct reconstruction methods in ultrasound imaging[END_REF][START_REF] Ammari | Photo-acoustic imaging for attenuating acoustic media[END_REF] have performed sharp analysis of these problems both from the modelling and numerical poins of view.

In this paper we propose to investigate this inverse problem with an alternative formulation. We propose a control approach: in our model the function to be recovered is the control function while the pressure is the state function which satisfies a wave equation. We stress that our objective is to give modelling hints, and not to provide new results. Otherwise expressed, the originality of this paper lies essentially in the way the problem is stated.

There are many ways to define the state equation and the cost functional (observation) in the optimal control or controllability approach. We state a few open problems for future investigation, which cover various fields such as controllability theory, optimal control both from the theoretical and numerical viewpoints, the relationship between Fourier analysis and the control approach (Fourier aliasing effects and numerical approximations for the wave equation), shape optimization and so on.

The paper is organized as follows: in the next section, we present the physical problem and the model, together with the integral formulation (involving the spherical Radon transform). In Section 3, we present the inverse problem and usual variational techniques to address the ill-posedness. We propose an optimal control approach in section 4, present approximated and regularized models and derive general optimality conditions. In Section 5, we present simple numerical tests to compare the proposed approach to some classical methods.

The direct TAT model

A precise description of the physical model of TAT (or PAT) can be found in [30]. We outline here a simplified model. Let us denote by v(x, t) the fluid velocity at position x and time t, and make the following physical assumptions:

1. v(x, t) is small; 2. the fluid is non viscous and non turbulent; 3. there is no external force.

We also assume that the density ρ = ρ(x, τ ) and pressure p = p(x, τ ) undergo small variations: ρ = ρ 0 + δρ with |δρ| << ρ 0 and p = p 0 + δp with |δp| << p 0 .

This allows for linearizing the mass conservation and momentum conservation equations, which then read ∂ρ ∂τ = -ρ 0 ∇ • v and ρ 0 ∂v ∂τ = -∇p, respectively. In the latter equations and from now on, ρ and p stand for the variations δρ and δp.

Combining these equations yields

∂ 2 ρ ∂t 2 -∆p = 0. (2)
In order to link ρ and p with the absorbed electromagnetic radiation, one makes use of the thermal expansion equation, whose linearized form reads

∂ρ ∂t = 1 c 2 ∂p ∂t - β c p r. (3)
Here, β, c p and c denote respectively the thermal expansion coefficient, the specific heat capacity and the (adiabatic) speed of sound, while r = r(x, t) denotes the absorbed electromagnetic power. From ( 2) and (3), we easily obtain

1 c 2 ∂ 2 ∂t 2 -∆ p = β c p ∂r ∂t . (4)
Now, the absorbed electromagnetic power r(x, t) is related to the absorption coefficient ψ = ψ(x) by the equation r(x, τ ) = I(x, τ )ψ(x), in which I(x, τ ) is the radiation intensity. Finally, due to the high magnitude of the speed of light, one can assume that I(x, t) takes the separated form I(x) = J(x)j(t). Equation (4) then reads

∂ 2 ∂t 2 -c 2 ∆ p = βc 2 Jψ c p ∂j ∂t .
The function f (x) := βc 2 (x)J(x)ψ(x) c p is referred to as the energy deposition function. Finally, we can assume that the pressure increment, together with its time derivative, is zero at the initial time (right before illumination). We are then led to consider the following system:

(T AT )              ∂ 2 ∂t 2 (x, t) -c 2 (x)∆p(x, t) = f (x) ∂j ∂t (t), p(x, 0) = 0, ∂p ∂t (x, 0) = 0.
To simplify the problem, it is usually assumed that the sound speed c is constant (normalized to 1). An alternative formulation of the above linear direct problem is an integral formulation. Recall that the direct problem consists in determining p(x, t) for all (x, t) ∈ R 3 × R + from the knowledge of f and j. It is well known (see e.g. [19]) that the solution is given by

p(x, t) = dj dt ⊛ tRf (x, t). (5)
Here, the convolution operation ⊛ is defined by

g ⊛ h (x, t) = t 0 g(t -s)h(x, s) ds
and the operator R, referred to as the spherical Radon transform, is defined by

(Rf )(x, t) := 1 4π S 2 f (x + tω) dω.
In TAT, one is rather confronted to the inverse problem:

Recover the energy deposition function f (x) from measurements of p(x, t) for x over a surface S outside the illuminated fluid.

We now show that additional approximations give rise to an integral formulation of the latter inverse problem, under the assumption that the intensity profile j(t) is nearly a Dirac centered at the origin. Assuming that j has support in [0, T ] (see Figure 2 (left)), one may write, for every smooth function h(x, t),

d dt j ⊛ h (x, t) = dj dt ⊛ h (x, t). If j ≈ δ, then dj dt ⊛ h ≈ dh dt , so that dj dt ⊛ tRf (x, t) ≈ d dt (tRf ). Consequently, Equation (5)
can be approximated by

p(x, t) = d dt (tRf ). ( 6
)
We see that, after integrating and dividing by t, measurements of p(x, t) on S × R + give access to approximate values of (Rf )(x, t) on the same domain. However, this integral formulation is purely linear and cannot be genelarized to non-linear models. From the numerical point of view, the integral formulation leads to the so-called filtered back-projection method.

Variational techniques for the inverse problem

The integral formulation of the problem leads to a linear inverse problem of standard form, that is:

Recover f from an approximate knowledge g of Rf , where R is a compact linear operator (modelling the data acquisition process) from a norm space into another.

In TAT, R is the spherical Radon transform. As already mentionned, there are several methods to solve the inverse problem using exact or approximate formulas to estimate the inverse of R.

A classical approach to solve (ill posed) inverse problems is the least-square formulation using a regularization process (a penalization term for example) to get the well posedness. The general formulation stands

min L(g) -A(f ) + αH(f ) ,
where g is the data and f the function to recover. We describe A, L operators thereafter and H is a regularizing operator. The functional framework has to be made precise as well. In the sequel we focus on two approaches:

• The first one is the regularization by mollification which is inspired by the integral formulation : here A = R (the spherical Radon transform) and the operator L is a linear operator that provides a preprocessing transformation for the data g.

We call regularization by mollification a set of reconstruction methodologies deriving from the idea that the original ill-posed problem of reconstructing the unknown object f should be replaced by that of recovering a smooth version of it. In this approach, the target object is no longer f but φ * f , where φ is a convolution kernel and * denotes the convolution operation.

This idea has been developed in two different ways, independently. One of them, which bears the name of approximate inverse, was introduced by Louis and Maaß [26,27,28]. In the context of TAT, it was developed by Altmeier, Schuster and Scherzer [19]. The other way goes back to Lannes et al. [23], who introduced a particular variational regularization in the context of deconvolution and aperture synthesis. This approach was further analyzed and developed more recently by Alibaud et al.

[1] and Bonnefond et al. [11] (see also [START_REF] Alibaud | Analyse asymptotique de la régularisation par mollification, Actes du colloque: Mathématiques pour l'image[END_REF]).

• The second approach is an optimal control approach that derives from the PDE formulation. It seems to be more general than the previous techniques since the PDE can be non linear what the integral formulation does not allow. In this case, L is the identity operator and A associates the solution f := A(g) to equation (1) (with g instead of u 0 ) to the data g. The regularizing term H(f ) comes from a priori knowledge of the solution f .

We first briefly present the regularization by mollification

Approximate inverses

Given a continuous and nonnegative function φ with φ = 1, one defines the family

φ β (x, y) := 1 β n φ x -y β , β > 0,
called an approximation of unity. The mollified version of f , denoted f β , can expressed as the scalar product of f by φ β (•, y):

f β (y) := f (x)φ β (x, y) dx = f, φ β (•, y) .
Now, it is well-known that, in several senses,

f (x)φ β (x, •) dx → f as β ↓ 0.
Assuming that φ β ∈ ran R * , let ψ β be defined by

ψ β (•, y) = (R * ) -1 φ β (•, y). (7)
Then, one has:

f β (y) = f, φ β (•, y) = f, R * ψ β (•, y) = Rf, ψ β (•, y) .
We see that if one can compute explicitly ψ β , f β is readily obtained from the data g = Rf by taking the scalar product with the so-called reconstruction kernel ψ β . In the case where R * is not invertible, one may replace (R * ) -1 in equation ( 7) by the pseudo-inverse of R * (which corresponds to the least-square solution of the equation R * ψ(•, y) = φ β (•, y)). Clearly, the difficulty in this method lies in general in equation [START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF] which is an ill-posed problem as well. However, as mentioned above this nice duality trick has been successfully applied in the context of TAT by Haltmeier et al. in [19].

Regularization by mollification

We assume here that the original object f 0 has support in a bounded domain Ω 0 in R d , and that the mollified object φ * f 0 has most of its support in Ω ⊃ Ω 0 . Notice that, if φ is chosen to have compact support, the definition of Ω is straightforward. We also assume that R maps L 2 (Ω) into some infinite dimensional separable Hilbert space G, and that R is compact and injective.

The convolution kernel is regarded as a member of the family

φ β (x) := 1 β d φ x β , with φ ∈ L 1 (R d ) and R d φ(x) dx = 1. ( 8
)
The latter family is an approximation of unity, and the choice of β determines the target resolution.

Clearly, in this approach, the relevant regularization parameter is β.

Let C β denote the convolution by φ β . We now aim at approximating the new target object C β f = φ β * f 0 , which is performed by solving the following optimization problem:

(P β ) Minimize 1 2 Φ β g -Rf 2 G + α 2 (I -C β )f 2 L 2 s.t. f ∈ L 2 (Ω),
in which I denotes the identity and Φ β : G → G is itself a solution to:

(Q β ) Minimize 1 2 RC β -XR ⌊E 2 L(E,G) s.t. X ∈ L(G), X = 0 on (ran R ⌊E ) ⊥ .
Here, E is a subspace of L 2 (R d ), L(E, G) denotes as usual the space of continuous linear mappings from E to G and L(G) := L(G, G).

The choice of the functional to be minimized in (P β ) can be explained as follows. Any object f is the sum of its low frequency component C β f and of its high frequency component (I -C β )f The functional to be minimized acts of each component as independently as possible: the fit term introduces constraints on the low frequency component while the regularization term deals with the high frequency component. The complementarity of the filters φβ et 1 -φβ allows for a smooth transition between the experimental and regularization constraints.

We emphasize that, unlike most regularization techniques, the fit term requests adequacy to the regularized data Φ β g. The rationale for this is that, in some sense, Φ β g is optimally consistent with the new target object. Ideally, Φ β should be such that

Φ β R = RC β , (9)
which may not be feasible. Therefore, Φ β is chosen to minimize the discrepancy in Equation ( 9).

A control formulation

In this section, the pressure p will be denoted y, in order to fit the traditional notation in optimal control.

The state equation

We have seen in the previous section that the equation that describes the behavior of the system is (in R n with n = 3):

         ∂ 2 y ∂t 2 (t, x) -c 2 (x)∆y(t, x) = u o (x) ∂j ∂t (t), (t, x) ∈ [0, T ] × R n , y(0, x) = 0, x ∈ R n , ∂y ∂t (0, x) = 0, x ∈ R n . ( 10
)
where j is the electromagnetic pulse intensity profile and c the sound speed.

More generally we set u(t, x) := u o (x) ∂j ∂t (t). The function to be recovered u is regarded as a control function while the pressure y is regarded as the state function. Therefore Equation ( 10) is the so-called state equation.

Remark 1 By using Duhamel's principle [13], one may rewrite this state equation as

         ∂ 2 y ∂t 2 (t, x) -c 2 (x)∆y(t, x) = 0 in [0, T ] × R n , y(0, x) = u o (x), in R n , ∂y ∂t (0, x) = 0 in R n . (11)
So, in this context, the control function can be a distributed one or an initial one. However, this equivalence is only valid in the linear case and we focus on the distributed case.

Recall that, if V is a norm space endowed with • V , then for every p ∈ [1, ∞]

L p (0, T ; V ) := {y : [0, T ] × R n → R | y(t, •) ∈ V a.e. t ∈ (0, T ) and T 0 y(t, •) p V dt < +∞}.
and

C(0, T ; V ) := {y : [0, T ] × R n → R | ∀t ∈ (0, T ) y(t, •) ∈ V and t → y(t, •) V is continuous }.
The original problem domain is R n . This is consistent with the fact that one usually uses classical Fourier analysis (Radon transform, spherical means,...) since the uncertainty principle do not allow compact supports both in the time and frequency domains. The control viewpoint requires the definition of boundary conditions (to infinity if the domain is R n ) to use duality methods and Green formulae. By contrast, the advantage of this approach is that we may choose an open bounded subset B as the domain and use finite differences or finite elements to compute the solution. A challenging application is to set a model which allows for using both spectral (Fourier) analysis and finite difference methods. We recall a regularity and existence abstract result which gives indication about the functional framework to be used (with appropriate regularity assumptions on c and c ≥ c 0 > 0 ):

Theorem 4.1 ([25] Vol. 1 p 286) If u ∈ L 2 ([0, T ] × R n ), then (10) has a unique solution in { y ∈ L 2 (0, T, H 1 (R n )) , y t ∈ L 2 ([0, T ] × R n ) . Moreover y ∈ C(0, T ; H 1 (R n )) and y t ∈ C(0, T ; L 2 (R n )) ,
and the solution y continuously depends on u.

In the sequel, we shall avoid using the whole space R n . Therefore, we consider that infinity means far enough and that the problem domain can be an open ball B ⊂ R n with large radius. Indeed the pressure y quickly vanishes and we may consider that it is equal to 0 on the boundary ∂B of the ball. The ball B contains Ω as well as the sensors location ω ε (we shall justify the notation later on). If we impose Dirichlet conditions on the boundary, we get:

             ∂ 2 y ∂t 2 (t, x) -c 2 (x)∆y(t, x) = u(t, x), in ]0, T [×B := Q, y(0, x) = 0, in B ∂y ∂t (0, x) = 0 in B y = 0 on S (12)
where S :=]0, T [×∂B is the lateral boundary of Q.

Remark 2 We may of course consider Neumann boundary conditions or mixed conditions: this requires additional investigation on the physics. We may also consider R n as the domain using appropriate decay properties at infinity instead of the ball: the techniques are quite similar.

We have a generic classical result (see [14] p. 632 for example) which can be sharpened of course. To simplify the presentation, we now assume from now that the sound speed c is equal to 1.

Theorem 4.2 Assume that u ∈ L 2 (0, T ; H -1 (B)). Then, Equation (12) has a unique solution y

[u] ∈ L 2 (0, T ; H 1 o (B)) ∩ C 0 ([0, T ]; L 2 (B)
). In the sequel we denote

W := L 2 (0, T ; H 1 o (B)).
Open problem 1 We have set u(t, x) := u o (x) ∂j ∂t (t) where j is a smooth approximation of the Dirac measure. This gives the natural regularity of the control function. However, an interesting challenge is to perform the following analysis using less smooth control functions.

The control problem

Following Scherzer [19] we assume that u has a compact support contained in [0, T ] × Ω, where Ω is a bounded open subset of R n containing the support of the object to be recovered. The pressure y is measured (observed) on Σ = [0, T ] × Γ where Γ ⊂ ∂Ω is contained (most of time strictly) in ∂Ω. The inverse problem which consists in recovering u from measurements of y on Σ can be formulated as the following minimization problem:

Minimize Σ (y(u) -y d ) 2 dt dγ s.t. u ∈ U,
where y d is the measured pressure and U is the admissible control function space. Many difficulties then arise:

• The problem is known to be ill-posed: we cannot even ensure the existence of a solution.

One must add a regularization term for u. A standard way to perform the regularization consists in adding to the functional a Tychonov term:

α Q u 2 dx dt with α > 0,
or rather term inspired by the regularization by mollification method

α Q (B β u) 2 dx dt where α > 0, B β (u) = u -φ β * x u, (13)
and φ β is defined as in subsection 3.2. Here * x denotes the spatial convolution (with respect to x). More precisely

Q (B β u) 2 dx dt = T 0 B (u -φ β * x u) 2 (x, t)dx dt .
Note that B β = I gives the classical Tychonov term.

Open problem 2 A crucial point is to investigate the physical model and the qualitative properties of the control in order to choose an appropriate regularization term: we may choose the solution with minimal energy, but it may not be the best choice. We choose the above regularization by mollification term for the moment. This choice has to be considered precisely in the future.

Open problem 3 A natural question is that of describing the asymptotic behavior of the solution as the regularization parameter α goes to 0 (this is well-known for Tychonov regularization process for example) and perform sensitivity analysis with respect to β.

• In order to define the observation term in the cost functional, the trace of y[u] on Σ has to belong to L 2 (Σ). This depends on the regularity of u. However, in the case of point detectors, the observation term would take the form

N i=1 [0,T ] (y[u] -y d ) 2 (t, x i ) dt = N i=1 [0,T ]×R n (y[u] -y d ) 2 (t, x) dt dδ xi (x),
in which dδ xi denotes the Dirac measure at x i . Notice that one may relax the above functional and write instead

N i=1 [0,T ]×R n (y[u] -y d ) 2 (t, x) dt dµ i (x),
in which dµ i is a relaxed version of dδ xi , e.g. a gaussian measure centered at x i . Clearly, the width of each µ i should be less than, say, one half of the expected wavelength of the pressure wave, so that the requested fit is physically acceptable.

Open problem 4 Study the asymptotic behavior of the solution as the fit measures µ i approach the Dirac measures δ xi .

• Even if the cost functional makes sense, using duality methods and integrating by parts is not possible since Γ is not a closed surface (roughly speaking, there is no inside and outside).

To overcome the last difficulty we regard Γ as the "limit" of a family of open subsets ω ε : lim

ε→0 ω ε = Γ.
This convergence must be clarified of course (using tools from shape optimization, as Γ-convergence for example). A simple way to perform this approximation is to set

ω ε = x∈Γ B(x, ε),
where B(x, ε) is the open ball of radius ε centered at x (see Figure 3). In the case where there is a finite number of sensors located at

x i , 1 ≤ i ≤ N, x i ∈ Γ one may choose ω ε = N i=1 B(x i , ε) . ( 14
)
So, we replace a boundary observation by a distributed observation. We set Ω ε = Ω ∪ ω ε and Γ ε = ∂ω ε .

Open problem 5

• More generally Ω ε is defined as a perturbation of Ω: there are many questions related to the asymptotic behavior of the ε-solution as ε → 0. One may use tools as the Γ-convergence for example.

• Another interesting question is to look for an optimal ω ε with respect to a criterion given by the experimental setup for example. Finally, the regularized optimal control problem reads:

(P ε )      min J ε (y, u) := 1 2 [0,T ]×ωε (y -y d ) 2 (t, x) dt dx + α 2 [0,T ]×Ω (B β u) 2 (t, x) dt dx y = y[u] solution to (10), u ∈ U
where α > 0, y d ∈ L 2 (Q) and U is defined for example as

U := {u ∈ L 2 (]0, T [×B) | u(t, x) = 0 a.e. (t, x) ∈]0, T [×(B\Ω) } .
Remark 3 On can choose

U := {u ∈ L 2 (]0, T [×B) |u(t, x) = u o (x) ∂j ∂t (t) a.e. (t, x) ∈]0, T [×Ω, supp(u o ) ⊂ Ω } .
In the case where the control is an initial one the regularization term could be

α 2 B (B β u o ) 2 (x) dx = α 2 Ω (B β u o ) 2 (x) dx .
Theorem 4.3 For every ε > 0 problem (P ε ) has a unique solution (y ε , u ε )

Proof.-The proof is standard.

Open problem 6 One has to precisely study the asymptotic behavior of (y ε , u ε ) with respect to ε and if possible give error estimates.

Optimality conditions

Once the optimal control problem has been formulated it is easy to use the classical machinery (see [24]) introducing the adjoint state p ε :

                 ∂ 2 p ε ∂t 2 -∆p ε = y ε -y d in [0, T ] × ω ε 0 elsewhere p ε (T, x) = 0 in B, ∂p ε ∂t (T, x) = 0, in B, p ε = 0, on S
We get the following optimality system Theorem 4.4 For every ε > 0 a necessary and sufficient condition for (y ε , u ε ) ∈ W × U to be the solution to

(P ε ) is                  ∂ 2 y ∂t 2 -∆y = u, in Q, y(0, x) = 0, in B, ∂y ∂t (0, x) = 0, in B, y = 0, on S (15a)                  ∂ 2 p ε ∂t 2 -∆p ε = y ε -y d in [0, T ] × ω ε 0 elsewhere p ε (T, x) = 0 in B, ∂p ε ∂t (T, x) = 0, in B, p ε = 0, on S (15b) ∀u ∈ U, Q (p ε + αB * β B β u ε )(u -u ε ) dt dx ≥ 0 , (15c)
where B * β denotes the adjoint operator of B β . Note that we have embedded the TAT model in a more general class of problems where the control function is no longer "decoupled". In the (TAT) case where u(t, x) = f 0 (x)f (t) we may use the very specific structure of the control function. A similar source identification problem has already been studied for example in [33] in a different framework, where the measured data was the normal derivative.

The controllability point of view

Another way to consider the problem, which seems to be the most natural one is a controllability approach. Roughly speaking, we look for a control function u with support in ]0, T [×Ω such that y[u] is given on the set Σ =]0, T [×Γ. It is clear that the solution is not unique if there is some. However the existence is not obvious both is the cases of exact and approximate controllability.

There is a huge literature on controllability/observability theory and we cannot mention all the papers. However, Bardos, Lebeau and Rauch have shown in a famous paper [9] that for the observation or control of solutions of second-order hyperbolic equation a necessary and sufficient condition for controllability is that the region of control meet every ray of geometric optics that has, at worst, transverse reflection at the boundary. More precisely, for multidimensional problems, the region of control must meet each ray in a non diffractive point. This condition is not ensured in our case since the observability (stabilization) region ω ε is not necessarily connected to the control region Ω.

Open problem 7 What could be a good choice of ω ε or Γ to give some existence results for the controllability problem: Find u with compact support in Ω such that y[u] = y d on Γ or ω ε ? Which model would be appropriate for using controllability techniques for the wave equation (HUM method, Carleman estimates) ?

Numerical experimentation

We now present some numerical simulations to compute a simple TAT model with a closed observation surface (a sphere) and a constant sound speed c ≡ 1. More precisely we want to recover the source u 0 which drives the following equation

         ∂ 2 y ∂t 2 -∆y (t, x) = 0, (t, x) ∈ [0, T ] × Ω, y(t, x) = 0, (t, x) ∈ [0, T ] × ∂Ω, y(0, x) = u 0 , ∂y ∂t (0, x) = 0, x ∈ Ω (16)
from measurements y obs on the boundary of Ω.

Here Ω is the 0 -centered ball of radius 1. Thanks to symmetry we reduce the problem to a 2D problem.

We investigate three methods

• The basic filtered back-projection method;

• The Time Reversal method;

• The resolution of system ( 15) with a conjugate-gradient method.

These methods have been tested using Shepp-logan phantom. To perform a rough comparison, we do not use the recent improvement of these methods (see for example [19,[START_REF] Ammari | Photo-acoustic imaging for attenuating acoustic media[END_REF][START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF][START_REF] Qian | An efficient Neumann-series based algorithm for Thermoacoustic and Photoacoustic Tomography with variable sound speed (with )[END_REF]). In that spirit we do not compare precisely the convergence speeds and the CPU run time, since we have not used optimized codes and methods. Note that we have chosen an experimental setting where the comparison of these techniques was relevant (constant sound speed and closed observation surface).

5.1

The filtered back-projection

The filtered back-projection (FBP) method remains one of the most commonly used techniques to invert integral operators such as the classical Radon transform or the spherical Radon transform R, for it provides an exact inversion formula in some specific acquisition geometries. Concerning the thermo-acoustic tomography, since these formulae allow the inversion of the spherical Radon transform, they require the strongest approximation we have discussed: a constant speed of sound, a perfect Dirac pulse and a complete data set (most of the time, a whole sphere). The first known FBP formula was derived in [16] for odd dimensions. Later on, in [19], the authors used the approximate inverse concept to regularize this formula. Here, we implement an algorithm proposed by Finch in [15] for the two-dimensional case, which is based on the closed formula:

f (x) = 1 2π S 1 2 0 ∂r(r∂rR 2d )(ω, r) log r 2 -|x -ω| 2 drds(ω),
for smooth functions f supported in the unit ball. The derivations are approximated with symmetric finite differences and the numerical integrations are performed by means of linear interpolations and of the trapezoidal rule. This algorithm is very efficient and quite stable when one deals with a sufficiently large amount of detectors (Figures 5 and7). Nevertheless, the implementation becomes unstable when the angular sampling of the detectors around the unit circle doesn't allow the trapezoidal rule to be accurate (Figure 8).

Time reversal

The idea of this method is the following (see [20]): thanks to the Huygens' principle (in odd dimension), the wave leaves the unit disc in a finite time, such that there exists a time T such that the pressure vanishes for any t ≥ T . Consequently, as the wave standard equation is reversible, we may solve the backward problem starting from time T with boundary conditions equal to the measured data to recover the initial condition. This method has been adapted to the case where space dimension is even (see [START_REF] Qian | An efficient Neumann-series based algorithm for Thermoacoustic and Photoacoustic Tomography with variable sound speed (with )[END_REF]) and in attenuating acoustic media [START_REF] Ammari | Photo-acoustic imaging for attenuating acoustic media[END_REF][START_REF] Ammari | Time reversal in attenuating acoustic media[END_REF]. The backward problem was computed with a classical finite differences scheme.

The control approach

There are many classical tools to make this algorithm converge (splitting, relaxation). The more delicate part is the computation of Π U since the support constraint is a pointwise constraint and we perform a L 2 -projection. An alternative is to keep the constraint as a control constraint and deal with the associate Lagrange multiplier. This will be addressed in future work. Nevertheless, in a linear framework, Problem (P ε ) reduces to:

(P ε )      min 1 2 W u -y d 2 L 2 ([0,T ]×ωε) + α 2 B β u 2 L 2 ([0,T ]×Ω) u ∈ L 2 ([0, T ] × Ω),
where the linear operator W is defined by W u := y[u], so that the optimality condition (15c) reads:

(W * W + αB * β B β )u ε = W * y d . (17)
Here, B β is defined by ( 13) and the kernel φ β is a gaussian function with standard deviation β ≈ 3.10 -3 . We have set α = 0.4. Finally, ω ε is given by ( 14), where ε is set to a single pixel size.

Since the self-adjoint operator W * W + αB * β B β is positive definite, we chose to solve Equation ( 17) by means of the conjugate gradient algorithm. It is straightforward that W * W u ε is nothing but the solution p ε of Equation ( 15b) with y d replaced by 0. Moreover, W * y d is obtained as the solution of Equation (15b) with y ε -y d replaced by y d . The forward and backward problems (15a) and (15b) are solved by means of a leapfrog discretization scheme on a staggered grid.

In order to avoid handling large grids (due to the size of B), we use an appropriate PML (Perfectly Matched Layer) technique (see [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]).

Numerical tests

Each method is tested on the Shepp-Logan phantom: We investigate the cases where • data are not corrupted;

• a white gaussian noise (SNR=0.15) is added to the simulated data;

• the number of detectors is either 50 or 360 (with a uniform angular sampling on the whole circle).

Influence of the mollification

In order to show the relevance of the chosen regularization process, we present thereafter the solutions obtained with the control approach for α = 0 and α = 0.4. If α = 0, the existence of a solution of Problem (P ε ) is not ensured. However, it has been possible to compute a minimizing sequence. The solutions presented in Figure 9 are obtained with 20 iterations of a conjugate gradient algorithm.

Conclusion

We hope we have clearly presented many open questions about a challenging imaging process. We can apply these remarks to photo-acoustic tomography (PAT) as well: the equations are the same (wave equation related to acoustic part) but the frequency scale is different. There are many others questions that arise as for example the optimal location of sensors (shape-optimization problem), variable speed sound and/or non linear equations, filtering effect on u and so on.

The numerical tests we have presented are to be performed with more accuracy. Indeed, the sensitivity to parameters (α, β) has to be deeply investigated. Moreover, it is known that for time semi-discrete systems, due to high frequency spurious components, the exponential decay property may be lost as the time step tends to zero (see Zuazua [34,35]) so that numerical schemes are highly unstable. In [17], Ervedoza and Zuazua use a decoupling argument of low and high frequencies, the low frequency observability property for time semi-discrete approximations of conservative linear systems and the dissipativity of the numerical viscosity on the high frequency components. All these questions will be discussed in a future work.

Figure 2 :

 2 Figure 2: Physical principle of Thermo-acoustic Tomography: intensity time profile of the EM pulse (left) and scheme showing the pressure wave generation and propagation (right).

Figure 3 :

 3 Figure 3: Observation surface relaxation

Figure 4 :

 4 Figure 4: The Shepp-Logan phantom (256×256 pixels).

Figure 5 :

 5 Figure 5: Noiseless data, 360 detectors. Left: filtered back-projection solution. Center: time reversal solution. Right: control approach solution.

Figure 6 :

 6 Figure 6: Noiseless data, 50 detectors. Left: filtered back-projection solution. Center: time reversal solution. Right: control approach solution.

Figure 7 :

 7 Figure 7: Noisy data, 360 detectors. Left: filtered back-projection solution. Center: time reversal solution. Right: control approach solution.

Figure 8 :

 8 Figure 8: Noisy data, 50 detectors. Left: filtered back-projection solution. Center: time reversal solution. Right: control approach solution.

Figure 9 :

 9 Figure 9: Left (resp. Right): control approach solution without (resp. with) regularizationnoiseless data, 50 detectors.

Figure 10 :

 10 Figure10: Convergence rates for α = 0 (continuous line) and α = 0.4 (dashed line)-The relative error is computed as u c -u orig 2 / u orig 2 where u c is the computed solution and u orig the (original) object to recover.


Hakiri Akram
email: hakiri@laas.fr

Berthou Pascal
email: berthou@laas.fr

Gayraud Thierry
email: gayraud@laas.fr

Controlled Stochastic Petri Net Model for End-to-End Network QoS Provisioning in Middleware-based Multimedia and Real-Time Systems

Keywords: Controlled SPN, Quality of Service, Multimedia System, Heterogeneous Networks

End-to-end quality of service (QoS) is central to the objectives of the today's networks requirements of middleware based distributed real-time and embedded (DRE) systems. Any middleware based QoS system should be totally oriented to this goal, and in the scope of this purpose several mechanisms, components and approaches were, are being and will be developed in order to achieve it. In this paper, we show how controlled behavior of such QoS-aware systems can be developed based on stochastic Petri Nets. Afterwards, We show how to obtain, using such an interpreted formal model, powerful numerical analysis for the management of the network QoS.

INTRODUCTION

Distributed real-time and embedded (DRE) middlewarebased systems refer to a distributed platform of interfaces and services that reside 'between' the application and the operating system and aim to facilitate the development, deployment and management of distributed applications. They often comprise multiple end-to-end application flows that may require various QoS properties affecting CPU, memory, and network resources [START_REF] Schmidt | R&D Advances in Middleware for Distributed, Real-time and Embedded Systems[END_REF]. Examples include, among others, unmanned air vehicle, experimental aircraft for disaster recovery, data-centers, space mission systems, video surveillance, real-time and on demand video transmission, homeland security, on line stock trading, and weather monitoring. In most of these systems, the right answer delivered too late becomes the wrong answer, i.e achieving the end-to-end QoS is essential. Most challenging DRE systems will operate in large scale, distributed and heterogeneous environment that take data from large number of sensors, multimedia sources, realtime streaming systems, and remote control systems [START_REF] Schmidt | Middleware R&D Challenges for Distributed Real-time and Embedded Systems[END_REF].

To address the run-time QoS demands, research efforts focus on the QoS-enabled components and emerging technologies including middleware; the basic mechanisms for end-toend QoS provision are located at the data plane (e.g. packet scheduling) and control plane (e.g., admission control and resources reservation). Signalling is carried out at various levels: end-to-end, at application level, through the use of SIP [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF]; hop-by-hop, at inter-domain level, complemented by routing and/or traffic engineering; locally, at intra-domain level, for resource management and provision. Likewise, the middleware is a distributed software layer which shields the appli-cation from the complexity and heterogeneity of the underlying distributed environment. The middleware was designed to help the development and management of distributed systems. Examples of middleware include the High Level Architecture (HLA), Java Message Service (JMS), CORBA Event Notification, OMG Data Distribution Service (DDS), and Grid Computing.

Recently, significant research has been made towards making software as industrial reusable components, evolving access to information, computer and network resources over middleware-based distributed systems. In order to assess the adequacy of our middleware model-based approach to support effective QoS policies, our approach consists of the following:(1) identification of the QoS policies and the analysis of the QoS constraints, (2) proposal of Stochastic Petri Net model-based approach for middleware-based Network QoS provisioning. To the best of our knowledge there is not any other published paper dealing with this specific topic of QoS provisioning including together the application, the middleware and the network resource allocation. In this paper, we show how controlled behavior of such QoS-aware systems can be developed on the control plane and how it invokes the middleware services using some kind of Petri Net formalism.

The remainder of this paper is organized as follows: after this introduction, Section 2 briefly describes the theory of Generalized Stochastic Petri Nets (GSPN). Section 3 uses case studies to motivate common requirements of middleware-based multimedia and DRE systems associated with the network QoS provisioning. Section 4 analysis the QoS policies that should be managed for such system and addresses the key challenges associated with meeting these requirements. GSPN middleware-based model used in our approach is described in section 5, as well as the results from simulation. Section 6 compares our work with related researches; finally, some concluding remarks and future works are presented in Section 7.

GENERALIZED STOCHASTIC PETRI NETS

In most QoS research area, model driven approach [START_REF] Sandeep | Generators for Synthesis of QoS Adaptation in Distributed Real-Time Embedded Systems[END_REF], state machine [START_REF] Cavalli | New approaches for passive testing using an Extended Finite State Machine specification[END_REF], dataflow analysis [START_REF] Sandeep | Constraint-based design-space exploration and model synthesis[END_REF] are used to explore possible solution for QoS assurance between components; they focus on deterministic process analysis and how QoS parameters can interoperate between each other. However, these approaches do not suffice for QoS analysis, because when analyzing the QoS provisioning, additional information's are needed: when an application is sharing differ-ent flows, each flow may have specific QoS parameters. For example, fire sensor should not have the same data priority as camera surveillance; data should be sent in reliable transport service with high priority. Hence, providing some kind of urgency, probability of events occurring and precedence between camera flow and fire sensor flow is very important. All these characteristics show the difficulty for the previous approaches to ensure the QoS requirements of DRE systems.

Generalized Stochastic Petri Nets (GSPN) [START_REF] Diaz | Petri Nets Fundamental Models, Verification and Applications[END_REF] are wellknown interpreted extensions of autonomous Petri Net (PN) [START_REF] Diaz | Petri Nets Fundamental Models, Verification and Applications[END_REF]. Petri Nets are similar to those approaches but have additional benefits in modeling. Since graphs in Petri Nets allow the representation of possible assembled DRE systems, it consider the performance aspect in the design of distributed and complex concurrent systems. Probabilistically time delay is introduced in transitions to allow performance evaluation and resolving scheduling problems in DRE. In addition to the PN's capability of structural and functional validation, GSPN provides numerical solution for efficient algorithms with Markov Chain (MC) [START_REF] Diaz | Petri Nets Fundamental Models, Verification and Applications[END_REF] to get the exact performance parameters such as throughput and delay. M 0 :P -→ N + , is the initial marking: a function that assigns a non negative integer to each place Π:P-→N, is the priority function which associates lowest priority (0) with timed transitions and higher priorities (≥ 1) with immediate transitions:

Π(t) = 0, if t is timed ≥ 1, if t is immediate (1) W:T-→ R

MOTIVATING MIDDLEWARE-BASED NET-WORK QOS PROVISIONING

Before detailing the GSPN model, a brief example is given to illustrate how and why the QoS GSPN-based model allows the design, specification and evaluation of QoS requirements of DRE systems.

Challenge 1: Alleviating the ambiguity in QoS Specification

Context. Network architecture having more complexity and quality requirements will need more precise specification for better control and design [START_REF] Schmidt | R&D Advances in Middleware for Distributed, Real-time and Embedded Systems[END_REF]. Cooperative enterprizes often use IP network over high speed network. Traffic can be grouped into several classes, including email, video conference, fire/smoke sensors data, camera surveillance, imagery traffic, etc. To provide transmitting the right data at the right place in the right time, each flow must specify its requirements to the underlying transport network service. As pointed out in [START_REF] Liu | QoS-UniFrame: a Petri net-based modeling approach to assure QoS requirements of distributed real-time and embedded systems[END_REF], DRE systems can be specified, modeled, analyzed, and evaluated with the aid of the structural and behavioral analysis of Petri Nets formalism. Problem. When constructing DRE systems, developers need to make decisions about the cooperative aspect of different DRE components. Each decision may induce diverse execution orders, execution time, and events. Therefore, there are huge number of possible of cases generated bases on different decisions [START_REF] Liu | QoS-UniFrame: a Petri net-based modeling approach to assure QoS requirements of distributed real-time and embedded systems[END_REF]. To reduce the complexity of exploring all the possible cases by evaluating their QoS requirements, the evaluation should expresses the specific time before/after the events affect the QoS parameters, the kind of parameters that should be evaluated, and how to carry out the relationship between them. Hence, there are no mature approach, principle or solution to enable large scale DRE systems to be QoS predictable, repeatable, validated, and enhanced. Solution Approach. In our approach, we show that the modeling begins with a representation of adequate mechanisms for each application, and is, after each step, followed by the validation. Then, different behavior can be represented and their global behavior evaluated. The specification will be considered correct, when, after validation the general properties of the model are correct.

Challenge 2. Alleviating the complexity of resource

reservation and deployment Context. The manner on which the reservation strategies in presence of middleware depends on the number and the type of calls received from different hosts [START_REF] Kinnebrew | Intelligent Resource Management and Dynamic Adaptation in a Distributed Real-time and Embedded Sensor Web System[END_REF]. For example, Homeland security sensors and video surveillance camera are more sensible for delays, bandwidth. Basically, each kind of traffic should specify its QoS requirements [START_REF] Roy | Toward Effective Multi-capacity Resource Allocation in Distributed Real-time and Embedded Systems[END_REF] to the network during the call admission control phase with the aid of some kind of signaling protocol. Problem. This area represents the next great wave of evolution for middleware. Despite the ease connectivity provided by middleware, however, the DRE traffic has always commanded to respect some kind of Service Level Agreement (SLA) [START_REF] Nanda | Providing end-to-end guaranteed quality of service over the Internet: a survey on bandwidth broker architecture for differentiated services network[END_REF] to avoid the situations where penalties may be paid. The premises underlying the push towards end-to-end resources reservation mediated by middleware are that different level of service are possible under different conditions, costs, and constraints. A tradeoff between the levels of service is used to respect the constraints. Solution Approach. The service level involved between the middleware and the control plane is fulfilled, when, after validation of its properties, the following problems are solved: accessibility (how to grant and control the access to the service), controllability (how to control service instances), and resources allocation (how to locate and define the admission control strategy and the resources reservation). The model developed in this approach provides one way to translate the QoS requirements from one level to another (application, middleware, transport, and network). For example, as traffic from video application is expected to be a substantial portion of the traffic carried by emerging wired and wireless networks, video source should use variable bit rate at the encoder to manage the errors, packet lost (especially in wireless networks) and delay (delay is annoying to the viewer) at the application level and notify the middleware about these changes. Middleware should, at its turn, notify the control plane to adapt the resource reservation strategy to changes.

Challenge 3. Reducing the complexity of the multi-

domain end-to-end QoS guarantee Context. To support end-to-end QoS for DRE applications, the complete chain of involved systems and the network should support it. Ideally, applications make use of end-toend signaling to signal its end-to-end QoS requirements and agreed QoS between each other [START_REF] Manghwani | End-to-End Quality of Service Management for Distributed Real-Time Embedded Applications[END_REF]. For instance, Video on demand applications should include traffic description (delay, bandwidth, jitter, packet lost,...) describing which kind of treatments by the network element. Since the Differentiated Service (DiffServ) and the Integrated Service (IntServ) [START_REF] Nanda | Providing end-to-end guaranteed quality of service over the Internet: a survey on bandwidth broker architecture for differentiated services network[END_REF] are the most used architecture for end-to-end QoS provisioning, flows having different priorities, delays, and bandwidth requirements, have to obtain the required network resources fulfilling these requirements. The network backbone (access and core network) should be able to allow the transmission of each flow without affecting other flows, and permit the call admission control of incoming demands. Problem. Despite the advances of the currently available network architecture and models for resources handling (Diff-Serv, IntServ), however, it is not sufficient to deploy them within routers to guarantee QoS between two entities in the Internet. End-to-End QoS system should be able to coordinate the routers in the network according to the users' demands. For example, communication provided the nowadays Internet infrastructure has, by nature, multi-domains and multitechnologies architecture. Data packet travel across many domains, including diverse Internet Service Providers (ISP) in different countries, each has its own deployed network backbone, what makes hard to provide end-to-end QoS guarantee [START_REF] Schmidt | R&D Advances in Middleware for Distributed, Real-time and Embedded Systems[END_REF]. Solution Approach. The concept of policy control developed in our approach aims to provide QoS SIP-based [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF] Petri Nets model to allow application expressing their QoS requirements using text-based exchanged messages. Like shown in Figure 1, the SIP/SDP [START_REF] Camarillo | Integration of Resource Management and Session Initiation Protocol (SIP)[END_REF] model coordinates the application signaling, the resources reservation and provides to the Bandwidth Broker (BB) [START_REF] Nanda | Providing end-to-end guaranteed quality of service over the Internet: a survey on bandwidth broker architecture for differentiated services network[END_REF] the facilities to coordinate the routers in the network, that applies to both IntServ and DiffServ infrastructure. Under the proposed Petri Net model we consider an admission control policy fulfilling both per-flow and Classbased end-to-end guaranteed service to allow QoS-enabled network. Stochastic priorities assigned to transitions in GSPN model allow the better scheduling of different incoming calls to satisfy the QoS requirements of each flow. DRE applications can select the necessary firing time depending on the probabilistic value of the timed transition.

GLOBAL MODEL DESCRIPTION

It is important to define a common modeling system, using a level of abstraction that will permit the quantitative analysis of non-trivial system, then we will detail the architecture of the proposed model including the application, the middleware and the network. Figure 1 describes the global architecture of architecture adaptively control the video transmission captured from cameras, then send them to middleware-based distributed process to be distributed over QoS-Aware network to one or more remote receivers (end-systems) including displayers and image processing software. Application/Middleware Interface: it aims providing users the ability to specify high-level management policies: the registration mechanisms and controlling how flows are exchanged and delivered. Middleware: when considered in the GSPN Model, it is considered as a Message-Oriented Middleware (MOM). Moreover, the middleware is composed of Message Manager and Queue Manager to process incoming flow. This abstraction facilitates the modeling of the whole system. Middleware/Network Interface: it uses a QoS allocator to notify the Control Plane about the middleware data and the messages that will be exchanged during the reserved section. Service Plane: Application will use SIP to reserve the communication path and channel signaling. For example, Multimedia application use it to the renegotiate the codec, but it uses classical transport protocol for data transmission. Control Plane: This common control plane promises to simplify network operation and management by automating endto-end provisioning of connections, managing network resources and providing the level of QoS. It includes two mechanisms: the admission control unit and the resource reservation. The transport plane is an abstraction of the underlying interconnection technology.

GSPN MODEL FOR QOS PROVISIONING

The architecture of the GSPN model for QoS provisioning is presented in Figure 2. We assume, for clarification seeks, that we have 2 types of flow, but the model can be easily extended as well as for great number of streams.

Modeling

The main components are the Application (App), the Service Plane (SIP), The Middleware Interface (MI), the middleware (Midd), and the Control Plane (CP):

The application: is composed of the sender application (AppS) and the receiver application (AppR), and it is responsible for sending and receiving packets to and from the middleware, respectively. Those packets represented by tokens in the places Sender and receiver to show that application is ready to send/receive data. When modeling the packet lost of multimedia stream, the permissible delay between encoding and playback may be too low to allow retransmission or channel coding. A practical and widely used approach to this problem is considering that the codec can process the Packet-Loss Concealment (PLC) [5] schemes and masks the loss of a packet by using information from the last good packet. Hence, in the sender side, the codec self correction mechanism is represented in transition T1 by the stochastic arrival process with 'pdf' distribution "λ codec ". Since most modern decoders are capable of performing error concealment to mitigate losses, feedback can be used to notify the encoder of the losses, in particular applications using Session Description Protocol (SDP) need to agree on what voice codec to use for transmission and reception, the place "codec buffer" represent the way to provide a loss-aware rate-distortion scheme that can be used by both SIP and SDP protocols [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF]. The Service Plane (SIP): In this model we consider that we have two User agents at both the sender side and the receiver side described by the places 'SIP UA1' and 'SIP UA2', respectively. Each Agent is capable of generating requests or sending a response to a request. It is the user agents that react to one another. When a dialog is created, it is a user agent communicating with another user agent. The challenge 1 can be to agree the session end verified, when, looking at Table 1 and Table 2, the semantics of the SIP bloc using PN formalism, describing the behavioral and structural proprieties of this model is bounded, deadlock free and liveliness. We check whether or not a UA1 will wait infinitely. The net does not fulfill this property: starting from the initial marking, We can reach this marking at any time of the simulation (The reachability graph is not represented here, but it can be found intuitively). The model is well suitable for end-2-end signaling between the end applications. More, the effect of the firing sequence is determined by the incidence matrix and the vector of transition sequence, and has an effect on the behavior of the model. Since the number of marked places at the same time is equal to 1 at both sides of the SIP bloc, the protocol is conservative and only one operation is permitted at any given time.

Middleware Interface: in addition to the Concealment scheme provided by the coded, the middleware provides, via place P3, a notification interface to adapt the application bit rate to allow better user quality of experience. Further, when a token is present in place P3, this means that the network states is not able to transmit all the data incoming from the application. This case can occur in Internet network, where the number of client being admitted is greater. Hence, the quality of Reservation finished and end-to-end path is no yet required, so release resources by UA1 T20

Reservation finished and end-to-end path is no yet required, so release resources by UA1 service can be decreased at the application to avoid congestion.

The Middleware: since the signaling path is established by the SIP protocol, the application starts sending the data to the middleware. Place P0 represents the input buffer of the middleware used to marshal data. The place P1 describes the main internal buffer middleware and it takes a number of tokens (N) representing the maximum number of messages that this queue can process at a given time. In many practical cases, the middleware retrieves the data from the upper layer and processes them into its proper queue then gives them to the underlying network layer. Therefore, transition T4, which has rate α 2 , represents messages sent to the network via place P13 which represents the interface between the middleware and the network. Since we considered alleviating complexity of the service deployment, this hierarchical approach gives answer to challenge 2. We should note that at the receiver side the same abstraction of the middleware (main buffer, queue) are represented. Since this bloc acts just as a receiver the notification interface is not provided for simplicity seeks'. Therefore place P01 represent the incoming buffer of the receiver side middleware and places P11 and P21 show its main queue. Messages incoming from the network and traverse the middleware and, via transition T41, reaches the receiver application with distribution rate α 3 . This pdf function models the delay taken from the middleware to process the data at the upper layer, and this rate depends on the type of the stream (multimedia, real-time,...) and the underlying operating system (real-time, not real-time).

The Control Plane: the management at this plane is mainly implemented by two components: the resource manager (called here Call Admission Control ) and the resource allocator. Figure 2 shows the control plane in our model: the Call Admission control is located in bloc CAC, and the resource allocator is given in bloc RA. Table 3 and Table 4 illustrate the specifications and the semantics of both the places and the transition of the control plane: at the beginning of the simulation, available resources for the stream processing, presented by tokens in place P14 representing the length of the CAC queue. When a token is present in places P4, PC5 and/or PC6 with timeout mean duration 1 λ 1 and 1 λ 2

(we assume a pdf distribution), a call stream is coming to the CAC bloc and a stream will check for the availability of its resources requirements. In fact, the middleware component is responsible, via transition T4, for the specification of the QoS requirements of each flow incoming from the application. In a real case, An Service Level Specification (SLS) [START_REF] Nanda | Providing end-to-end guaranteed quality of service over the Internet: a survey on bandwidth broker architecture for differentiated services network[END_REF] contains exclusively the technical details specified by an SLA. It is essentially a translation of a SLA into the appropriate information necessary to configure network devices. Once a domain has agreed to honor the conditions set out in a SLS, it must be responsible for giving the guaranteed service to traffic specified in the SLS for the duration of the agreement. Hence, when a token arrives at place P4 by firing the timed transition T8 with timeout duration equal to 1 λ 3 , and then a call is incoming from a middleware with 'pdf' distribution and a Resource Allocation Request (RAR) is sent for resources in place P15 as outlined in the SLS. Indeed, the length of CAC queue in place P14 is updated and a Resource Allocator Answer (RAA) is returned to the host requesting to confirm whether the network has been reconfigured successfully according to the SLS in place P16 without exceeding the maximum length of queue indicated by the 'M' on the inhibitor arc between P14 and T8. Once the call is admitted via transition T9 and T10 with an immediate priority, and there is enough available resources for a call, transition T11 and T12 are enabled. The reservation strategy need to be determined somewhere by the corresponding pdf rates of λ 4 and λ 5 , respectively to decide whether it is advantageous to admit an incoming call or not. Hence, the decision process for the resource allocation starts at places P8 and P9 after have been admitted by the CAC, and the place P16 is updated and the necessary resources are reserved. Finally, packets are sent to the buffer of the next DiffServ boundary nodes (places P10 and P11), and the traffic is transferred between the different nodes of the network with the respect of each traffic type. Transition T21 represents the bounder node which is directly connected to the host receiving packets.

Once data are received by the middleware, the inverse process starts and the application is notified. For example, in Figure 2 the receiver is notified to receive multimedia stream including (audio and video) via the codec buffer. The information about the audio and video layers with stream identification and synchronization information are essential to the decoding and the subsequent rendering for each of them. So, before being presented at the upper layer, represented by the place Receiver, the synchronization have to be achieved.

Analyzes and discussion

As we outlined before in this paper, concise specification of new marking probabilities can be assisted by numerical analysis. These probabilities determine the mechanisms by which a transition removes token from random set of its input places and deposit tokens in random set of its output places when it fires. In addition to definition given in section 2, we introduce the following: Let define by G the finite or infinite set of marking in our GSPN model. For

s ∈ G;s = (s 1 , s 2 , ..., s n ) (2)
where s j is the number of tokens in place P j ∈ P. The GSPN model of Figure 2 is called K-bounded k ≥ 1 and verify the following condition:

max(s 1 , s 2 , ..., s n ) ≤ k (3)
the marking is 's', that is, we note: Thus, the token count in places P j is never exceeds k. Let T(s) a set of the enabled transitions in Figure 2, when T (s) = {t ∈ T, s j ≥ 1; f or P j ∈ I(t) and s j = 0 f or P j ∈ L(t)} (4) In Figure 2, a transition t ∈ T-T(s) is disabled where the marking is 's', and:

G(t) = {s ∈ G : t ∈ T (s)} (5)
We define S and S' the timed and the immediate markings, respectively:

S = s ∈ G : T (s) ∩ T = 0 (6) S = G -S = s ∈ G : T (s) ∩ T = 0 (7)
Since our GSPN model verify these conditions, any one of its markings changes when one or more enabled transitions fire. Let:

T * ⊆ T (s), let p(s , s, T * ) ( 8
)
where p is the probability of new marking s' given the marking s and the transition set T * , we can easily note the following:

s j -∑ t∈T * 1 I(t * ) (P j ) ≤ s ≤ s j + ∑ t∈T * 1 I(t * ) (P j ) (9)
where j= 1..n, and 1 is the indicator function of the set A.

Global analysis: the PN model is K-bounded and deadlock free: this model works properly, and it takes into account the losses that can occur during communication, and the resulting behavior is live.

Local analysis: consider the middleware bloc in Figure 2, a token in place P0 corresponds to the job (message) in the system, the firing of the transition T0 corresponding to the event "arrival of message" and firing of timed transition T4 corresponds to the event "completion of the service" within the internal message queue of the middleware, and message is scheduled to be sent to the network. The middleware bloc (MIDD) can be modeled as in Figure 3:

10) p((s 1 -1, s 2 -1); s 1 t 0 = ρ (
f or s = (s1, s2) ∈ G

In this manner, the GSPN scheme modeling the message queue capturing the feedback mechanisms in the network queue, fulfilling the stability condition in equation 13 can be assumed:

α = ρ 1 -ρ ≤ 1 (13)
Furthermore, We have shown the SIP bloc is able to fulfill the end-to-end QoS signaling and providing the required QoS as needed by the protocol. When considering the control plane using the same formalism, it is easy to analyze the reservation strategy chooser by the CAC: preparing resources to be used by an active call take some time to be done, corresponding to the firing rates in transitions T6, T7 and T8 given in Table 4. As a result, the PN model is k-bounded, deadlock free, and live.

RELATED WORKS

This Section compare our activities on Middleware-based QoS provisioning with related works. QoS management in reflective middleware: current research issues include adaptive and reflective middleware, and middleware for mobile and ubiquitous systems. For example, [START_REF] Capra | Reflective Middleware Solutions for Context-Aware Applications[END_REF] presents a programmable networking approach to provide QoS as component based architecture. In our approach, we focus on how to use existing middleware and mapping them to provide the QoS for existing application. For instance, HLA and DDS-based existing applications used in large scale networks suffer from the lack of QoS provisioning in multidomain networks like Internet.

Network QoS broker in middleware:Previous works focusing on integrating the signaling process into the QoS provisioning mechanisms. [START_REF] Teodora | A Session Initiation Protocol based Middleware for Multi-Application Management[END_REF] focuses on message based signaling middleware for the control plane to offer per-class QoS. In real word, it does not suffice to deploy router to guarantee the QoS between two hosts at the control plane. likewise, [START_REF] Chi | A user-centric network communication broker for multimedia collaborative computing[END_REF] presented a network communication broker to provide per-class QoS for multimedia collaborative applications. This work, however, add an interface to the application and the middleware for the QoS notification when an event occur in the network and the application should be adapted to this modification. This approach has the advantage of working in plug and play like approach: interface can be used to adapt the application to the network context, and with a minimum of modification this model using simple signaling approach can provide also the 'On-demand' QoS. Middleware-based Network QoS management: Earlier works [START_REF] Dasarathy | Network QoS Assurance in a Multi-Layer Adaptive Resource Management Scheme for Mission-Critical Applications using the CORBA Middleware Framework[END_REF][9] on integrating QoS with the middleware focused on how to add layer3 and layer2 services for CORBA based communication. This approach can be easily deployed in mono-domain network, where only one administrative domain can manage the whole network. If extending it to the internet, the Class of traffic specified at each host will be dropped by the network infrastructure (core and access routers) of other domains. Hence, the QoS provisioning should be enhanced by SIP signaling protocol to enforce resources reservation process during the signaling phase and before the data transmission. The goal of the model presented in this paper is twofold: at the service plane, application signaling using SIP allows the sender to contact the receiver to obtain its IP address and to agree the codec to be used; at the control plane, the network QoS provisioning aims to translate the application requests embedded within the SIP messages , via the middleware, to the network layer and coordinate the endto-end path management and to enforce the end-to-end network resource reservation. Table 5 of the related works presented in this Section vs. our approach. R1 refers to the QoS management in reflective middleware, R2 refers to the Network QoS broker in middleware, and R3 refers to the middleware-based network QoS management. Despite the adaptive and reflective approaches in R1 and R2 work well today when they receive the all the resources required, however, it has been shown that over-provisioning resources has several drawbacks in scalability and fail completely under the slightest anomaly [START_REF] Schmidt | Middleware R&D Challenges for Distributed Real-time and Embedded Systems[END_REF][START_REF] Schmidt | R&D Advances in Middleware for Distributed, Real-time and Embedded Systems[END_REF]. Further, in R3 the QoS provisioning is considered in mono-domain architec-ture. From the expectation of this Table, our approach aims to put all approaches together to provide adaptive framework which enables the end-to-end QoS provisioning over multidomains independently for the underlying transport technologies. More, the stochastic delays in timed transitions allow the diverse DRE flows selecting their QoS requirements. Although the architecture generally works with standard SIP, it does not explicitly deal with QoS provisioning in core networks.

CONCLUSION

We have proposed using GSPN modeling a new model for QoS provisioning based middleware for multimedia and real time application. The following is the summary of the lessons learned: (1) A PN model based approach provides a high level view and the key concepts actions and protocol to be used to coordinate the different technologies for better QoS provisioning in multi-domain infrastructure, and probabilistically time delay is introduced transitions to allow performance evaluation and resolving scheduling problems in DRE. (2) the functional description of this model is based on three concepts: Provisioning, Invocation and operation. (3) For middleware based communication, both per-class and per-flow are needed to support the QoS on demand: after signaling and control phases, the data should be transferred using the appropriate transport protocol. However, this approach still need further research. for example, when including virtual data path with defined QoS parameters, which is the future direction of this work.

  Definition GSPN=(P, T, I, O, H, M 0 ,Π,W) is the Petri Net comprising the following: P is a set of places and T is a set of Transitions such that P ∩ T= . Let introduce a (possibly empty) set T' ⊂ T of immediate transitions. Transitions in T-T' are called timed transition. I(t), O(t), H(t): T×P-→ N, are the input, output and inhibitor function, respectively, for each transition t ∈ T.

Figure 1 .

 1 Figure 1. Global system architecture

Figure 2 .

 2 Figure 2. GSPN Model architecture

Figure 3 .

 3 Figure 3. Middleware Queue model

Table 1 .

 1 Specification of the places in SIP bloc of Figure2

	Place 1 Semantics
	PSIP1	Wait for UA2 rcv replay
	PSIP5	Wait for UA2 rcv replay
	PSIP2	Wait for ACK from UA2
	PSIP6	Wait for ACK from UA1
	PSIP3	Wait for Session connection with UA2
	PSIP7	Wait for Session connection with UA2
	PSIP4	Wait for Termination with UA2
	PSIP8	Wait for termination with UA1
	PSIP9	Invite() initialized from UA1
	PSIP10 Ringing() replay from UA2
	PSIP11 OK response from UA2
	PSIP12 ACK response from UA1, session initiated,
	PSIP13 BYE message send from UA2,
		it decides to end the session
	PSIP14 200 OK, the replay from UA1

Table 2 .

 2 Specification of the transitions in SIP bloc

	Transitions 1 Semantics 2
	T221	UA1 sends connection Request to UA2
	T22	UA2 receive the Request
	T222	UA2 sends Replay
	T223	UA1 receives Replay
	T224	UA2 sends response to establish
		the connection
	T225	UA1 sends ACK message
	T18	Connection established and end-2-end
		path established, data will be transmitted
		by the data plane
	T227	BYE message sent from UA1
	T227	BYE message received by UA2
	T19

Table 3 .

 3 Specification of the places in CAC bloc

	Transitions Semantic
	P4	Call admission
	P5	Call from flow 1
	P6	Call from flow 2
	P14	CAC Queue
	P15	Available resources
	P16	Reserved resources

Table 4 .

 4 Specification of the transitions in CAC bloc

	Transitions Semantics	Rates
	T6	mean duration of timeout of flow 1 λ 1
	T7	mean duration of timeout of flow 1 λ 2
	T8	mean duration of timeout of CAC	λ 3
	T9	Flow1 Call Forward	1
	T10	Flow2 Call Forward	1
	T11	Reservation Process of Flow 1	λ 4
	T12	Reservation Process of Flow 2	λ 5

Table 5 .

 5 illustrates a comparison Comparison of the different approaches

	Features	R1 R2 R3 Our approach
	Application Adaptability
	QoS specification
	Resource allocation
	QoS guarantee
	Multi-domain QoS

Acknowledgment

This research is supported by the French FUI-DGE (Single Inter-Ministerial Fund of the Directorate General for Enterprise) program within the network simulation Platform (PLATSIM).