H Boudaoud

E M Daya

S Belouettar
email: salim.belouettar@tudor.lu

L Duigou

M Potier-Ferry

Damping analysis of beams submitted to passive and active control

In this paper an analytical method is proposed for damping analysis of sandwich beams with piezoelectric and viscoelastic layers. Based on the classical zig-zag model and some assumptions about the electric field, this method leads to an analytical expression of the modal loss factor and frequency. Considering two feedback control laws, the obtained hybrid damping of the sandwich beam is characterized. Numerical finite element applications are considered, in order to examine the efficiency and limitations of the presented method.

Introduction

Vibration suppression of structures is a problem in engineering science that has occupied researchers for a while. The methodology used can be categorized into two groups, namely passive and active controls. In passive control, the material properties of the structure, such as damping and stiffness, are modified so as to change the response of the structure. More than 40 years ago, it was well established that an efficient way to increase passive damping is the use of a sandwich construction with alternating elastic and viscoelastic layers [START_REF] Kerwin | Damping of flexural waves by a constrained viscoelastic layer[END_REF][START_REF] Ross | Damping of plate flexural vibrations by means of viscoelastic laminae[END_REF]. Therefore, many analytical and numerical methods have been proposed to determine the damping parameters (loss factor and frequency) [START_REF] Ambartsumian | On a theory of bending of anisotropic plates[END_REF][START_REF] Carrera | Historical review of zig-zag theories for multilayered plates and shells[END_REF][START_REF] Hu | Review and assessment of various theories for modeling sandwich composites[END_REF]. So, these modal parameters can be obtained in different ways, i.e., by the direct frequency response method [START_REF] Soni | Finite element analysis of viscoelastically damped sandwich structures[END_REF], the complex eigenvalue method [7], the modal strain energy method [START_REF] Soni | Finite element analysis of viscoelastically damped sandwich structures[END_REF][START_REF] Rikards | Finite element analysis of damping the vibrations of laminated composites[END_REF], and the nonlinear complex eigenvalue techniques [START_REF] Chen | Damping prediction of sandwich structures by orderreduction-iteration approach[END_REF][START_REF] Daya | A numerical method for nonlinear eigenvalue problems: Application to vibrations of viscoelastic structures[END_REF][START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. Application to vibrations of viscoelastic shells[END_REF].

In active control, the structural response is controlled by adding an external effort to the structure. For two decades many researchers tried to gain benefit from coupling both control types in the so-called ''hybrid vibration control''. A detailed review and assessment of hybrid damping treatments of structures, with focus on smart structures with piezoelectric control actuation, has been given by Rao [START_REF] Rao | Recent advances in sensing and control of flexible structures via piezoelectric material technology[END_REF]. Since the 1950 s, there have been a lot of papers published on sandwich structures using passive damping. Excellent reviews on the evolution of modeling passive damping were given by Nakra [START_REF] Nakra | Vibration control with viscoelastic materials[END_REF][START_REF] Nakra | Vibration control with viscoelastic materials-ii[END_REF][START_REF] Nakra | Vibration control with viscoelastic materials-iii[END_REF], Torvik [START_REF] Torvik | The analysis and design of constrained layer damping treatments[END_REF], Mead [START_REF] Mead | A comparison of some equations for the flexural vibration of damped sandwich beams[END_REF], and Hu [START_REF] Hu | Review and assessment of various theories for modeling sandwich composites[END_REF]. More recently, many researchers, as Lee [START_REF] Lee | Theory of laminated piezoelectric plates for the design of distributed sensors/actuators. Part i: Governing equations and reciprocal relationships[END_REF], Crawley [START_REF] Crawley | Induced strain actuation of isotropic and anisotropic plates[END_REF], and Preumont [START_REF] Preumont | Vibration control of active structures[END_REF], have provided governing equations, reciprocal relationships and experimental verification of the induced-strain actuation, for smart structures dealing with piezoelectric laminates and active control. Nevertheless, a widely made assumption is the use of a linearly varying electric field [START_REF] Kapuria | An efficient coupled zigzag theory for dynamic analysis of piezoelectric composite and sandwich beams with damping[END_REF][START_REF] Maurini | On a model of layered piezoelectric beams including transverse stress effect[END_REF][START_REF] Liew | Finite element method for the feedback control of fgm shells in the frequency domain via piezoelectric sensors and actuators[END_REF]. Concerning hybrid control, a large part of studies have highlighted structures with a hybrid arrangement called ACLD (SCLD) for active (smart) constrained layer damping [START_REF] Park | Vibration control of bending modes of plates using active constrained layer damping[END_REF][START_REF] Varadan | Closed loop finite-element modeling of active/passive damping in structural vibration control[END_REF][START_REF] Balamurugan | Finite element formulation and active vibration control study on beams using smart constrained layer damping (scld) treatment[END_REF][START_REF] Ray | Active control of laminated cylindrical shells using piezoelectric fiber reinforced composites[END_REF]. In many theoretical and experimental investigations the active system increases the passive transverse shear deformation of the viscoelastic constrained layer and appears to be an effective means to control vibration [START_REF] Vasques | Combined feedback/feedforward active control of vibration of beams with ACLD treatments: Numerical simulation[END_REF]. From an analytical point of view, two-and three-dimensional solutions, in the case of simply supported sandwich beam and plates with piezoelectric layers, have been given by Heyliger [START_REF] Heyliger | Exact solutions for laminated piezoelectric plates in cylindrical bending[END_REF], Vel [START_REF] Vel | Analytical solution for the cylindrical bending vibration of piezoelectric composite plates[END_REF], Cupial [START_REF] Cupial | Vibration and damping analysis of three-layered composite plate with viscoelastic mid-layer[END_REF], and Franco Correia [START_REF] Correia | Modelling and design of adaptive composite structures[END_REF].

From this literature review there is no method for modal parameter calculations of structures with hybrid damping treatments. Indeed, when the material properties are viscoelastic and piezoelectric, the stiffness becomes complex and depends non-linearly on the vibration frequency and on the active control. This method is similar to the one developed for piezoelectric/elastic/piezoelectric laminates by Duigou [START_REF] Duigou | Equivalent stiffness and damping of sandwich piezoelectric beams submitted to active control[END_REF] in which the electrode effect has not been taken into account. In this paper an analytical analysis is proposed when the starting model is a piezoelectric/elastic/viscoelastic/elastic/piezoelectric sandwich beam. The structure is modeled with the zig-zag theory. Next the electromechanical problem is reduced to a purely mechanical problem using elastic and viscoelastic equivalent stiffness properties. The control strategy is based on a proportional direct feedback law or a velocity derivative feedback law. The effect of viscoelastic material damping and the effect of the piezoelectric material properties on modal vibration reduction are studied. The results are validated by comparison with a two-dimensional finite element simulation using the code ABAQUS [34].

An electromechanical model for adaptive laminated beams

In this part, a model for a piezoelectric/elastic/viscoelastic/ elastic/piezoelectric beam is presented. The geometry of the five-layer beam is depicted in Fig. 1. In view of analytical studies, this model is designed to be as simple as possible, including Euler-Bernoulli and Timoshenko kinematics and the assumption of transverse electric field. More precisely, each pair of elastic-piezoelectric layers is represented as an Euler-Bernoulli beam and the whole structure by the zig-zag theory to account for the shear in the central layer [START_REF] Carrera | Historical review of zig-zag theories for multilayered plates and shells[END_REF]. Various levels of approximations are discussed, as well as the treatment of the electrostatic problem in the sensor.

Kinematic models

To describe the properties of a sandwich piezoelectric/elastic/ viscoelastic/elastic/piezoelectric beam, Euler-Bernouilli assumptions are considered for the stiff layers and Timoshenko ones for the central one. The five layers have the same deflection w and the stiff layers have the same rotations. So, the axial displacements in the different layers can be written in the following form:

U i (x, z) = u i (x) -(z -z i )w (x) i = s, e 1 , e 2 , a U v (x, z) = u(x) + zβ(x). (1)
The origin of z axis is the medium plane and z i indicates the medium axis ordinate of each layer:

z s = -z a = (h v + 2h e + h p )/2, z e 1 = -z e 2 = (h v + h e )/2. ( 2
)
u i represents the axial displacement of the ith layer and u that of the central viscoelastic layer. β is its rotation.

The displacement continuity conditions at the interfaces allow us to obtain the axial displacement field in terms of the axial displacement at the center u, of the rotation β and of the deflection w. So, the axial strains of the different layers can be expressed:

                     ε 1s = u + h v 2 β -Hw -(z -z s )w ε 1e 1 = u + h v 2 β - h e 2 w -(z -z e 1 )w ε 1v = u + zβ ε 1e 2 = u - h v 2 β + h e 2 w -(z -z e 2 )w ε 1a = u - h v 2 β + Hw -(z -z a )w (3)
where

H = h p 2 + h e .

Layer generalized constitutive equations

Let us assume that the elastic and viscoelastic materials are linear, homogenous and isotropic. Each layer is considered as an Euler-Bernoulli beam, except the central layer, where Timoshenko kinematics and constant shear stress assumption hold. According to Hooke's law and to the kinematical relations (Eqs. (1) and ( 2)), one defines the axial force N, the bending moments M of the elastic or viscoelastic layers, and the shear force T of the core, as follows:

N e 1 = E e S e u + h v 2 β - h e 2 w M e1 = E e I e 1 w N e 2 = E e S e u - h v 2 β + h e 2 w M e2 = E e I e 2 w N v = S v Y * u M v = I v Y * β T = S v 2(1 + ν v ) Y * ( ẇ + β) (4)
where * denotes the convolution product, Y is the relaxation function, E i is the Young's modulus, ν i the Poisson ratio, I i the quadratic moment, and S i the cross-section area of the considered ith layer. Piezoelectric materials have the capacity to act or react to an electric field or a strain. This property results in the following constitutive equation:

σ = cε -e t E D = eε + E (5)
where σ, ε, D, and E are the stress and strain tensors, electric displacement vector, and field vector, respectively. c, e, and are the elasticity matrix, the piezoelectric matrix and the dielectric permittivity. The coefficients of the elasticity matrix are measured at a constant electric field and the components of the permittivity matrix are measured at constant strain.

The orthotropic piezoelectric materials are poled in the thickness direction and the extension mechanism is considered. The constitutive equations can be represented in terms of a reduced number of components with respect to 1-2-3 axes of the piezoelectric material. We suppose that the piezoelectric layers are thin, the stress tensor is uniaxial, and the electric fields D and E are parallel to the transverse direction. So, the constitutive equations are reduced as in De [35] to

σ 1 D 3 = c * 11 -e * 31 e * 31 * 33 ε 1 E 3 ( 6
)
where .

             * 33 = 33 +
The relation between the electric field and the electric potential can be reduced to the following:

E 3 = - ∂φ(x, z) ∂z . (7)
In this analysis, we consider that the central layer is conductive with a uniform potential, fixed to zero [START_REF] Gopinathan | A review and critique of theories for piezoelectric laminates[END_REF]. Often in the literature, the electric field is assumed to be linearly varying across the thickness, in accordance with the same variation about stress and strain [START_REF] Kapuria | An efficient coupled zigzag theory for dynamic analysis of piezoelectric composite and sandwich beams with damping[END_REF][START_REF] Maurini | On a model of layered piezoelectric beams including transverse stress effect[END_REF][START_REF] Liew | Finite element method for the feedback control of fgm shells in the frequency domain via piezoelectric sensors and actuators[END_REF]. Here, we assume thin piezoelectric layers. This yields a constant electric field through the thickness of each piezoelectric layer:

E 3 = - φ(x) h ( 8
)
where φ(x) is the superficial potential. h is the thickness of the piezoelectric layer. So, from the expression (Eq. ( 3)), one defines the axial forces N and the bending moments M in the piezoelastic layers as follows:

N s = c * 11 S p u + h v 2 β -Hw + e * 31 S p h p φ s (x) M s = c * 11 I s w N a = c * 11 S p u - h v 2 β + Hw -e * 31 S p h p φ a (x) M a = c * 11 I a w with I i = S i (z -z i ) 2 dS i = s, a (9)
where S i is the section area of each ith layer.

Electromechanical constitutive equations of the sandwich beam

In accordance with zig-zag kinematics, the beam stresses can be represented by a global normal force N, by global moments M β and M w and by a shear force. Thus the virtual work of internal forces is in the following form: [START_REF] Daya | A numerical method for nonlinear eigenvalue problems: Application to vibrations of viscoelastic structures[END_REF] where N, M β , M w and T are defined as

δW int = L 0 {Nδu + M w δw + M β δβ + T (δβ + δw )}dx
N = N s + N e1 + N v + N e2 + N a , M β = h v 2 (N e 1 -N e 2 ) + M v + h v 2 (N s -N a ), M w = M s + M e 1 + M e 2 + M a - h e 2 (N e 1 -N e 2 ) -H(N s -N a ) T = S v 2(1 + ν v ) Y * ( ẇ + β). (11)
with

N i = S i σ i dS, M i = S i (z -z i )σ i dS. ( 12
)
By neglecting the rotational inertia, the Principle of d'Alembert can be written as in [START_REF] Duigou | Modélisation numérique de l'amortissement passif et actif des tôles sandwichs comportant des couches viscoélastiques ou piezoélectriques[END_REF] L

0 Nδu + M w δw + M β δβ + T (δβ + δw ) + (ρS) eq ∂ 2 u ∂t 2 δu + ∂ 2 w ∂t 2 δw dx = 0 ∀δu, δw, δβ (13)
where

(ρS) eq = ρ s S s + ρ e1 S e1 + ρ v S v + ρ e2 S e2 + ρ a S a . (14)
ρ i and S i are respectively the density and the section area of layer i. We consider a viscoelastic core that is thin and soft. So, M v and N v can be neglected [START_REF] Rao | Frequency and loss factor of sandwich beams under various boundary conditions[END_REF]. By combining Eqs. ( 11), ( 4) and ( 9), the global electromechanical constitutive equations of the laminate are deduced:

                     N = (ES) mec u + e * 31 S p h p (φ s (x) -φ a (x)) M β = (EI) β β + γ w w + e * 31 h v 2 S p h p (φ s (x) + φ a (x)) M w = γ w β + (EI) mec w -e * 31 H S p h p (φ s (x) + φ a (x)) T = S v 2(1 + ν v ) Y * ( ẇ + β) (15)
with

               (
(EI) β = h 2 v 2 (E e S e + c * 11 S p ). (16)
Notice that there are no membrane/bending coupling terms because of the symmetrical multi-layer construction.

Condensation of electric sensor potential

We suppose no electric displacement in the piezoelectric sensor and the displacement is supposed constant through-the-thickness, leading to

D 3 (z s ) = e * 31 * ε 1s (z s ) - * 33 * φ s (x) h p = 0. ( 17
)
The sensor potential φ s (x) is obtained as a function of the kinematic variables:

φ s (x) = e * 31 * 33 h p u + h v 2 β -Hw . ( 18
)
Considering the relation (Eq. ( 18)) one eliminates the sensor electric potential in Eq. ( 15). This leads to a modified constitutive law:

             N = (ES) eq u + α h v 2 β -αHw -e * 31 S p h p φ a (x) M β = α h v 2 u + (EI) βeq β + γ weq w + e * 31 h v 2 S p h p φ a (x) M w = -αHu + γ weq β + (EI) eq w -e * 31 H S p h p φ a (x) (19)
with

             (ES) eq = (ES) mec + α (EI) eq = (EI) mec + αH 2 γ weq = γ w - h v 2 Hα (EI) βeq = (EI) β + h v 2 2 α α = e 2 31 * 33 * S p . (20)
Two coupling terms appear in Eqs. [START_REF] Crawley | Induced strain actuation of isotropic and anisotropic plates[END_REF] and [START_REF] Preumont | Vibration control of active structures[END_REF]. The α term in γ weq represents the bending/rotation coupling term of electric origin. The membrane/bending coupling term (-αH in N and M w ) is due to the condensation of the sensor electric potential. Formulas (Eq. ( 19)) can be expressed in the following matrix form:

N M β M w = [C]   u β w   + e * 31 S p h p φ a (x)[V a ] (21) with [C] =     (ES) eq α h v 2 -αH α h v 2 (EI) βeq γ weq -αH γ weq (EI) eq     and [V a ] =    -1 h v 2 -H    . (22)
The [C] matrix is real and symmetric. It brings out the membrane/bending coupling term due to the condensation of the electrostatic equation in the sensor electric potential.

Accounting for electrode effect

We consider now a practical case which supposes that two electrodes are bonded on the upper and lower piezoelectric layers. These electrodes lead to an additional equipotentiality boundary condition. Considering as formerly that there is no electric displacement in the piezoelectric sensor, and using the equipotentiality condition, one gets a new system of equations:

L 0 D 3 (z s )dx = L 0 e * 31 * ε 1s (z s ) - * 33 * φ s (x) h p dx = 0 ∂φ s (x) ∂x = 0. ( 23
)
The sensor potential φ s (x) is obtained here as a mean of the kinematic variables:

φ s = e * 31 * 33 h p L 0 u + h v 2 β -Hw dx. ( 24
)
Considering the relation (Eq. ( 24)), one eliminates the sensor electric potential in Eq. [START_REF] Nakra | Vibration control with viscoelastic materials-iii[END_REF]. This leads to a modified constitutive law:

                             N = (ES) mec u + α L 0 u + h v 2 β -Hw dx -e * 31 S p h p φ a M β = α h v 2 L 0 u + h v 2 β -Hw dx + (EI) β β + γ w w + e * 31 h v 2 S p h p φ a M w = -αH L 0 u + h v 2 β -Hw dx + γ w β + (EI) mec w -e * 31 H S p h p φ a . (25)

Accounting for feedback control

We will consider two different feedback control laws. The first one is a direct proportional control law:

φ a (x) = G d φ s (x) (26)
where G d is the direct control gain value. Using this control law leads to the following form for the constitutive law:

N M β M w = ([C] + g[V a ].[V s ] t )   u β w   (27) with [V s ] =    1 h v 2 -H    and g = G d α. ( 28
)
The second feedback law is a negative velocity feedback control law:

φ a (x) = -G v φs (x) (29)
where G v is the velocity control gain value. Considering (Eq. ( 18)), the feedback control law leads to the following form for the constitutive law (here g = G v α):

N M β M w = [C]   u β w   -g[V a ].[V s ] t   u β ẇ   . (30)
Finally, equivalent beam constitutive models (Eqs. ( 27) and ( 30)) are derived, that account first for an analytical solution of the electric problem in the sensor layer and then for the control law in the actuator layer. The first equation (Eq. ( 21)) is a sort of elastic law, but without any symmetry of the matrix since V a is different from V s . This lack of symmetry is due to the nonconservative character of the control law (Eq. ( 26)). Eq. ( 30) is similar to a Kelvin viscoelasticity law. Nevertheless, the damping matrix g[V a ][V s ] t does not have a definite sign. As an outcome, the dynamical problem may be not well posed as in [START_REF] Duigou | Equivalent stiffness and damping of sandwich piezoelectric beams submitted to active control[END_REF]. The latter consequence does not constitute a shortcoming for law frequency analysis.

Moreover, if one considers the practical case with an electrode effect one gets a similar expression for the proportional control law:

N M β M w = [C elastic ]   u β w   + α[V equi ] L 0 u + h v 2 β -Hw dx (31)
where

[C elastic ] = (ES) eq 0 0 0 (EI) β γ w 0 γ w (EI) mec and [V equi ] = 1 -G d 1 + G d -1 -G d ( 32
)
and for the velocity feedback control law

N M β M w = [C elastic ]   u β w   + [A] L 0 u + h v 2 β -Hw dx + [Γ ] L 0 u + h v 2 β -H ẇ dx (33)
with

[A] = α 1 1 -1 and [Γ ] = αG v -1 1 -1 . ( 34
)

An analytical analyses of damping properties

In order to define the damping properties of the sandwich beam we consider free harmonic vibrations. The unknowns u, w, β, N, M β , M w , and T are assumed in harmonic form (ω is the frequency and can be complex):

u = ue iωt , w = we iωt , β = βe iωt N = Ne iωt , M w = M w e iωt , M β = M β e iωt , T = T e iωt . ( 35
)
In this case, the viscoelastic material constitutive equation can be defined by a complex Young's modulus E v (ω), that is frequency dependent. So, the shear force T is written in the following form:

T = G v (ω)S v (w + β) (36)
where

G v (ω) = E v (ω) 2(1+ν v )
is the complex shear modulus, which depends on ω. The modal problem, which is defined from (Eq. ( 13)), becomes

L 0 N M β M w   δu δβ δw   + T (δβ + δw ) -(ρS) eq ω 2 (uδu + wδw)dx = 0 ∀δu, δw, δβ. (37)
The modal loss factor η and the frequency can be obtained in a classical way [START_REF] Rao | Frequency and loss factor of sandwich beams under various boundary conditions[END_REF]:

ω 2 = 2 (1 + iη). ( 38
)
To get an analytical solution, a simply supported beam is considered (N = 0, w = 0, M w = 0, M β = 0 at the extremities). In this case the closed forms of the modes are

w = W sin(kx) β = B cos(kx) u = U cos(kx) k = nπ L (39)
where W , B, and U are the complex amplitudes and k is the mode number (n = 1, 2, . . .).

Analysis with the direct proportional control feedback

If one considers the direct proportional control law (Eq. ( 26)), the constitutive law (Eq. ( 21)) becomes

N M β M w = [C] + g[V a ].[V s ] t U B kW (-k sin(kx)).
(40)

Condensation of the axial displacement u and the rotation β

As is usual, the axial inertia effect is neglected in this analysis, allowing us to reduce the number of generalized displacements (u, β, w) to one w. From the virtual work principle, after integrating by parts Eq. ( 37), the mechanical conservation laws are derived as

M β -T = 0, (a) M w -T -ω 2 w(ρS) eq ẅ = 0. (b) (41)
By inserting Eqs. ( 15) and (36) into Eqs. (41)(a) and (41)(b), the following equations are obtained in which the unknowns are reduced to U, B, and W .

               [(ES) eq -g]U + h v 2 (α -g)B = H(α -g)kW (a) - h v 2 (α + g)k 2 U - (EI) βeq + g h v 2 2 k 2 + G v (ω)S v × B = G v (ω)S v + γ weq -g h v 2 H k 2 kW . (b) (42)
From Eq. ( 42), the amplitudes U, B, and W can be expressed only in terms of W :

U B kW = [R mode ] kW (43) with [R mode ] = U mode B mode 1 (44)
where

                       U mode = 2 H hv (EI) βeq + g( hv 2 ) 2 k 2 + G v (ω)S v + G v (ω)S v + γ weq -g hv 2 H k 2 (ES)eq-g hv 2 (α-g) (EI) βeq + g( hv 2 ) 2 k 2 + G v (ω)S v -hv 2 (α + g)k 2 B mode = G v (ω)S v + γ weq + hv 2 H α 2 -g 2 (ES)eq-g -g k 2 -G v (ω)S v + hv 2 2 α 2 -g 2 (ES)eq-g -g -(EI) βeq k 2 . ( 45
)
Since the shear modulus G v (ω) appears in its expression, [R mode ] is complex and depends on ω.

Complex eigenfrequency

To derive the eigenfrequency, the variational Eq. ( 37) is used, in conjunction with the condensation of u and β. The following virtual displacements are considered:

δu = u = U cos(kx), δw = w = W sin(kx), δβ = β = B cos(kx) ( 46
)
where X is the conjugate complex of X . So, introducing Eqs. ( 43), (44), ( 36) and (46) into Eq. (37) provides the following expression:

ω 2 = k 4 [R mode ] t [C] + g[V a ].[V s ] t [R mode ] + k 2 S v G v (ω)|B mode + 1| 2 (ρS) eq
. (47)

Approximation of the complex frequency

The formula (Eq. ( 47)) does not provide an explicit expression of the eigenfrequency since ω depends on the mode itself and the shear modulus. For simplification and clarity, we assume that the shear modulus of the central layer is constant. These assumptions, often used in the open literature, are not realistic, but could be used to derive the first estimation of the eigenfrequency. Another approximation can be performed when neglecting the influence of the material damping η v in the vector [R mode ], i.e. by replacing G v (ω) by G v R in Eqs. ( 45) and ( 44) and then the mode is real and one gets the following expression of the eigenfrequency and of the loss factor:

       2 = k 4 [R mode ] t [C] + g[V a ].[V s ] t [R mode ] + k 2 S v G v R |B mode + 1| 2 (ρS) eq η η v = k 2 S v G R v (1 + B mode ) 2 2 (ρS) eq . ( 48
)
Notice here that the modal damping is proportional to the viscoelastic damping and that the control gain g can modify this damping through the mode (B mode ) and the pulsation .

Accounting for electrode effect

If one considers the electrode effect with the direct proportional control law (Eq. ( 26)), the constitutive law (Eq. ( 31)) becomes this time

N M β M w = [C elastic ] U B kW (-k sin(kx)) + α[V equi ]Φ (49) with Φ = U + h v 2 B -HkW ((-1) n -1). ( 50
)
The condensation of the axial displacement u and the rotation β lead this time to a new form of the amplitude U, B, and W defined by Eqs. ( 43) and ( 44) and with U mode , B mode defined as

     U mode = bf -ec ea -bd B mode = af -dc bd -ae (51) with                                  a = 1 2 Lk 2 (ES) mec + (1 -G c )α((-1) n -1) 2 L b = h v (1 -G c )α((-1) n -1) 2 2L c = -H(1 -G c )α((-1) n -1) 2 L d = h v (1 + G c )α((-1) n -1) 2 2L e = 1 2 L(k 2 (EI) β + G r S v ) + h 2 v (1 + G c )α((-1) n -1) 2 4L f = 1 2 L(k 2 γ w + G r S v ) - Hh v (1 + G c )α((-1) n -1) 2 2L . ( 52
)
Using a similar procedure as in Section 3.1.2, one can derive the eigenfrequency by introducing Eqs. ( 43), ( 44), ( 36), ( 51) and (46) into Eq. [START_REF] Duigou | Modélisation numérique de l'amortissement passif et actif des tôles sandwichs comportant des couches viscoélastiques ou piezoélectriques[END_REF]. One gets the following expression:

ω 2 = 2k 2 (gU mode + hB mode + i) (ρS) eq L ( 53
)
where

             g = -H(1 + G c )α((-1) n -1) 2 L h = 1 2 L(k 2 γ w + G r S v ) - Hh v (1 + G c )α((-1) n -1) 2 2L i = 1 2 L(k 2 (EI) mec + G r S v ) + H 2 (1 + G c )α((-1) n -1) 2 L .
(54)

Analysis with velocity feedback control

Considering the velocity derivative control law (Eq. ( 29)), the constitutive law (Eq. ( 30)) becomes

N M β M w = [C] -iωg[V a ].[V s ] t U B kW (-k sin(kx)).
(55)

Condensation of the axial displacement u and the rotation β

Considering the formerly defined constitutive law, and using a similar procedure as defined in Section 3.1.1, one gets the new form of amplitude U, B, and W defined by Eq. ( 44), (43) and with U mode , B mode defined by

                 U mode = 2 H hv (EI) βeq -iωg( hv 2 ) 2 k 2 + G v (ω)S v + G v (ω)S v + γ weq + iωg hv 2 H k 2 2 (ES)eq hv α (EI) βeq -iωg( hv 2 ) 2 k 2 + G v (ω)S v -hv 2 (α -iωg)k 2 B mode = G v (ω)S v + γ weq + hv 2 H (α 2 -iωαg) (ES)eq + iωg k 2 -G v (ω)S v + hv 2 2 (α 2 -iωαg) (ES)eq + iωg -(EI) βeq k 2 .
(56)

Complex eigenfrequency

To derive the eigenfrequency, a similar procedure as defined in Section 3.1.2 is used. So, introducing Eqs. ( 43), ( 44) and (46) into Eq. ( 37), and replacing [N, M β , M w ] by Eq. ( 55) and T by Eq. ( 36), provides the following expression of ω. U mode , B mode are defined by Eq. (56):

ω 2 = k 4 [R mode ] t [C] -iωg[V a ].[V s ] t [R mode ] + k 2 S v G v (ω)|B mode + 1| 2 (ρS) eq . (57)
From this equation, three major structural effects are noticeable:

-The elastic effect is defined essentially by the term ''

[R mode ] t ([C]) [R mode ]''
-The passive damping introduced by the viscoelastic layer is expressed by the imaginary part of ''

S v G v (ω)|B mode + 1| 2 ''. -The active damping is expressed by the term ''[R mode ] t (ωg[V a ] .[V s ] t )[R mode ]''.
For small values of passive and active damping, the total damping is the sum of active and passive damping.

Accounting for the electrode effect

If one considers the electrode effect with the velocity feedback control law (Eq. ( 29)), the constitutive law (Eq. ( 33)) becomes this time

N M β M w = [C elastique ] U B kW (-k sin(kx)) + [A] Φ + [Γ ] Φ (58)
where [A] and [Γ ] are defined by Eq. ( 34) and Φ by Eq. (50).

The condensation of the axial displacement u and the rotation β lead to a new form of the amplitude U, B, and W defined by Eqs. (43), ( 44) and (51) with this time in the last equation

                                   a = 1 2 Lk 2 (ES) mec + (1 + iωG v )α((-1) n -1) 2 L b = h v (1 + iωG v )α((-1) n -1) 2 2L c = -H(1 + iωG v )α((-1) n -1) 2 L d = h v (1 -iωG v )α((-1) n -1) 2 2L e = 1 2 L(k 2 (EI) β + G r S v ) + h 2 v (1 -iωG v )α((-1) n -1) 2 4L f = 1 2 L(k 2 γ w + G r S v ) - Hh v (1 -iωG v )α((-1) n -1) 2 2L . (59)
Using a similar procedure as in Sections 3.1.2 and 3.1.4, one can derive the eigenfrequency. Introducing Eqs. ( 43), (44), ( 36), ( 51) and (46) into Eq. ( 37) provides the following expression:

ω 2 = 2k 2 (gU mode + hB mode + i) (ρS) eq L (60)
where

             g = -H(1 -iωG v )α((-1) n -1) 2 L h = 1 2 L(k 2 γ w + G r S v ) - Hh v (1 -iωG v )α((-1) n -1) 2 2L i = 1 2 L(k 2 (EI) mec + G r S v ) + H 2 (1 -iωG v )α((-1) n -1) 2 L . (61)

Results and discussion

For the numerical studies, the geometrical data of the sandwich beam are fixed as follows: h = 10 mm, L = 50 h, h e = 6 h 15 ,

h v = h s = h a = h 15 .
The densities of materials are, for the elastic layers, ρ e = 2040 kg/m 3 , for the viscoelastic core ρ v = 1200 kg/m 3 , and for the piezoelectric layers ρ p = 7500 kg/m 3 .

The elastic shear modulus is fixed to G e = 25 GPa. As for the material properties of the viscoelastic layer, several values will be considered. The constitutive equations of the PZT4 orthotropic piezoelectric layers may be summarized as followed in classical matrix form:

                     σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 D 1 D 2 D 3                      =            C 11 C 12 C 13 0 0 0 0 0 -e 31
           ×                      1 2 3 γ 23 γ 13 γ 12 E 1 E 2 E 3                      (62)
where

             C 11 = C 22 =

Numerical validation of the model

In order to validate the different assumptions of the present analytical model (Section 2), the direct frequency response method with a 2D finite element (FE) in ABAQUS Standard Code is applied to obtain the first modal properties (flexural frequencies and associated damping ratios) of the sandwich beam subjected to a direct proportional feedback control law. Note that this kind of control is the only one available in the ABAQUS code. The sandwich, Fig. 2, is composed of 240 8-node plane strain elements with a reduced integration. A steady-state linear dynamic analysis procedure ABAQUS [34] to monitor the linear response of beam structure subjected to continuous harmonic excitation has been used.

Table 1 presents the values of frequency and loss factor obtained from ABAQUS FE code and those from the present model using the complex [R mode ] (Eq. ( 47)) for different values of g (gain) and G v (shear modulus of the core). In all presented situations, the analytical results (Eq. ( 47)) are very close to the numerical ones.

Validity range of approximated frequency

In order to analyse the impact of considering [R mode ] real (Eq. ( 48)) or complex (Eq. ( 47)) we compare the first eigenvalue, using the software mathematic package Maple. The effects of different viscoelastic cores (shear moduli) as well as different direct control gains were analyzed.

Analysis with no control effect (g = 0)

In this part, the direct piezoelectric control gain is fixed (g = 0).

Results with three different viscoelastic shear modulus G R v and two different viscoelastic material damping (η v = 0.5 and η v = 1 in Table 3) are considered.

It appears clearly that ignoring the complex part of [R mode ]

leads to unsatisfactory results in some identified cases. The layer discrepancies are obtained in the case of a relatively stiff core (G R v /G e = 10 -3 ) or of a large material damping (η v = 1.0). These differences concern the loss factor, while the eigenfrequencies are not very sensitive to the assumption of a real mode. This assumption of a real mode leads to a good approximation if G R v /G e ≤ 10 -4 , η v ≤ 0.5.

Analysis with various direct control gains

In this part, the viscoelastic shear modulus (G v = 2.5 MPa) and the elastic shear modulus (G e = 25 GPa) are fixed. This value of the viscoelastic shear modulus is more or less representative of real materials. For the tested values of the material loss factors (η v = 0.1, 0.5, 1), the assumption of a real mode [R mode ] induces at most an error less than 2%, in the range -20 ≤ g ≤ 50. Out of this range, the system tends to be unstable.

From Fig. 4, one observes that the frequency varies more or less linearly, with respect to the control gain, in the range -10 ≤ g ≤ 40. Indeed, this control law looks like an elastic law.

It modifies mainly the stiffness of the structure and therefore the eigenfrequencies. The modal loss factors η depend on the control only via the changes in the frequencies, and its variation becomes significant only for large values of the control gain g, Fig. 5.

In a second part one considers the equipotentiality condition imposed by an electrode. One can see for the tested values of the material loss factors (η v = 0.1, 0.5, 1) that the effect of the active control is weak. Indeed if one considers value of the control gain g in the range [-1, 1] the modal properties with and without electrode are similar (Tables 1 and2). Then one consider more important values for the control gain g and we test for a fixed value from the loss factor (η v = 0.5) the importance of taking into account the equipotentiality. In Fig. 3 one can clearly see that if |g| is important the equipotentiality condition leads to smaller value of the eigenfrequency.

Analysis with various velocity feedback control gains

In this part, the same moduli (G v = 2.5 MPa and G e = 25 GPa) as in 4.2.2 are considered. From Figs. 6 and7, one sees that this control law has a strong influence on the structural damping properties, but not on the eigenfrequencies. In the considered range 0 ≤ g ≤ 0.025, the loss factor increases linearly with the control gain. Moreover, the curves of Fig. 7 are more or less parallel. This means that the contributions of active and passive damping are nearly additive.

In a second part one considers the equipotentiality condition imposed by an electrode. We tested different values of the material loss factors (η v = 0.1, 0.5, 1) and of the gain control g but we could not observe any stability range.

Conclusion

In this paper, simple models and analytical modal solutions for piezoelectric/elastic/viscoelastic/elastic/piezoelectric beams have been presented. Direct proportional feedback as well as a velocity feedback control law have been considered, and constitutive equations of the whole controlled sandwich beam have been established. The eigenfrequencies and the loss factors have been defined for this hybrid damped system. As expected, the direct control law modifies the stiffness and the frequencies while the second control law influences mainly the damping properties of the system. The influence of the real mode assumption on frequency and modal damping. First vibration mode.

Fig. 1 .

 1 Fig. 1. Notations used for the piezoelectric/elastic/viscoelastic/elastic/piezoelectric beam.

C 12 C 22

 1222

Fig. 2 .

 2 Fig. 2. Dimensions of the five-layer sandwich structure.

Fig. 3 .

 3 Fig. 3. Variation of the frequency with and without equipotentiality condition η = 0.5 direct proportional control g varying. G v = 2.5 MPa, G e = 25 GPa.

Fig. 4 .

 4 Fig. 4. First eigenfrequency as a function of the direct proportional control gain g. G v = 2.5 MPa, G e = 25 GPa.

Fig. 5 .

 5 Fig. 5. The damping η of the first mode as a function of the direct proportional control gain g. G v = 2.5 MPa, G e = 25 GPa.

Fig. 6 .

 6 Fig. 6. First eigenfrequency as a function of the velocity feedback control gain g. G v = 2.5 MPa, G e = 25 GPa.

  ES) mec = 2E e S e + 2c *

11 S p (EI) mec = E e I e1 + I e2 + h 2 e 2 S e + c * 11 (I a + I s + 2H 2 S p ) γ w = -h v 2 (h e S e E e + 2c * 11 S p H)

  139.0 GPa; C 33 = 115.3 GPa C 13 = C 23 = 74.3 GPa; C 12 = 77.8 GPa C 44 = C 55 = 25.6 GPa; C 66 = 30.6 GPa e 31 = e 32 = -5.21 C/m 2 ; e 33 = 15.1 C/m 2 ; e 15 = e 24 = 12.7 C/m 2 11 = 22 = 1.31e -8 F/m; 33 = 1.15e -8 F/m.

	(63)

Table 1

 1 Comparison of the beam model (analytical) with a 2D finite element analysis (numerical)

	g	G v (MPa)	f (Hz)		η
			Analytical	Numerical	Analytical	Numerical
		2.5	67.2	67.2	0.171	0.171
	1	25	86.4	86.5	0.038	0.038
		2500	90.2	90.3	0.00043	0.00042
		2.5	66.2	66.1	0.169	0.169
	0	25	84.8	84.7	0.038	0.039
		2500	88.2	88.1	0.00042	0.00045
		2.5	65.1	65	0.167	0.168
	-1	25	82.9	82.9	0.037	0.037
		2500	86.2	86.2	0.00041	0.00043

η v = 0.5. Direct proportional control law.

Table 2

 2 Comparison

			Analytical	Numerical	Analytical	Numerical
		2.5	66.9	66.8	0.170	0.176
	1	25	85.9	85.9	0.038	0.038
		2500	89.5	89.6	0.00043	0.00042
		2.5	66.0	66.0	0.169	0.172
	0	25	84.4	84.5	0.037	0.037
		2500	87.9	87.8	0.00042	0.00045
		2.5	65.1	65.1	0.167	0.168
	-1	25	83.0	83.0	0.037	0.037
		2500	86.3	86.2	0.00041	0.00042
	η

of the beam model (analytical) with a 2D finite element analysis (numerical) g G v (MPa) f (Hz) η v = 0.5. Direct proportional control law and equipotentiality condition.

Table 3

 3 Purely passive damping(g = 0) R mode ] real [R mode ] complex (%) [R mode ] real [R mode ] complex (%)

	η v	G R v (MPa)	f (Hz) with	f (Hz) with	η with	η with
		0.25	44.0	43.9 (0.2)	0.108	0.107 (0.9)
	0.5	2.5	65.0	66.1 (1.7)	0.167	0.169 (1.2)
		25	84.1	84.8 (0.8)	0.032	0.038 (15.8)
		0.25	43.9	44.2 (0.7)	0.214	0.212 (0.9)
	1.0	2.5	65.6	69.1 (5.1)	0.26	0.273 (4.8)
		25	84.8	85.9 (1.3)	0.028	0.049 (42.9)

[

Acknowledgements

The authors acknowledge the financial support of the Luxembourgish National Research Fund (FNR-MAT-08-01 ADYMA) and the European FP6 ERA-NET Materials (MATERA).


Bernard Berthomieu

Didier Le Botlan

Silvano Dal

Petri Net Reductions for Counting Markings

We propose a method to count the number of reachable markings of a Petri net without having to enumerate these first. The method relies on a structural reduction system that reduces the number of places and transitions of the net in such a way that we can faithfully compute the number of reachable markings of the original net from the reduced net and the reduction history. The method has been implemented and computing experiments show that reductions are effective on a large benchmark of models.

Introduction

Structural reductions are an important class of optimization techniques for the analysis of Petri Nets (PN for short). The idea is to use a series of reduction rules that decrease the size of a net while preserving some given behavioral properties. These reductions are then applied iteratively until an irreducible PN is reached on which the desired properties are checked directly. This approach, pioneered by Berthelot [START_REF] Berthelot | Checking properties of nets using transformations[END_REF][START_REF] Berthelot | Transformations and decompositions of nets[END_REF], has been used to reduce the complexity of several problems, such as checking for boundedness of a net, for liveness analysis, for checking reachability properties [START_REF] Jensen | TAPAAL and reachability analysis of P/T nets[END_REF] or for LTL model checking [START_REF] Esparza | Net reductions for LTL model-checking[END_REF].

In this paper, we enrich the notion of structural reduction by keeping track of the relation between the markings of an (initial) Petri net, N 1 , and its reduced (final) version, N 2 . We use reductions of the form (N 1 , Q, N 2 ), where Q is a system of linear equations that relates the (markings of) places in N 1 and N 2 . The reductions are tailored so that the state space of N 1 (its set of reachable markings) can be reconstructed from that of N 2 and equations Q. In particular, when N 1 is totally reduced (N 2 is then the empty net), the state space of N 1 corresponds with the set of non-negative integer solutions to Q. Then Q acts as a symbolic representation for sets of markings, in much the same way one can use decision diagrams or SAT-based techniques.

In practice, reductions often lead to an irreducible non-empty residual net. In this case, we can still benefit from an hybrid representation combining the state space of the residual net (expressed, for instance, using a decision diagram) and the symbolic representation provided by linear equations. This approach can provide a very compact representation of the state space of a net. Therefore it is suitable for checking reachability properties, that is whether some reachable marking satisfies a given set of linear constraints. However, checking reachability properties could benefit of more aggressive reductions since it is not generally required there that the full state space is available (see e.g. [START_REF] Jensen | TAPAAL and reachability analysis of P/T nets[END_REF]). For particular properties, it could be more efficient to only approximate the set of reachable markings, and therefore to use more aggressive reduction rules (see e.g. [START_REF] Jensen | TAPAAL and reachability analysis of P/T nets[END_REF]).

At the opposite, we focus on computing a (symbolic) representation of the full state space. A positive outcome of our choice is that we can derive a method to count the number of reachable markings of a net without having to enumerate them first.

Computing the cardinality of the reachability set has several applications. For instance, it is a straightforward way to assess the correctness of tools-all tools should obviously find the same results on the same models. This is the reason why this problem was chosen as the first category of examination in the recurring Model-Checking Contest (MCC) [START_REF] Kordon | Complete Results for the 2017 Edition of the Model Checking Contest[END_REF][START_REF] Kordon | MCC'2017 -The seventh model checking contest[END_REF]. We have implemented our approach in the framework of the TINA toolbox [START_REF] Berthomieu | The tool TINAconstruction of abstract state spaces for Petri nets and Time Petri nets[END_REF] and used it on the large set of examples provided by the MCC (see Sect. 7). Our results are very encouraging, with numerous instances of models where our performances are several orders of magnitude better than what is observed with the best available tools.

Outline. We first define the notations used in the paper then describe the reduction system underlying our approach, in Sect. [START_REF] Berthelot | Transformations and decompositions of nets[END_REF]. After illustrating the approach on a full example, in Sect. 4, we prove in Sect. 5 that the equations associated with reductions allow one to reconstruct the state space of the initial net from that of the reduced one. Section 6 discusses how to count markings from our representation of a state space while Sect. 7 details our experimental results. We conclude with a discussion on related works and possible future directions.

Petri Nets

Some familiarity with Petri nets is assumed from the reader. We recall some basic terminology. Throughout the text, comparison (=, ≥) and arithmetic operations (-, +) are extended pointwise to functions.

A marked Petri net is a tuple N = (P, T, Pre, Post, m 0 ) in which P , T are disjoint finite sets, called the places and transitions, Pre, Post : T → (P → N) are the pre and post condition functions, and m 0 : P → N is the initial marking. Figure 1 gives an example of Petri net, taken from [START_REF] Stahl | Decomposing Petri net State Spaces[END_REF], using a graphical syntax: places are pictured as circles, transitions as squares, there is an arc from place p to transition t if Pre(t)(p) > 0, and one from transition t to place p if Post(t)(p) > 0. The arcs are weighted by the values of the corresponding pre or post conditions (default weight is 1). The initial marking of the net associates integer 1 to place p 0 and 0 to all others.

A marking m : P → N maps a number of tokens to every place A firing sequence σ over T is a sequence t 1 , . . . , t n of transitions in T such that there are some markings m 1 . . . , m n+1 with m 1

t 1 → m 2 ∧ . . . ∧ m n tn → m n+1 . This can be written m 1 σ → m n+1 . Its displacement, or marking change, is ∆(σ) = Σ n i=1 (Post(t i ) -Pre(t i ))
, where ∆ : P → Z, and its hurdle H(σ) is the smallest marking (pointwise) from which the sequence is firable. Displacements (or marking changes) and hurdles are discussed in [START_REF] Hack | Decidability questions for Petri Nets[END_REF], where the existence and uniqueness of hurdles is proved. As an illustration, the displacement of sequence t 2 t 5 t 3 in the net Figure 1 is {(p 2 , -1), (p 5 , 1), (p 9 , 1)} (and 0 for all other places, implicitly), its hurdle is {(p 2 , 1), (p 3 , 1), (p 7 , 1)}.

The postset of a transition t is t

• = { p | Post(t)(p) > 0 }, its preset is • t = { p | Pre(t)(p) > 0 }. Symmetrically for places, p • = { t | Pre(t)(p) > 0 } and • p = { t | Post(t)(p) > 0 }.
A net is ordinary if all its arcs have weight one; for all transition t in T , and place p in P , we have Pre(t)(p) ≤ 1 and Post(t)(p) ≤ 1. Otherwise it is said generalized.

A net N is bounded if there is an (integer) bound b such that m(p) ≤ b for all m ∈ R(N ) and p ∈ P . The net is said safe when the bound is 1. All nets considered in this paper are assumed bounded.

The net in Figure 1 is ordinary and safe. Its state space holds 14 markings.

The Reduction System

We describe our set of reduction rules using three main categories. For each category, we give a property that can be used to recover the state space of a net, after reduction, from that of the reduced net.

Removal of Redundant Transitions

A transition is redundant when its effects can always be achieved by firing instead an alternative sequence of transitions. Our definition of redundant transitions slightly strengthens that of bypass transitions in [START_REF] Recalde | Improving the decision power of rank theorems[END_REF]. It is not fully structural either, but makes it easier to identify special cases structurally.

Definition 3.1 (Redundant transition). Given a net (P, T, Pre, Post, m 0 ), a transition t in T is redundant if there is a firing sequence σ over T \ {t} such that ∆(t) = ∆(σ) and H(t) ≥ H(σ).

There are special cases that do not require to explore combinations of transitions. This includes identity transitions, such that ∆(t) = 0, and duplicate transitions, such that for some other transition t and integer k, ∆(t) = k.∆(t ). Finding redundant transitions using Definition 3.1 can be convenient too, provided the candidate σ are restricted (e.g. in length). Figure 2 (left) shows some examples of redundant transitions. Clearly, removing a redundant transition from a net does not change its state space.

Removal of Redundant Places

A place is redundant if it never restricts the firing of its output transitions. Removing redundant places from a net preserves its language of firing sequences [START_REF] Berthelot | Transformations and decompositions of nets[END_REF]. We wish to avoid enumerating marking for detecting such places, and further be able to recover the marking of a redundant place from those of the other places. For these reasons, our definition of redundant places is a slightly strengthened version of that of structurally redundant places in [START_REF] Berthelot | Transformations and decompositions of nets[END_REF] (last clause is an equation).

Definition 3.2 (redundant place). Given a net (P, T, Pre, Post, m 0 ), a place p in P is redundant if there is some set of places I from P \ {p}, some valuation v : (I ∪ {p}) → (N -{0}), and some constant b ∈ N such that, for any t ∈ T :

1. The weighted initial marking of p is not smaller than that of I: b = v(p).m 0 (p) -Σ q∈I v(q).m 0 (q) 2. To fire t, the difference between its weighted precondition on p and that on I may not be larger than b: v(p).Pre(t)(p) -Σ q∈I v(q).Pre(t)(q) ≤ b 3. When t fires, the weighted growth of the marking of p is equal to that of I:

v(p).(Post(t)(p) -Pre(t)(p)) = Σ q∈I v(q).(Post(t)(q) -Pre(t)(q))
This definition can be rephrased as an integer linear programming problem [START_REF] Silva | Linear algebraic and linear programming techniques for the analysis of place/transition net systems[END_REF], convenient in practice for computing redundant places in reasonably sized nets (say when |P | ≤ 50). Like with redundant transitions, there are special cases that lead to easily identifiable redundant places. These are constant places-those for which set I in the definition is empty-and duplicated places, when set I is a singleton. Figure 2 (right) gives some examples of such places.

From Definition 3.2, we can show that the marking of a redundant place p can always be computed from the markings of the places in I and the valuation function v. Indeed, for any marking m in R(N ), we have v(p).m(p) = Σ q∈I v(q).m(q) + b, where the constant b derives from the initial marking m 0 . Hence we have a relation k p .m(p) = ρ p (m), where k p = v(p) and ρ p is some linear expression on the places of the net. Theorem 3.2. If N is the result of removing some redundant place p from net N , then there is an integer constant k ∈ N * , and a linear expression ρ, such that, for all marking m: m ∪ {(p,

(1/k).ρ(m))} ∈ R(N ) ⇔ m ∈ R(N ).

Place Agglomerations

Conversely to the rules considered so far, place agglomerations do not preserve the number of markings of the nets they are applied to. They constitute the cornerstone of our reduction system; the purpose of the previous rules is merely to simplify the net so that agglomeration rules can be applied. We start by introducing a convenient notation.

Definition 3.3 (Sum of places).

A place a is the sum of places p and q, written a = p q, if: m 0 (a) = m 0 (p) + m 0 (q) and, for all transition t, Pre(t

)(a) = Pre(t)(p) + Pre(t)(q) and Post(t)(a) = Post(t)(p) + Post(t)(q).
Clearly, operation is commutative and associative. We consider two categories of place agglomeration rules; each one consisting in the simplification of a sum of places. Examples are shown in Fig. 3. Given a net (P, T, Pre, Post, m 0 ), a pair of places p, q in P can be chain agglomerated if there is some t ∈ T such that:

X p t q → X a t t1 t0 t2 Y p1 Z p2 X p0 → t1 t0
• t = {p}; t • = {q}; Pre(t)(p) = Post(t)(q) = 1;
• q = {t}; and m 0 (q) = 0. Their agglomeration consists of replacing places p and q by a place a equal to their sum: a = p q.

Definition 3.5 (Loop agglomeration). A sequence of n places (π i ) n-1 i=0 can be loop agglomerated if the following condition is met:

(∀i < n)(∃t ∈ T )(Pre(t) = {(π i , 1)} ∧ Post(t) = {(π (i+1)(mod n) , 1)}) .
Their agglomeration consists of replacing places π 0 , . . . , π n-1 by a single place, a, defined as their sum: a = n-1 i=0 π i .

Clearly, whenever some place a of a net obeys a = p q for some places p and q of the same net, then place a is redundant in the sense of definition 3.2. The effects of agglomerations on markings are stated by Theorem 3.3. Theorem 3.3. Let N and N be the nets before and after agglomeration of some set of places A as place a. Then for all markings m over (P \ A) and m over A we have:

(m ∪ m ) ∈ R(N ) ⇔ m ∪ {(a, Σ p∈A m (p)} ∈ R(N ).
Proof. Assume N is a net with set of places P . Let us first consider the case of the chain agglomeration rule in Fig. 3 (top). We have to prove that for all marking m of P \ {p, q} and for all values x, y in N:

m ∪ {(p, x), (q, y)} ∈ R(N ) ⇔ m ∪ {(a, x + y)} ∈ R(N )
Left to right (L): Let N + be net N with place a = p q added. Clearly, a is redundant in N + , with v(a) = v(p) = v(q) = 1. So N and N + admit the same firing sequences, and for any m ∈ R(N + ), we have m(a) = m(p) + m(q). Next, removing places p and q from N + (the result is net N ) can only relax firing constraints, hence any σ firable in N + (and thus in N ) is also firable in N , which implies the goal.

Right to left (R): we use two intermediate properties (∀m, x, u, v implicit). We write m ∼ m when m and m agree on all places except p, q and a, and

m ≈ m when m ∼ m ∧ m(p) = m (a) ∧ m(q) = 0. Property (1): m ∪ {(a, x)} ∈ R(N ) ⇒ m ∪ {(p, x), (q, 0)} ∈ R(N ).
Since ∆(t) = 0, any marking reachable in N is reachable by a sequence not containing t, call these sequences t-free. Property (1) follows from a simpler relation, namely (Z):

whenever m ≈ m (m ∈ R(N ), m ∈ R(N )
) and m δ → w , (δ t-free), then there is a sequence ω such that m ω → w and w ≈ w . Any t-free sequence firable in N but not in N can be written σ.t .γ, where σ is firable in N and transition t is not firable in N after σ. Let w, w be the markings reached by σ in N and N , respectively. Since σ is firable in N , we have w ≈ w , by (L) and the fact that σ is t-free (only t can put tokens in q). That t is not firable at w but firable at w is only possible if t is some output transition of a since w ∼ w and the preconditions of all other transitions of N than a are identical in N and N . That is, t must be an output transition of either or both p or q in N . If t has no precondition on q in N , then it ought to be firable at w in N since w(p) = w (a). So t must have a precondition on q; we have w(q) ≥ Pre(t )(q) in N and w (a) ≥ Pre (t )(a) in N . Therefore, we can fire transition t n times from w in N , where n = Pre(t )(q), since w (a) = w(p) and t is enabled at w , and this leads to a marking enabling t . Further, firing t at that marking leaves place q in N empty since only transition t may put tokens in q. Then the proof of Property (1) follows from (Z) and the fact that Definition 3.4 ensures m 0 ≈ m 0 .

Property (2): if m ∪ {(p, x), (q, 0)} ∈ R(N ) and (u + v = x) then m ∪ {(p, u), (q, v)} ∈ R(N ).
Obvious from Definition 3.4: the tokens in place p can be moved one by one into place q by firing t in sequence v times.

Combining Property (1) and ( 2) is enough to prove (R), which completes the proof for chain agglomerations. The proof for loop agglomerations is similar.

The Reduction System

The three categories of rules introduced in the previous sections constitute the core of our reduction system. Our implementation actually adds to those a few special purpose rules. We mention three examples of such rules here, because they play a significant role in the experimental results of Sect. 7, but without technical details. These rules are useful on nets generated from high level descriptions, that often exhibit translation artifacts like dead transitions or source places.

The first extra rule is the dead transition removal rule. It is sometimes possible to determine statically that some transitions of a net are never firable. A fairly general rule for identifying statically dead transitions is proposed in [START_REF] Esparza | Net reductions for LTL model-checking[END_REF]. Removal of statically dead transitions from a net has no effects on its state space.

A second rule allows us to remove a transition t from a net N when t is the sole transition enabled in the initial marking and t is firable only once. Then, instead of counting the markings reachable from the initial marking of the net, we count those reachable from the output marking of t in N and add 1. Removing such transitions often yields structurally simpler nets.

Our last example is an instance of simple rules that can be used to do away with very basic (sub-)nets, containing only a single place. This is the case, for instance, of the source-sink nets defined below. These rules are useful if we want to fully reduce a net. We say that a net is totally reduced when its set of places and transitions are empty

(P = T = ∅). Definition 3.6 (Source-sink pair). A pair (p, t) in net N is a source-sink pair if • p = ∅, p • = {t}, Pre(t) = {(p, 1)} and Post(t) = ∅. Theorem 3.4 (Source-sink pairs). If N is the result of removing a source-sink pair (p, t) in net N then (∀z ≤ m 0 (p))(∀m)(m ∪ {(p, z)} ∈ R(N ) ⇔ m ∈ R(N )).
Omitting for the sake of clarity the first two extra rules mentioned above, our final reduction system resumes to removal of redundant transitions (referred to as the T rule) and of redundant places (R rule), agglomeration of places (A rules) and removal of source-sink pairs (L rule).

Rules T have no effects on markings. For the other three rules, the effect on the markings can be captured by an equation or an inequality. These have shape v p .m(p) = q =p v q .m(q)+b for redundant places, where b is a constant, shape m(a) = Σ p∈A m(p) for agglomerations, and shape m(p) ≤ k for source-sink pairs, where k is some constant. In all these equations, the terms m(q) are marking variables; variable m(q) is associated with the marking of place q. For readability, we will often use the name of the place instead of its associated marking variable. For instance, the marking equation 2.m(p) = 3.m(q) + 4, resulting from a (R) rule, would be simply written 2.p = 3.q + 4.

We show in Sect. 5 that the state space of a net can be reconstructed from that of its reduced net and the set of (in)equalities collected when a rule is applied. Before considering this result, we illustrate the effects of reductions on a full example.

An Illustrative Example -HouseConstruction

We take as example a model provided in the Model Checking Contest (MCC, http: //mcc.lip6.fr), a recurring competition of model-checking tools [START_REF] Kordon | Complete Results for the 2017 Edition of the Model Checking Contest[END_REF]. This model is a variation of a Petri net model found in [START_REF] Lyle | Petri Net Theory and the Modeling of Systems[END_REF], which is itself derived from the PERT chart of the construction of a house found in [START_REF] Levy | Introduction to the criticalpath method[END_REF]. The model found in the MCC collection, reproduced in Fig. 4, differs from that of [START_REF] Lyle | Petri Net Theory and the Modeling of Systems[END_REF] in that it omits time constraints and a final sink place. In addition, the net represents the house construction process for a number of houses simultaneously rather than a single one. The number of houses being built is represented by the marking of place p 1 of the net (10 in the net represented in Fig. 4).

We list in Fig. 5 a possible reduction sequence for our example net, one for each line. To save space, we have omitted the removal of redundant transitions. For each reduction, we give an indication of its kind (R, A, . . . ), the marking equation witnessing the reduction, and a short description. The first reduction, for instance, says that place p 19 is removed, being a duplicate of place p 20 . At the second step, places p11 and p7 are agglomerated as place a1, a "fresh" place not found in the net yet.

Each reduction is associated with an equation or inequality linking the markings of the net before and after application of a rule. The system of inequalities gathered is shown below, with agglomeration places a i eliminated. We show in the next section that the set of solutions of this system, taken as markings, is exactly the set of reachable markings of the net. This example is totally reduced using the sequence of reductions listed. And we have found other examples of totally reducible net in the MCC benchmarks. In the general case, our reduction system is not complete; some nets may be only partially reduced, or not at all. When a net is only partially reducible, the inequalities, together with an explicit or logic-based symbolic description of the reachability set of the residual net, yield a hybrid representation of the state space of the initial net. Such hybrid representations are still suitable for model checking reachability properties or counting markings.

p19 = p20 p4 = p6 + p15 + p11 + p7 p12 = p10 + p8 p9 + p5 = p6 + p8 p13 = p10 + p9 p21 + p18 = p22 + p20 + p23 p27 = p23 + p26 p25 + p16 + p15 + p11 + p7 = p26 + p23 + p22 + p20 + p17 + p14 + p10 + p8 p26 + p23 + p22 + p20 + p17 + p14 + p10 + p8 + p6 + p3 + p2 + p1 ≤ 10
Order of application of reduction rules. Our reduction system does not constrain the order in which reductions are applied. Our tool attempts to apply them in an order that minimizes reduction costs.

The rules can be classified into "local" rules, detecting some structural patterns on the net and transforming them, like removal of duplicate transitions or places, or chain agglomerations, and ''non-local" rules, like removal of redundant places in the general case (using integer programming). Our implementation defers the application of the nonlocal rules until no more local rule can be applied. This decreases the cost of non-local reductions as they are applied to smaller nets.

Another issue is the confluence of the rules. Our reduction system is not confluent: different reduction sequences for the same net could yield different residual nets. This follows from the fact that agglomeration rules do not preserve in general the ordinary character of the net (that all arcs have weight 1), while agglomeration rules require that the candidate places are connected by arcs of weight 1 to the same transition.

An example net exhibiting the problem is shown in Fig. 6(a). Agglomeration of places p3 and p4 in this net, followed by removal of identity transitions, yields the net in Fig. 6(b). Place a1 in the reduced net is the result of agglomerating p3 and p4; this is witnessed by equation a1 = p3+p4. Note that the arcs connecting place a1 to transitions t0 and t1 both have weight 2.

Next, place p2 in the reduced net is a duplicate of place a1, according to the definitions • If p2 is removed instead, then the resulting net cannot be reduced further: places p0, a1 and p1 cannot be agglomerated because of the presence of arcs with weight larger than 1.

t 0 t 1 p 2 p 3 p 4 t 2 t 3 p 0 p 1 (a) t 0 t 1 p 2 a 1 p 0 p 1 2 2 (b)
Confluence of the system could be easily obtained by restricting the agglomeration rules so that no arcs with weight larger than 1 could be produced. But it is more effective to favour the expressiveness of our reduction rules.

Alternatively, agglomeration rules could be generalized to handle arbitrary weights on the arcs linking the agglomerated places; this is a scheduled improvement of our reduction system.

Correctness of Markings Reconstruction

We prove that we can reconstruct the markings of an (initial) net, before application of a rule, from that of the reduced net. This property ensues from the definition of a Case (R) : By Th. 3.2 we have : R(

N 2 ) = R(N 3 ) ↑ P 2 ∩ k.p = l ↑ P 2 . replacing R(N 2 ) by this value in (H) yields R(N 1 ) = ((R(N 3 ) ↑ P 2 ∩ k.p = l ↑ P 2 ) ↑ V ∩ Q ↑ V) ↓ P 1 .
Since P 2 ⊆ V, we may safely lift to V instead of P 2 , so: R(

N 1 ) = (R(N 3 ) ↑ V ∩ k.p = l ↑ V ∩ Q ↑ V) ↓ P 1 . Which is equivalent to: R(N 1 ) = (R(N 3 ) ↑ V ∩ Q; k.p = l ↑ V) ↓ P 1 ,
and equal to (G) since P 3 ⊆ V and Q = (Q; k.p = l).

Case (A): Let S p denotes the value Σ p∈A (p). By Th. 3.3 we have: R(N 2 ) = (R(N 3 ) ↑ (P 2 ∪ P 3 ) ∩ a = S p ↑ (P 2 ∪ P 3 )) ↓ P 2 . Replacing R(N 2 ) by this value in (H) yields: R(N 1 ) = (((R(N 3 ) ↑ (P

2 ∪ P 3 ) ∩ a = S p ↑ (P 2 ∪ P 3 )) ↓ P 2 ) ↑ V ∩ Q ↑ V) ↓ P 1 . Instead of V, we may lift to U since U = V ∪ {a}, a ∈ V(Q) and a ∈ P 1 , so: R(N 1 ) = (((R(N 3 ) ↑ (P 2 ∪ P 3 ) ∩ a = S p ↑ (P 2 ∪ P 3 )) ↓ P 2 ) ↑ U ∩ Q ↑ U) ↓ P 1 .
Projection on P 2 may be omitted since P 2 ∪ P 3 = P 2 ∪ {a} and a ∈ V(Q), leading to:

R(N 1 ) = ((R(N 3 ) ↑ (P 2 ∪ P 3 ) ∩ a = S p ↑ (P 2 ∪ P 3 )) ↑ U ∩ Q ↑ U) ↓ P 1 . Since P 2 ∪ P 3 ⊆ U, this is equivalent to: R(N 1 ) = (R(N 3 ) ↑ U ∩ a = S p ↑ U ∩ Q ↑ U) ↓ P 1 . Grouping equations yields: R(N 1 ) = (R(N 3 ) ↑ U ∩ Q; a = S p ↑ U) ↓ P 1 , which is equal to (G) since Q = (Q; a = S p ).

case (L):

The proof is similar to that of case (R) and is based on the relation R(N 2 ) = R(N 3 ) ↑ P 2 ∩ p ≤ k ↑ P 2 , obtained from Th. 3.4. Theorem 5.1 states the correctness of our reduction systems, since we can compose reductions sequentially and always obtain a net-abstraction. In particular, if a net N is fully reducible, then we can derive a system of linear equations Q such that (N, Q, ∅) is a net-abstraction. In this case the reachable markings of N are exactly the solutions of Q, projected on the places of N . If the reduced net, say N r , is not empty then each marking m ∈ R(N r ) represents a set of markings Q m ⊂ R(N ): the solution set of Q in which the places of the residual net are constrained as in m, and then projected on the places of N . Moreover the family of sets

{ Q m | m ∈ R(N r )} is a partition of R(N ).

Counting Markings

We consider the problem of counting the markings of a net N from the set of markings of the residual net N r and the (collected) system of linear equations Q. For totally reduced nets, counting the markings of N resumes to that of counting the number of solutions in non negative integer variables of system Q. For partially reduced nets, a similar process must be iterated over all markings m reachable in N r (a better implementation will be discussed shortly).

Available methods. Counting the number of integer solutions of a linear system of equations (inequalities can always be represented by equations by the addition of slack variables) is an active area of research.

A method is proposed in [START_REF] Behle | 0/1 vertex and facet enumeration with BDDs[END_REF], implemented in the tool Azove, for the particular case where variables take their values in {0, 1}. The method consists of building a Binary Decision Diagram for each equation, using Shannon expansion, and then to compute their conjunction (this is done with a specially tailored algorithm). The number of paths of the BDD gives the expected result. Our experiments with Azove show that, although the tool can be fast, its performances on larger system heavily depend on the ordering chosen for the BDD variables, a typical drawback of decision diagram based techniques. In any case, its usage in our context would be limited to safe nets.

For the general case, the current state of the art can be found in the work of De Loera et al. [START_REF] Jesús | Algebraic and Geometric Ideas in the Theory of Discrete Optimization[END_REF][START_REF] Jesús | Effective lattice point counting in rational convex polytopes[END_REF] on counting lattice points in convex polytopes. Their approach is implemented in a tool called LaTTe; it relies on algebraic and geometric methods; namely the use of rational functions and the decomposition of cones into unimodular cones. Our experiments with LaTTe show that it can be conveniently used on systems with, say, less than 50 variables. For instance, LaTTe is strikingly fast (less than 1s) at counting the number of solutions of the system computed in Sect. 4. Moreover, its running time does not generally depend on the constants found in the system. As a consequence, computing the reachability count for 10 or, say, 10 12 houses takes exactly the same time.

An ad-hoc method. Though our experiments with LaTTe suffice to show that these approaches are practicable, we implemented our own counting method. Its main benefits over LaTTe, important for practical purposes, are that it can handle systems with many variables (say thousands), though it can be slower than LaTTe on small systems. Another reason for finding an alternative to LaTTe is that it provides no builtin support for parameterized systems, that is in the situation where we need to count the solutions of many instances of the same linear system differing only by some constants.

Our solution takes advantage of the stratified structure of the systems obtained from reductions, and it relies on combinatorial rather than geometric methods. While we cannot describe this tool in full details, we illustrate our approach and the techniques involved on a simple example.

Consider the system of equations obtained from the PN corresponding to the dashed zone of Fig. 4. This net consists of places in the range p 18 -p 27 and is reduced by our system to a single place, a 13 . The subset of marking equations related to this subnet is: Assume place a 13 is marked with n tokens. Then, by Th. 5.1, the number of markings of the original net corresponding with marking a 13 = n in the reduced net is the number of non-negative integer solutions to system (Q, a 13 = n). Let us define the function A 13 : N -→ N that computes that number. We first simplify system (Q). Note that no agglomeration is involved in the redundancy (R) equations for p 19 and p 27 , so these equations have no effects on the marking count and can be omitted. After elimination of variable a 5 and some rewriting, we obtain the simplified system (Q ): . Now consider the system consisting of the first two equations and the redundancy equation a 4 = a 11 . If a 11 = x, its number of solutions is ((2))(x) × ((3))(x) (the variables in both equations being disjoint). Finally, by noticing that a 11 can take any value between 0 and n, we get:

A 13 (n) = n a 11 =0 ((2))(a 11 ) × ((3))(a 11 )
This expression is actually a polynomial, namely 1 8 n 4 + 11 12 n 3 + 19 8 n 2 + 31 12 n + 1. By applying the same method to the whole net of Fig. 4, we obtain an expression involving six summations, which can be reduced to the following 18th-degree polynomial in the variable X denoting the number of tokens in place p 1 . In the general case of partially reduced nets, the computed polynomial is a multivariate polynomial with at most as many variables as places remaining in the residual net. When that number of variables is too large, the computation of the final polynomial is out of reach, and we only make use of the intermediate algebraic term.

Computing Experiments

We integrated our reduction system and counting method with a state space generation tool, tedd, in the framework of our TINA toolbox for analysis of Petri nets [START_REF] Berthomieu | The tool TINAconstruction of abstract state spaces for Petri nets and Time Petri nets[END_REF] (www. laas.fr/tina). Tool tedd makes use of symbolic exploration and stores markings in a Set Decision Diagram [START_REF] Thierry-Mieg | Hierarchical set decision diagrams and regular models[END_REF]. For counting markings in presence of agglomerations, one has the choice between using the external tool LaTTe or using our native counting method discussed in Sect. 6.

Benchmarks. Our benchmark is constituted of the full collection of Petri nets used in the Model Checking Contest [START_REF] Kordon | Complete Results for the 2017 Edition of the Model Checking Contest[END_REF][START_REF] Hillah | Petri Nets Repository: a tool to benchmark and debug Petri Net tools[END_REF]. It includes 627 nets, organized into 82 classes (simply called models). Each class includes several nets (called instances) that typically differ by their initial marking or by the number of components constituting the net. The size of the nets vary widely, from 9 to 50 000 places, 7 to 200 000 transitions, and 20 to 1 000 000 arcs. Most nets are ordinary (arcs have weight 1) but a significant number are generalized nets. Overall, the collection provides a large number of PN with various structural and behavioral characteristics, covering a large variety of use cases. Reduction ratio and prevalence. Our first results are about how well the reductions perform. We provide two different reduction strategies: compact, that applies all reductions described in Sect. 3, and clean, that only applies removal of redundant places and transitions. The reduction ratios on number of places (number of places before and after reduction) for all the MCC instances are shown in Fig. 7, sorted in descending order. We overlay the results for our two reduction strategies (the lower, in light color, for clean and the upper, in dark, for compact). We see that the impact of strategy clean alone is minor compared to compact. Globally, Fig. 7 shows that reductions have a significant impact on about half the models, with a very high impact on about a quarter of them. In particular, there is a surprisingly high number of models that are totally reducible by our approach (about 19% of the models are fully reducible).

Computing time of reductions. Many of the reduction rules implemented have a cost polynomial in the size of the net. The rule removing redundant places in the general case is more complex as it requires to solve an integer programming problem. For this reason we limit its application to nets with less than 50 places. With this restriction, reductions are computed in a few seconds in most cases, and in about 3 minutes for the largest nets. The restriction is necessary but, because of it, we do not reduce some nets that would be fully reducible otherwise.

Impact on the marking count problem. In our benchmark, there are 169 models, out of 627, for which no tool was ever able to compute a marking count. With our method, we could count the markings of at least 14 of them.

If we concentrate on tractable nets-instances managed by at least one tool in the MCC 2017-our approach yields generally large improvements on the time taken to count markings; sometimes orders of magnitude faster. Table 1 (top) lists the CPU time (in seconds) for counting the markings on a selection of fully reducible instances. We give the best time obtained by a tool during the last MCC (third column) and compare it with the time obtained with tedd, using two different ways of counting solutions (first with our own, native, method then with LaTTe). We also give the resulting speed-up. These times also include parsing and applying reductions. An absent value (-) means that it cannot be computed in less than 1 hour with 16 Gb of storage.

Concerning partially reducible nets, the improvements are less spectacular in general

Related Work and Conclusion

Our work relies on well understood structural reduction methods, adapted here for the purpose of abstracting the state space of a net. This is done by representing the effects of reductions by a system of linear equations. To the best of our knowledge, reductions have never been used for that purpose before. Linear algebraic techniques are widely used in Petri net theory but, again, not with our exact goals. It is well known, for instance, that the state space of a net is included in the solution set of its so-called "state equation", or from a basis of marking invariants. But these solutions, though exact in the special case of live marked graphs, yield approximations that are too coarse. Other works take advantage of marking invariants obtained from semiflows on places, but typically for optimizing the representation of markings in explicit or symbolic enumeration methods rather than for helping their enumeration, see e.g. [START_REF] Schmidt | Using Petri net invariants in state space construction[END_REF][START_REF] Wolf | Generating Petri net state spaces[END_REF]. Finally, these methods are only remotely related to our.

Another set of related work concerns symbolic methods based on the use of decision diagrams. Indeed they can be used to compute the state space size. In such methods, markings are computed symbolically and represented by the paths of some directed acyclic graph, which can be counted efficiently. Crucial for the applicability of these methods is determining a "good" variable ordering for decision diagram variables, one that maximizes sharing among the paths. Unfortunately, finding a convenient variable ordering may be an issue, and some models are inherently without sharing. For example, the best symbolic tools participating to the MCC can solve our illustrative example only for p 1 ≤ 100, at a high cost, while we compute the result in a fraction of a second for virtually any possible initial marking of p 1 .

Finally, though not aimed at counting markings nor relying on reductions, the work reported in [START_REF] Stahl | Decomposing Petri net State Spaces[END_REF] is certainly the closest to our. It defines a method for decomposing the state space of a net into the product of "independent sets of submarkings". The ideas discussed in the paper resemble what we achieved with agglomeration. In fact, the running example in [START_REF] Stahl | Decomposing Petri net State Spaces[END_REF], reproduced here in Figure 1, is a fully reducible net in our approach. But no effective methods are proposed to compute decompositions.

Concluding remarks. We propose a new symbolic approach for representing the state space of a PN relying on systems of linear equations. Our results show that the method is almost always effective at reducing computing times and memory consumption for counting markings. Even more interesting is that our methods can be used together with traditional explicit and symbolic enumeration methods, as well as with other abstraction techniques like symmetry reductions for example. They can also help for other problems, like reachability analysis.

There are many opportunities for further research. For the close future, we are investigating richer sets of reductions for counting markings and application of the method to count not only the markings, but also the number of transitions of the reachability graph. Model-checking of linear reachability properties is another obvious prospective application of our methods. On the long term, a question to be investigated is how to obtain efficiently fully equational descriptions of the state spaces of bounded Petri nets.

Theorem 3 . 1 .Figure 2 :

 312 Figure 2: Some examples of redundant transitions (left) and places (right)

Figure 3 :Definition 3 . 4 (

 334 Figure 3: Agglomeration examples: chain (top), loop (for n = 3, bottom)

Figure 4 :

 4 Figure 4: HouseConstruction-10 example net .

Figure 5 :

 5 Figure 5: Reduction traces for net HouseConstruction-10

Figure 6 :

 6 Figure 6: Non confluence example

R p 19

 19 = p 20 R p 27 = p 23 + p 26 A a 4 = p 21 + p 18 A a 11 = p 23 + a 5 A a 5 = p 22 + p 20 A a 13 = a 11 + p 26 R a 4 = a 5 + p 23 (Q)

A a 4 =

 4 p 21 + p 18 R a 4 = a 11 A a 11 = p 23 + p 22 + p 20 A a 13 = a 11 + p 26 (Q ) Let ((k))(x) denote the expression x+k-1 k-1 , which denotes the number of ways to put x tokens into k slots. The first equation is a 4 = p 21 + p 18 . If a 4 = x, its number of solutions is ((2))(x) = x + 1. The second equation is a 11 = p 23 + p 22 + p 20 . If a 11 = x, its number of solutions is ((3))(x) = (x+2)(x+1) 2

Figure 7 :

 7 Figure 7: Distribution of reduction ratios (place count) over the 627 PN instances.

  . A transition t in T

	p 0	p 3	t 2	p 4	p 7	t 5	p 8
	t 0
							p
			t 1	t 3		t 4	t 6
		p 1		p 2	p 5		p 6

9 Figure 1: An example Petri net is said enabled at m if m ≥ Pre(t). If enabled at m, transition t may fire yielding a marking m = m -Pre(t) + Post(t). This is written m t →m , or simply m → m when only markings are of interest. Intuitively, places hold integers and together encode the state (or marking) of a net; transitions define state changes. The reachability set, or state space, of N is the set of markings R(N ) = { m | m 0 * → m }, where * → is the reflexive and transitive closure of →.

* This material is based upon work supported by the RTRA STAE project IFSE2: "technical engineering for embedded systems".

net-abstraction relation, defined below.

We start by defining some notations useful in our proofs. We use U, V, . . . for finite sets of non-negative integer variables. We use Q, Q for systems of linear equations (and inequalities) and the notation V(Q) for the set of variables occurring in Q. The system obtained by concatenating the relations in Q 1 and Q 2 is denoted (Q 1 ; Q 2 ) and the "empty system" is denoted ∅.

A valuation e of N V is a solution of Q, with V = V(Q), if all the relations in Q are (trivially) valid when replacing all variables x in V by their value e(x). We denote Q the subset of N V(Q) consisting in all the solutions of Q.

If E ⊆ N V then E ↓ U is the projection of E over variables U, that is the subset of N U obtained from E by restricting the domain of its elements to U. conversely, we use E ↑ U to denote the lifting of E to U, that is the largest subset E of N U such that E ↓ V = E. Definition 5.1 (Net-abstraction). A triple (N 1 , Q, N 2 ) is a net-abstraction, or simply an abstraction, if N 1 , N 2 are nets with respective sets of places P 1 , P 2 (we may have P 1 ∩ P 2 = ∅), Q is a linear system of equations, and:

Intuitively, N 2 is an abstraction of N 1 (through Q) if, from every reachable marking m ∈ R(N 2 ), the markings obtained from solutions of Q-restricted to those solutions such that x = m(x) for all "place variable" x in P 2 -are always reachable in N 1 . The definition also entails that all the markings in R(N 1 ) can be obtained this way. Proof. Property (1) is obvious from Definition 5.1. Property ( 2) is proved by case analysis. First, let V = V(Q) ∪ P 1 ∪ P 2 and U = V ∪ P 3 and notice that for all candidate

Case (T) : Q = Q. By Th. 3.1, we have 1: Computation times (in seconds) and speed-up for counting markings on some totally (top) and partially (bottom) reduced nets though still significant. Counting markings in this case is more expensive than for totally reduced nets. But, more importantly, we have to build in that case a representation of the state space of the residual net, which is typically much more expensive than counting markings. Furthermore, if using symbolic methods for that purpose, several other parameters come into play that may impact the results, like the choice of an order on decision diagram variables or the particular kind of diagrams used. Nevertheless, improvements are clearly visible on a number of example models; some speedups are shown in Table 1 (bottom). Also, to minimize such side issues, instead of comparing tedd with compact reductions with the best tool performing at the MCC, we compared it with tedd without reductions or with the weaker clean strategy. In that case, compact reductions are almost always effective at reducing computing times. Finally, there are also a few cases where applying reductions lower performances, typically when the reduction ratio is very small. For such quasi-irreducible nets, the time spent computing reductions is obviously wasted.