Romain Ravaud

Guy Lemarquand
email: guy.lemarquand@univ-lemans.fr

COMPARISON OF THE COULOMBIAN AND AMPERIAN CURRENT MODELS FOR CALCULATING THE MAGNETIC FIELD PRODUCED BY RADIALLY MAGNETIZED ARC-SHAPED PERMANENT MAGNETS

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

calculating the three magnetic field components produced by tile permanent magnets or ring permanent magnets radially magnetized. These four methods derived directly from the Maxwell's equations and the properties of the magnetic field shape. Consequently, some authors used the coulombian model or the magnetic potential [START_REF] Babic | Calculation improvement of 3D linear magnetostatic field based on fictitious magnetic surface charge[END_REF][START_REF] Babic | An improvement in the calculation of the magnetic field for an arbitrary geometry coil with rectangular cross section[END_REF][START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF][START_REF] Ravaud | Permanent magnet couplings: Field and torque threedimensional expressions based on the coulombian model[END_REF] while other authors used the amperian current model and the Biot-Savart Law [START_REF] Azzerboni | Magnetic field evaluation for disk conductors[END_REF][START_REF] Azzerboni | Threedimensional calculation of the magnetic field created by currentcarrying massive disks[END_REF][START_REF] Babic | New expressions for calculating the magnetic field created by radial current in massive disks[END_REF][START_REF] Akyel | Mutual inductance calculation for non-coaxial circular air coils with parallel axes[END_REF] for calculating the magnetic field created by currents or permanent magnets. However, these analytical methods are not suitable for the study of iron-core structures. This is why some alternative methods are based on series expansions [START_REF] Jian | Analytical calculation of magnetic field distribution in coaxial magnetic gears[END_REF], or finite-element methods are used for the design of permanent-magnet topologies [START_REF] Chau | Transient analysis of coaxial magnetic gears using finite element comodeling[END_REF].

The coulombian model implies the existence of magnetic charges located on the faces of a magnet and inside it. The amperian current model implies the determination of the vector potential produced by fictitious currents flowing around or inside a magnet. The two approaches are equivalent. The problem is thus to guess what model is the most appropriate for calculating the three components of the magnetic field produced by permanent magnets. It does not seem to be more difficult to use the amperian current model rather than the coulombian model for calculating the magnetic field created by parallelepiped magnets. Indeed, the coulombian model applied to a parallelepiped magnet whose polarization is directed along the normal direction implies two times the calculation of two surface integrals whereas the amperian current model implies four times the calculation of two surface integrals.

However, the two analytical methods lead to a fully analytical expression of the magnetic field [START_REF] Akoun | 3D analytical calculation of the forces exerted between two cuboidal magnets[END_REF]. This is why this configuration is often used by authors for optimization purposes: it has a very low computational cost [START_REF] Elies | Analytical optimization of the torque of a permanent-magnet coaxial synchronous coupling[END_REF][START_REF] Lemarquand | Nonsinusoidal torque of permanent-magnet couplings[END_REF][START_REF] Lemarquand | Annular magnet position sensor[END_REF][START_REF] Blache | New structures for linear displacement sensor with hight magnetic field gradient[END_REF][START_REF] Charpentier | Calculation of ironless permanent magnet coupling using semi-numerical magnetic pole theory method[END_REF]. However, as shown in [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF], this method can only be used in the near-field if parallelepiped magnets are used for representing arc-shaped permanent magnets. It is emphasized here that numerous approaches based on series expansions have been proposed [START_REF] Selvaggi | Computation of the three-dimensional magnetic field from solid permanent-magnet bipolar cylinders by employing toroidal harmonics[END_REF][START_REF] Selvaggi | Computation of the external magnetic field, near-field or far-field from a circular cylindrical magnetic source using toroidal functions[END_REF][START_REF] Zhilichev | Calculation of magnetic field of tubular permanent magnet assemblies in cylindrical bipolar coordinates[END_REF] while semi-analytical approaches have led to compute the magnetic field in all points in space [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF]. Moreover, other analytical methods allow the determination of the force between coil currents [START_REF] Conway | Noncoaxial inductance calculations without the vector potential for axisymmetric coils and planar coils[END_REF][START_REF] Kim | Mutual inductance of noncoaxial circular coils with constant current density[END_REF].

For arc-shaped permanent magnets, it seems to be more difficult to guess what model is the most appropriate. For arc-shaped permanent magnets whose polarization is radial, some new analytical expressions are derived from the vector potential, not only for ring but also for tile permanent magnets. These expressions are compared with the ones established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF], Babic [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF], Rakotoarison [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF] and Ravaud [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. Moreover, we discuss the validity and the computational cost of all these expressions. A Mathematica file is available online for showing the validity of the expressions obtained and the equivalence between the amperian current and couloumbian models [33].

The first part of this paper deals with the calculation of the magnetic field created by ring permanent magnets radially magnetized, and the second part of this paper deals with the magnetic field produced by tile permanent magnets radially magnetized. By using the amperian current model, we show how some expressions can be significantly improved and why some expressions cannot still be expressed entirely in terms of special functions.

RING PERMANENT MAGNETS WHOSE POLARIZATION IS RADIAL

Notation and Geometry

In this first part, we study the magnetic field produced by ring permanent magnets whose polarization is radial. For this purpose, let us consider the geometry shown in Figure 1. Figure 1(a) represents a ring permanent magnet radially magnetized according to the coulombian model. The ring inner radius is r 1 , and its outer one is r 2 . Its height is h = z 2 -z 1 . Its magnetic polarization J is directed towards 0. Some fictitious charge surface densities are located on the inner and outer faces of the ring, and a charge volume density appears inside the ring. Consequently, as mentioned in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF], up to three integrals must be determined for obtaining an analytical expression of the three magnetic field components created by this ring permanent magnet. The contribution of the charge surface densities can be expressed in terms of elliptic integrals of the first, second and third kind [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF] but the contribution of the charge volume density can be expressed in an analytical part based on elliptic functions and a numerical part [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. Consequently, the coulombian model applied to this configuration does not seem to be the best model for obtaining a fully analytical expression of the magnetic field.

In this paper, we use the amperian current model for calculating the magnetic field produced by a ring permanent magnet radially magnetized, as shown in Figure 1(b). The magnetic polarization J is transformed into a surface density currents k flowing on the upper and lower faces of the ring. Such a model has been employed for example by Lang [START_REF] Lang | Fast calculation method for the forces and stiffnesses of permanent-magnet bearings[END_REF] in the case of axially magnetized ring permanent magnets. The ratio between J and k is expressed as follows:

k = J µ 0 (1)
By using the amperian current model in the case of radially magnetized ring permanent magnets, the expression of the radial component can be improved significantly because it can be expressed only in terms of elliptic integrals, and no further numerical integrations are required. Consequently, the time necessary to calculate this component is lower with the amperian current model than the coulombian model. However, the two models are equivalent and give the same results as it is shown in the next section. We present now the new magnetic field expressions determined with the amperian current model.

The analytical expression of the radial component can be determined by using either the vector potential or directly the Biot-Savart law. Indeed, by using the analogy between the coulombian model and amperian current model, the expression of the radial component can be expressed in terms of elliptic integrals. By denoting r the observation point and r a point located on the charge distribution on the magnet, it is useful to define the reciprocal distance r -r as follows:

1

r -r = 1 r 2 + r2 -2rr cos( θ) + (z -z) 2 (2)
By using the amperian current model and as stated in [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF], the potential vector created by a ring permanent magnet whose polarization is radial can be expressed as follows:

A( r) = µ 0 4π   S up k( r ) r -r d Sup - S down k( r ) r -r d Sdown   (3)
where S up is the upper surface of the ring (z = z 2 ), and S down is the lower surface of the ring (z = z 1 ). Moreover, it is noted that in this analogy, the current volume density k v is 0 because the polarization is supposed to be perfectly radial.

k v = ∇ ∧ J µ 0 = 0 (4)
Thus, the magnetic field can be obtained by the following expressions:

H( r) = 1 µ 0 ∇ ∧ A( r) (5)
As there is a symmetry according θ in our configuration, the azimuthal component H θ (r, z) equals zero, and the magnetic field depends only on r and z, that is:

H( r) = H(r, z) = H r (r, z) u r + H z (r, z) u z (6)
For the rest of this paper, we use the operator

2 i=1 2 k=1 (-1) i+k [•]
for the case of ring permanent magnets whose polarization is radial. However, according to the convention chosen for calculating the three magnetic components, this operator can be replaced by

2 i=1 2 k=1 (-1) i+k+1 [•]
. The modulus of the magnetic field created is always the same, but the sign of the magnetic components depends directly on the way that we define the axes. Such a way of writing the magnetic field expressions has also been adopted by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] and Babic [START_REF] Babic | Improvement in the analytical calculation of the magnetic field produced by permanent magnet rings[END_REF]. The only influence of the choice of (-1)

i+k [•], (-1) i+k+1 [•] or (-1) i+j+k [•] and (-1) i+j+k+1 [•]
is the change in sign of the magnetic field component expressions.

Radial Component Determined with the Amperian Current Model

The radial component H r (r, z) can be determined by calculating the projection of H( r) along u r :

H r (r, z) = H( r) • u r = 1 µ 0 ∇ ∧ A( r) • u r (7)
The interest of using the amperian current model for this configuration lies in the fact that only two integrals must be calculated. Consequently, we obtain an analytical expression of the radial field component based on elliptic integrals. By using Mathematica, the arguments of the elliptic integrals used are defined as follows:

φ +,- 1 = (b + 2e)x bx ± √ 2 xe 2 (x -c) φ 2 = i sinh -1 -1 b + 2e b -2e cos( θ) φ 3 = b + 2e b -2e (8)
It is noted that we use mainly the elliptic integrals of the second and third kind for the calculation of the radial component. These special functions have also been used in [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF][START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF]. However, no further numerical integrations are required here. In short, the radial component H r (r, z) can be expressed as follows:

H r (r, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (g(i, k, 2π) -g(i, k, 0)) (9)
with

g i, k, θ = 2(z-z k )f r 2 +(z-z k ) 2 , r 2 +r 2 i +(z-z k ) 2 , rr i , -r 2 -2(z-z k ) 2 , θ (10)
where

f a, b, e, c, x, θ = η 2ξ 1 2ce 2 + ξ 2 F * [φ 2 , φ 3 ] +η -e 2 (c -x) bx √ 2 + 2ξ 1 Π * φ + 1 , φ 2 , φ 3 +η e 2 (c -x) bx √ 2 -2ξ 1 Π * φ - 1 , φ 2 , φ 3 -2ηax xe 2 -ce 2 √ 2 + bξ 1 Π * φ + 1 , φ 2 , φ 3 -2ηax -xe 2 +ce 2 √ 2+bξ 1 Π * φ - 1 , φ 2 , φ 3 (11)
where F * [x, y] and Π * [x, y, z] are the elliptic integrals of the second and third kind that have been used in previous papers [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF][START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF].

In addition, the parameters ξ 1 , ξ 2 , η are defined as follows:

ξ 1 = e 2 x(x -c) ξ 2 = x(b 2 -2e 2 ) η = i -e 2 sin( θ) 2 (b-2e) 2 csc θ 2 -1 b+2e xξ 1 (2ce 2 + ξ 2 ) (12)
Therefore, by using the amperian current model applied to a ring permanent magnet radially magnetized, it is possible to obtain a fully analytical expression of the radial field component. This result is of great importance because it can be used for optimizing the radial field in ironless loudspeakers [START_REF] Lemarquand | Ironless loudspeakers[END_REF][START_REF] Ravaud | Ironless loudspeakers with ferrofluid seals[END_REF].

Illustration of the Analogy between the Amperian Current and Coulombian Models for the Radial Component Expression H r (r, z)

A comparison has been performed for verifying the accuracy of the radial component expression with the expressions determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] and [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. It is noted that this expression can also be compared to the expression determined by Furlani by using θ s (2) -θ s (1) = 2π and by neglecting the current densities flowing in the axial direction (see [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]). We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, z = 0.0035 m. We represent in Figure 2 the radial component H r (r, z = 0.0035 m) versus the axial displacement z.

Figure 2 clearly shows that the analytical expression determined in this paper is accurate and can be used instead of the one established in [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF] because it is easier to use for optimization purposes. In fact, we can see that the choice of the model used is very important because it seems to be more difficult to obtain a fully analytical expression of the radial component H r (r, z) in the coulombian approach. On the other hand, this does not imply that the amperian current model is always more appropriate than the coulombian model, as shown in [START_REF] Ravaud | Analytical expression of the magnetic field created by tile permanent magnets tangentially magnetized and radial currents in massive disks[END_REF] in which radial currents in massive disks generate a magnetic field that can be deducted more easily from the coulombian model. In fact, more generally, for cylindrical geometries, the authors feel that the analysis of the geometry studied must receive particular attention according to the choice of the model used.

Axial Component H z (r, z)

The axial component H z (r, z) can be determined by calculating the projection of H( r) along u z :

H z (r, z) = H ( r ) • u z = 1 µ 0 ∇ ∧ A( r) • u z ( 13
)
This axial component H z (r, z) is thus given by:

H z (r, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (I 1 (r, z) + I 2 (r, z)) ( 14
)
with

I 1 (r, z) = -4r i (r -r i ) 2 + (z -z k ) 2 K * -4rr i (r -r i ) 2 + (z -z k ) 2 (15)
where K * [x] is the complete elliptic integral of the first kind [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF].

I 2 (r, z) = 2π 0 log r i -r cos( θ)+ r 2 +r 2 i + (z -z k ) 2 -2rr i cos( θ) d θ (16)
It is emphasized here that this expression is given in a more compact form than the one determined in a previous paper by the authors [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. Consequently, it seems to be more judicious to use the amperian current model for calculating the axial component of the magnetic field produced by a ring permanent magnet radially magnetized rather than the coulombian model. However, there is still a part of this expression that can be determined by using a numerical integration.

Illustration of the Analogy between the Amperian Current and Coulombian Models for the Radial Component Expression H z (r, z)

A comparison has been performed for verifying the accuracy of the axial component expression with the expressions determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] and [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. It is noted that this expression can also be compared to the expression determined by Furlani by using θ s (2) -θ s (1) = 2π. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.024 m. We represent in Figure 3 the axial component H z (r = 0.024 m, z) versus the axial displacement z. Figure 3 clearly shows that the analytical expression determined in this paper is accurate and can be used instead of the one established in [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF] because it has a lower computational cost, and it is given in a more compact form.

TILE PERMANENT MAGNETS WHOSE POLARIZATION IS RADIAL

Notation and Geometry

In this second part, we study the magnetic field produced by tile permanent magnets whose polarization is radial. For this purpose, let us consider the geometry shown in Figure 4. Figure 4(a) represents a tile permanent magnet radially magnetized according to the coulombian model. The ring inner radius is r 1 , and its outer one is r 2 . Its height is h = z 2 -z 1 . Its magnetic polarization J is directed towards 0. Its angular width is θ 2 -θ 1 . Some fictitious charge surface densities are located on the inner and outer faces of the ring, and a charge volume density appears inside the ring. Figure 4(b) represents a tile permanent magnet radially magnetized according to the amperian current model. As stated in the first section of this paper, we consider the current surface densities flowing around the tile permanent magnet.

In addition, we define the following operator for the rest of this paper:

℘ i,j,k [•] = J 4πµ 0 2 i=1 2 j=1 2 k=1 (-1) i+j+k [•] (17) u z u r h 1 2 1 r r 2 z 2 z 1 J u k u z u r h 1 2 1 r r 2 z 2 z 1 k u (a) (b)

Radial Component H r (r, θ, z)

We use the same equation as Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] for calculating the radial component of the magnetic field produced by a tile permanent magnet whose polarization is radial. [START_REF] Akyel | Mutual inductance calculation for non-coaxial circular air coils with parallel axes[END_REF] where 1 r -r (α, β)

H r (r, θ, z) = J 4πµ 0 z 2 z 1 r 2 r 1 r sin (θ -θ 1 ) r -r (θ 1 , z) 3 drdz - J 4πµ 0 z 2 z 1 r 2 r 1 r sin (θ -θ 2 ) r -r (θ 2 , z) 3 drdz + J 4πµ 0 r 2 r 1 θ 2 θ 1 (z -z 2 ) cos θ -θ r -r θ, z 2 3 rd θdr - J 4πµ 0 r 2 r 1 θ 2 θ 1 (z -z 1 ) cos θ -θ r -r θ, z 1 3 rd θdr
= 1 r 2 + r2 -2rr cos(θ -α) + (z -β) 2 (19)
It is noted that Eq. ( 18) can be decomposed in two parts. The contribution of the charges located on the curved planes (inner and outer faces of the ring) can be expressed in terms of elliptic integrals by using (9) in which

2 i=1 2 k=1 (-1) i+k (g(i, k, 2π) -g(i, k, 0)) is transformed into 2 i=1 2 k=1 (-1) i+k (g(i, k, θ -θ 2 ) -g(i, k, θ -θ 1 )
). Thus, the radial field component H r (r, θ, z) can be decomposed as follows:

H r (r, θ, z) = J 4πµ 0 2 i=1 2 k=1 (-1) i+k (g(i, k, θ -θ 2 ) -g(i, k, θ -θ 1 )) + J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 1 ) r -r 3 drdz - J 4πµ 0 z 2 z 1 r 2 r 1 r sin(θ -θ 2 ) r -r 3 drdz (20)
The last r and z integrations can be performed analytically, and we obtain the following expression for the radial field component:

H r (r, θ, z) = ℘ i,j,k [g(i, k, θ -θ j ) + k i,j,k + l i,j,k ] (21) k i,j,k = sin(θ -θ j ) log [Y i,j,k ] (22) l i,j,k = -cos(θ -θ j ) arctan [X i,j,k ] (23)
Y i,j,k = (z -z k ) + r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) X i,j,k = (z -z k )(r i -r cos(θ -θ j )) r sin(θ -θ j ) r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) (24)
As a conclusion, it is possible to obtain an analytical expression of the radial field component produced by a tile permanent magnet radially magnetized. This result is interesting for the optimization of radial fields in electrical machines.

Illustration of the Analogy between the Amperian Current and Coulombian Models for the Radial Component Expression H r (r, z) in the Case of a Tile Permanent Magnet

The expression of the radial field component produced by a tile permanent magnet radially magnetized is compared with the expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.02 m. We represent in Figure 5 H r (r = 0.024 m, z) versus the angular displacement θ. Representation of the radial component H r (r = 0.02 m, θ, z = 0.0035 m)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.02 m; (thick line = this work), (Points = Furlani's work [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]).

Azimuthal Component H θ (r, θ, z)

Let us now consider the azimuthal component H θ (r, θ, z) of the magnetic field produced by a tile permanent magnet radially magnetized. We use the same equation as Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF] for calculating this component. Therefore, by using the amperian current model, the azimuthal component is expressed as follows:

H θ (r, θ, z) = J 4πµ 0 2 j=1 (-1) j    r 2 r 1 θ 2 θ 1 (z -z) sin θ -θ r -r j r, θ 3    rdrd θ + J 4πµ 0 2 j=1 (-1) j    r 2 r 1 z 2 z 1 r -r cos(θ -θ j ) r -r (r, z) 3    drdz (25)
where 1

r -r j (α, β) = 1 r 2 + α 2 -2rα cos(θ -β) + (z -z j ) 2 1 r -r (α, β) = 1 r 2 + α 2 -2rα cos(θ -θ j ) + (z -β) 2 (26)
We find:

H θ (r, θ, z) = ℘ i,j,k [-cos(θ -θ j ) log [Y i,j,k ] -sin(θ -θ j ) arctan [X i,j,k ]] +℘ i,j,k (z -z k ) r (r i -r cos(θ -θ j ) + ξ i,j,k ) (27)
where

ξ i,j,k = r 2 + r 2 i + (z -z k ) 2 -2rr i cos(θ -θ j ) (28)
Therefore, it is possible to obtain a fully analytical expression of the azimuthal field component by using the amperian current model.

Illustration of the Analogy between the Amperian Current and Coulombian Models for the Azimuthal

Component Expression H θ (r, z) in the Case of a Tile Permanent Magnet

The expression of the azimuthal field component produced by a tile permanent magnet radially magnetized is compared with the expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions:

r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, -3 -2 -1 0 1 2 3 theta [rad]
r = 0.024 m, θ 2 -θ 1 = π 4 .
We represent in Figure 7 the axial component H z (r = 0.024, θ, z = 0.002) versus the angular displacement θ.

Figure 6 shows that the fully analytical expression obtained in this work for the azimuthal component gives the same azimuthal field as the one established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. However, the expression we give in this paper does not use any special functions or numerical integrals. Consequently, such an expression has a lower computational cost and is useful for optimization purposes.

Axial Component H z (r, θ, z)

The axial component H z (r, θ, z) is also determined with the same integral formulation of Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF].

H z (r, θ, z) = J 4πµ 0 2 j=1 (-1) j r 2 r 1 θ 2 θ 1 rr cos θ -θ r -r (α, γ) 3 drd θ - J 4πµ 0 2 j=1 (-1) j r 2 r 1 θ 2 θ 1 r2 r -r (α, γ) 3 drd θ (29) with 1 r -r (α, γ) = 1 r 2 + γ 2 -2rγ cos(θ -α) + (z -z j ) 2 (30)
After the two integrations, we obtain a compact form of the axial component:

H z (r, θ, z) = ℘ i,j,k -2r i (r -r i ) 2 + (z -z k ) 2 F * θ -θ j 2 , - 4rr i (r -r i ) 2 + (z -z k ) 2 + J 4πµ 0 θ 2 θ 1 log r i -r cos θ -θ + ξ θ d θ (31)
with

ξ θ = r 2 + r 2 i + (z -z k ) 2 -2rr i cos θ -θ (32)
where F * [φ, m] gives the elliptic integral of the first kind [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]. This expression has also been improved compared to the one determined in [START_REF] Ravaud | The three exact components of the magnetic field created by a radially magnetized tile permanent magnet[END_REF] because it is presented in a more compact form.

Illustration of the Analogy between the Amperian Current and Coulombian Models for the Axial Component

Expression H z (r, z) in the Case of a Tile Permanent Magnet

The expression of the axial field component produced by a tile permanent magnet radially magnetized is compared with the expression obtained by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. We take the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.015 m, θ 2 -θ 1 = π 4 . We represent in Figure 7 the axial component H z (r = 0.015, z = 0.002 m) versus the angular displacement θ. Figure 7 shows that the reduced expression obtained in this work for the axial component gives the same axial field as the one established by Furlani [START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF].

CONCLUSION

This paper has presented some improved analytical expressions of the magnetic field produced by ring and tile permanent magnets radially magnetized. The radial field is given in terms of elliptic integrals, and no further numerical integrations are required. Consequently, the time necessary for calculating this radial field is very low and allows us to easily carry out parametric studies. In the case of tile permanent magnets radially magnetized, the amperian current model leads to a fully analytical expression of the azimuthal component, as in the case of the coulombian model. However, only two integrals must be determined with the amperian current model, and we think that this model is more appropriate than the coulombian model for the study of this azimuthal field. The axial field produced by a tile permanent magnet or a ring permanent magnet can be expressed as in an analytical part based on elliptic integrals and a numerical part. This expression has been improved compared to the one established by Furlani, but we have not succeeded in analytically integrating the last term. Nevertheless, we have presented the expression of the axial field in a compact form, and this expression seems to be easier to use than the one determined by the authors in a previous paper with the coulombian model [START_REF] Ravaud | Discussion about the analytical calculation of the magnetic field created by permanent magnets[END_REF]. The amperian current model seems to be more appropriate than the coulombian model for calculating the magnetic field produced by radially magnetized arc-shaped permanent magnets.

Figure 1 .

 1 Figure 1. Representation of a ring permanent magnet radially magnetized: (a) coulombian model, (b) amperian current model.

Figure 2 .

 2 Figure 2. Representation of the radial component H r (r, z = 0.0035) versus the radial displacement r with the following values: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T; (thick line = amperian current model), (points = coulombian model).

Figure 3 .

 3 Figure 3. Representation of the axial component H z (r = 0.0249 m, z) versus the axial displacement z with the following values: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J 1 T; (thick line = amperian current model), (points = coulombian model).

Figure 4 .

 4 Figure 4. Representation of a tile permanent magnet radially magnetized: (a) coulombian model, (b) amperian current model.

Figure 5 .

 5 Figure 5.Representation of the radial component H r (r = 0.02 m, θ, z = 0.0035 m)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.02 m; (thick line = this work), (Points = Furlani's work[START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]).

Figure 6 .

 6 Figure 6. Representation of the azimuthal component H θ (r = 0.024, θ, z = 0.002)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.024 m, θ 2 -θ 1 = π 4 rad; (thick line = this work), (Points = Furlani's work [1]).

Figure 7 .

 7 Figure 7. Representation of the axial component H z (r = 0.015, θ, z = 0.002 m)) versus the angular displacement θ with the following dimensions: r 1 = 0.025 m, r 2 = 0.028 m, h = z 2 -z 1 = 0.003 m, J = 1 T, r = 0.015 m, θ 2 -θ 1 = π 4 ;(thick line = this work), (Points = Furlani's work[START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]).


Marie Badaroux
email: marie.badaroux@univ-grenoble-alpes.fr

Julie Dumas
email: julie.dumas@univ-grenoble-alpes.fr

Frédéric Pétrot
email: frederic.petrot@univ-grenoble-alpes.fr

Fast Instruction Cache Simulation is Trickier than You Think

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

INTRODUCTION

Hardware/Software co-simulation is a technology that is very useful during the design, implementation and even actual use of processor centric systems. It can be implemented in many different ways [START_REF] Jebali | Decoupling Processor and Memory Hierarchy Simulators for Efficient Design Space Exploration[END_REF][START_REF] Pétrot | On MPSoC Software Execution at the Transaction Level[END_REF], as does the simulator responsible for the execution of the software parts of the system. For this latter task, Dynamic Binary Translation (DBT) has proven to be very efficient and scalable. Its simulation strategy is based on the following consideration: a sequence of non-branch instructions ended by a branch can be executed atomically. Therefore, by translating a sequence once and reusing it for its many executions, the translation time is amortized, and the speed of simulation dictated mainly by the speed of execution of the translated code.

When it comes to evaluating performances, DBT is challenging. Speed is due to the execution of a lot of translated instructions out of the control of the simulation environment, and adding instrumentation will rapidly incur large slowdowns.

In this work, we take benefit from the block per block translation nature of the DBT to propose a fast yet accurate instruction cache model. We study the sources of inaccuracy of our model and detail the choices we made to mitigate them.

We assess simulation performance in terms of accuracy and simulation speed using the QEMU [START_REF] Bellard | QEMU, a fast and portable dynamic translator[END_REF] dynamic binary translator emulating the RISC-V instruction set architecture (ISA) running the PolyBench/C [START_REF] Pouchet | Polybench/C 4.1[END_REF] programs as bare-metal, and the multi-threaded PARSEC benchmark [START_REF] Bienia | The PARSEC benchmark suite: Characterization and architectural implications[END_REF] in user-mode. We use as baseline for comparison the current instruction level instrumentation available in vanilla QEMU. We show that performing a fair comparison is not simple, and propose a way to circumvent this issue.

The paper is organized as follows. We first outline in Section 2 the DBT process. Section 3 gives an overview of the most relevant related works on instruction cache simulation and instrumentation. Section 4 details the modeling approach we propose and its limitations and Section 5 details the implementation in QEMU. Section 6 presents the environment we put to work to run our experiments, and the results we obtained. Finally, we conclude in Section 7.

DYNAMIC BINARY TRANSLATION PRIMER

Figure 1 outlines the DBT process. In the frontend, the target instructions are fetched one at a time, decoded, translated into microoperations (µops), and sequentially added into a buffer. If the current instruction is a branch, we enter the backend in which the µops in the buffer are translated into optimized host instructions, and put in a Translation Block (or Tb). A prologue is added in front of the instructions for housekeeping, as is an epilogue to give back control to the simulation environment. The identifier of the sequence of code is the value of the program counter (pc) of the first instruction of the sequence. The Tb is then stored into a translation cache, and immediately executed. Once Tb execution finished, the control is returned to the simulator that seeks the next Tb to execute, identified by the address of the instruction that is the target of the branch ending the previous Tb.

The notion of translation block is very similar to the notion of basic-block coined by [START_REF] Backus | The FORTRAN automatic coding system[END_REF] in the static compilation context, but still, it has some differences. Indeed, the dynamic nature of the compilation makes it possible to have several Tb representing the exact same sequence of target instructions, if the sequence is the target of a goto that jumps into a larger sequence. Also, translation occurs within a known processor environment, including in particular the machine state, e.g. kernel or user-mode, page-table settings, etc. This means in particular that the addresses of the target instructions, be they physical or virtual, are statically known during Tb generation. Additionally, Tbs need to be deleted and retranslated under certain circumstances, such as when dynamically changing an instruction (when adding a breakpoint, patching a dynamically called function address, etc), or when the translation cache is full.

Last, but not least, it might happen that a Tb is only partially executed, when an exception is raised upon execution of an inner target instruction. For instance, a memory access might trigger a page-fault that will immediately end the current Tb to fetch the exception handler.

RELATED WORKS

A recent survey on cache simulators is available in [START_REF] Hadi Brais | A Survey of Cache Simulators[END_REF]. As explained in this necessarily non-exhaustive survey, about 30 years ago cache simulators were mainly studied for research purposes and were the alternative for when the corresponding hardware was not available.

The tool representative of this period is Dinero, whose version IV [START_REF] Edler | Dinero IV Trace-Driven Uniprocessor Cache Simulator[END_REF] is the latest. Nowadays, these simulators are also used to help optimize software stacks, including parallel ones.

Instrumentation is in important mechanism to analyse software but it can be challenging to implement. This topic has been addressed by [START_REF] Luk | Pin: building customized program analysis tools with dynamic instrumentation[END_REF][START_REF] Lyu | Dbill: An efficient and retargetable dynamic binary instrumentation framework using llvm backend[END_REF] that achieve good performance. The Dynamic Binary Instrumentation framework Valgrind [START_REF] Nethercote | Valgrind: a framework for heavyweight dynamic binary instrumentation[END_REF] proposes a more complete instrumentation for heavyweight Dynamic Binary Analysis. With the module Cachegrind, the analysis tool Valgrind is able to simulate classical cache hierarchies [START_REF] Weidendorfer | Sequential performance analysis with callgrind and kcachegrind[END_REF]. Cross-instrumentation, although based on the same principles, is less common [START_REF] Cota | Cross-ISA machine instrumentation using fast and scalable dynamic binary translation[END_REF][START_REF] Guillon | Program instrumentation with qemu[END_REF].

Adding representation of new structures in simulators can be challenging as it can lead to a non negligible execution time overhead. More than 20 years ago, an extension to the instruction set simulator SimICS was proposed by [START_REF] Peter | Efficient instruction cache simulation and execution profiling with a threaded-code interpreter[END_REF] to accurately simulate an instruction cache. Even if the resulting simulation was slowed down, it was deemed acceptable given the accuracy gained. In 2010, [START_REF] Castillo | Fast Instruction Cache Modeling for Approximate Timed HW/SW Co-Simulation[END_REF] presented a high-level instruction cache model. Their approach applied to co-simulation did not degrade simulation speed too much.

Thanks to its execution speed, DBT stands as one of the most used solution for functional high level simulation. However, it lacks microarchitectural details. Works like [START_REF] Böhm | Cycle-accurate performance modelling in an ultra-fast just-in-time dynamic binary translation instruction set simulator[END_REF] produce a DBT-based simulation with addition of architectural details and report a good trade-off between speed and precision.

Works to include cache simulation in QEMU have already been initiated long ago [START_REF] Gligor | Using Binary Translation in Event Driven Simulation for Fast and Flexible MPSoC Simulation[END_REF][START_REF] Van Dung | Cache simulation for instruction set simulator QEMU[END_REF], but were quite intrusive and architecture dependent. [START_REF] Carvalho | mcQEMU: Time-Accurate Simulation of Multi-core platforms using QEMU[END_REF] also extended the dynamic binary translator QEMU and claim time-accurate multiprocessor simulation. With the recent introduction of plugins into QEMU, cache simulation can be done non-intrusively. It is already available upstream [START_REF] Mandour | Cache Modelling TCG Plugin[END_REF] with a simple cache model per virtual CPU (i.e. emulated target CPU or vCPU in short). Our goal is to propose an instruction cache model specific to the DBT mechanism and to show that the overhead it incurs can be significantly reduced.

INSTRUCTION CACHE MODELING

A high-level instruction cache model generally implements only the directory holding the tags and a bit indicating if the tag at a given index is valid or not. The index of the line1 at which an address is stored in the directory, indistinctive of the exact cache geometry, is computed as a combination of the upper address bits (tag) and/or middle address bit (index). The lower bits of the address indicate the exact instruction to fetch within the line, but they are not relevant since the line is either valid as a whole, or invalid as a whole. Given the fact that Tb translation occurs with a known virtual to physical address translation context, virtually addressed, physically addressed, or virtually indexed physically tagged caches can be simulated.

Initial Intuition

Given the Tb per Tb execution principle at work in DBT, we have the warranty that once we have entered a Tb, all subsequent instructions of the Tb are at consecutive (either 2 or 4 bytes away in the RISC-V case) addresses. So, we need to check when entering the Tb whether or not the first instruction misses, as we have no prior knowledge. But from then on, we are sure that all following instructions that share the same index will hit. Then, when the index changes, we might have a miss, but again, for the instructions that follow, we are sure they will hit. So, because addresses are consecutive inside a Tb, we can a priori know which instructions will hit, and which might miss. So we know at Tb creation time that it is sufficient to check for this subset of instructions dynamically. Figure 2 illustrates this principle as an example. Without loss of generality, we assume a 16-byte cache line, which means that the line base addresses have their 4 least significant bits zeroed, and the RISC-V ISA. We jump into this Tb at address 0x800fa7bc, which lays in the middle of a line. We have to check if reading the instruction at this address misses, but we know ahead of time that the next address, 0x800fa7be, will hit. The next instruction is at 0x800fa7c0, which might be a miss as it stands on a new line, so we must check it. However, the 4 following instructions will for sure hit, since they belong to the same line. Overall, on this Tb, we have to run the cache simulator for 2 addresses, while 7 instructions are executed. As a result, we considerably reduce the overhead of doing cache simulation.

Error in Counting Instructions and How to Mitigate it

Unfortunately, the occurrence of exceptions is like a grain of sand seizing up this well oiled mechanism. When we enter a Tb, we make the assumption that all instructions it contains will be executed atomically, and therefore run our model on all addresses that might miss it contains. Alas, when an exception occurs, what is left of the Tb beyond the faulty instruction is not executed, and control is handed over to the simulation environment to fetch the instructions of the exception handler. Then on return from the handler, since the pc of the instruction that was following the access was not known, what was left of the Tb is retranslated into a new smaller Tb, and accounted for by our model. Figure 3 illustrates this with actual page-faults taking place during Linux boot, once the kernel page table has been set-up. We see that there are three re-translations, as the first two Tbs have only been partially executed, and that the two last Tbs are subsets of the first one. In that situation, we ran our cache simulator on 6 non actually executed instructions in the first Tb, and on 3 on the second Tb. Note that the faulty instruction is considered twice, once in the original Tb, and once in the retranslated Tb, but this is expected as the cache might have been trashed by the exception handler. Instructions such as wfi (that sleeps waiting for an interrupt) or pause (that yields back the processor) induce the same behavior, but they are easy to handle. Indeed, they occur so rarely that forcing them to end a Tb induces an insignificant slowdown. On the contrary, memory accesses occur very often and might represent from 30% to 50% of the instruction mix in some actual workloads [START_REF] Faravelon | Removing Load/Store Helpers in Dynamic Binary Translation[END_REF][START_REF] Kurian | Program balance and its impact on high performance RISC architectures[END_REF]. Using the same strategy as previously might incur too much of a slowdown for a negligible accuracy gain, but this needs to be checked, and it will be in Section 6.

Dependency on Simulator Runtime

Another, deeply annoying, unexpected behavior is the dependencies of the flow of executed target instructions on the time it takes to perform cache simulation for programs running on top of Linux 2 . We discovered this seemingly odd behavior while testing several cache implementations: to our amazement, the faster the simulator, the lower the number of executed instructions for a given program. To the best of our knowledge, that effect has not yet been reported in the literature. This makes the evaluations obtained in previous works (we specifically refer here to the ones using QEMU, but we see no reason why similar tools would behave differently), e.g. [START_REF] Gligor | Using Binary Translation in Event Driven Simulation for Fast and Flexible MPSoC Simulation[END_REF][START_REF] Kang | TQSIM: A fast cycleapproximate processor simulator based on QEMU[END_REF][START_REF] Van Dung | Cache simulation for instruction set simulator QEMU[END_REF], questionable, although no wrong per se, since the executed instructions produce a valid behavior. We investigated to understand the reasons behind this behavior, and we discovered they are due to repeated occurrences of timer interrupts because the simulator uses the host real-time clock (obtained, e.g. using rdtsc on x86 hosts) to trigger alarms. In this situation, the timer interruption is raised when a given wallclock Δ𝑡 has elapsed, which finally ends up calling the target Linux update_cfs_group scheduling function over and over. To exemplify the behavior, we have added in the cache model a delay using a for loop iterating 𝑛 times, which leads to the results given Figure 4. This can be mitigated by forcing QEMU to use a clock of its own 3 in which case only 1830 irqs are raised, for any added delay. However, in that case it cannot run vCPUs in parallel, which greatly hinders performance. As we are aiming to show that our cache model behaves as expected and compare its execution speed to an existing model or to vanilla QEMU, we must have the exact same instructions executed to be fair. To that end, we will use either bare metal software that does not program alarms or QEMU user-mode, both taking benefit of the host parallelism while having no dependency on simulation duration.

IMPLEMENTATION 5.1 QEMU TCG Plugins

The Tiny Code Generator (TCG) plugin is a feature introduced in QEMU version 4.2 [START_REF] Cota | Cross-ISA machine instrumentation using fast and scalable dynamic binary translation[END_REF]. It provides an API to write plugins in order to facilitate code instrumentation. Information is given through the plugin API at translation and execution times. It can be retrieved each time a block is translated or/and executed, or for every instruction and memory access that is executed. By instrumenting all instructions and memory accesses executed by the target through the TCG plugins, we have access to all the information we need to simulate a cache. Figure 5 shows the simplified plugin mechanism of QEMU's TCG. First, we have to indicate to which events proposed by the TCG plugin API we want to subscribe. It can be an instruction or a memory access execution or a translation block translation/execution. Then, during QEMU execution, each time the event occurs the main loop will send information to the plugin thanks to callback functions. Finally, inside the plugin, we can add our own code to do what we want. For this work, we do not consider the data-cache and solely focus on the instruction cache. TCG plugins are independent from the simulated architecture and thus can be used with all the different targets QEMU supports. In our case, we decided to work only with RISC-V but the approach we propose can be applied to all the other architectures as well.

Cache Simulation at Tb Granularity

The TCG plugin API is limited but offers the possibility to retrieve the target instructions belonging to a Tb during translation. During that stage, we have two options. The first one, implemented by the existing instruction cache plugin since 2021 [START_REF] Mandour | Cache Modelling TCG Plugin[END_REF], is to register a callback function for every instruction of the Tb. The function will be called right before instruction execution, once all its operand values known. Being at instruction granularity, the cache model is invoked each time an instruction is executed, which results in an important simulation time overhead. The second one is to register a callback each time a Tb is executed, which fits well the DBT principles, this is the option that we chose as introduced in Section 4. During the translation stage, thanks to the API of the TCG Plugins, we can parse the Tb and have access to the instructions it contains. Our plugin records in an array the addresses of the instructions that might produce a miss, and this array is given as argument to the function of the API that registers a callback for the Tbs. Then, just before a Tb is executed, the callback traverses the array and performs cache simulation on the sequence of addresses at once.

Using an array to store a subset of instructions requires dynamic memory allocation. However, thanks to the DBT mechanism, the Tbs are put in a cache to be reused. Even if it is application dependant, the percentage of Tb reuse is really high, so the per Tb memory allocation is negligible.

EXPERIMENTS

To evaluate our solution, we use two benchmarking suites: Poly-Bench/C for uniprocessor, and PARSEC for multiprocessors. The PolyBench/C suite is run with the MEDIUM inputs and the PARSEC suite with the LARGE inputs. We decided to work with the RISC-V target. We used QEMU in two different modes: user-mode, in which system calls and signals are handled by the host OS, and full-system, which is a Linux-capable target. QEMU, Linux and the machine we run on are not time deterministic. We therefore run our experiments 20 times to measure performances, which, by the virtue of the central limit theorem and given the standard deviation we obtain, gives a high confidence in the computed execution time mean.

In the following figures, cache refers to QEMU existing cache plugin, cacheTB refers to QEMU with a plugin that implements our solution and vanilla refers to QEMU without any plugin.

Error due to exceptions in Tbs

Because early internal Tbs exits can occur, we need to make sure the error it incurs is negligible. To that aim, we choose 2 programs that we run on a single vCPU: a boot of Linux that finishes just before printing the login prompt, and the LU decomposition of a 2048 × 2048 matrix from the parsec. We measure the total number of executed Tbs, the number of internal Tbs exits, the total number of instructions and the number of instructions that have been accounted for while they shouldn't. This is reported in Table 1. The average number of target instructions per Tb is 5.51 for Linux boot, and 16.61 for LU, which outlines the differences between both workloads. Although the error in number of instructions is negligible, we made a change in QEMU in which each memory access ends a Tb. For Linux boot, the number of executed Tbs almost quadruples, to 124,097,849, while the total simulation time increases by 50%. This is not a surprise, as in that workload, 44% of the instructions are memory accesses. In conclusion, our approach to consider all instructions that belong to a Tb as executed makes sense. It is not useful to mitigate it, since it costs a lot in simulation time while having zero impact on the cache statistics.

Statistics validation using QEMU existing cache plugin

We report in this section how we validated the statistics produced by our plugin. Although we thoroughly tested the behavior of our cache model against QEMU's current one, all the experiments were done with the following arbitrary instruction cache configuration: 8-way, 32-set, 64 bytes per line and LRU. . the existing cache plugin. We obtained the exact same number of instructions and instruction misses with both plugins for the whole benchmark. Thus we have been able to confirm the good behavior of our model in monoprocessor.

6.2.2 Multicore: PARSEC Suite. In order to validate our cache model without the issues raised in Section 4 when using Linux, we run the PARSEC suite in QEMU user-mode. Figure 6 shows the number of instruction misses for the number of virtual CPUs from 1 to 128 by integer power of 2 steps for the lu_cb PARSEC benchmark. The vertical black lines on the figure represent the range of statistics we did for 20 executions. We obtained for each of the vCPUs a mean of the statistics (number of instructions and number of instruction misses) that are really close from the ones of the cache plugin and the values are stable across the 20 executions. For all the PARSEC benchmarks that we used, we observed similar histograms than the one of lu_cb.

Simulation time: vanilla vs cache vs cacheTB

In this section, we compare the simulation time of QEMU cache plugin and our own cacheTB plugin. Figure 7 compares the simulation time of the benchmarks of the PolyBench/C suite in bare metal. Figure 8 compares the simulation time of the PARSEC program lu_cb run in user-mode with 1 to 128 virtual CPUs. Instead of showing the simulation time histograms for all PARSEC programs that we used and because they have all the same shape, we represent Figure 9 the simulation time using 128 virtual CPUs only for all programs (on a 2 GHz 64-core AMD EPYC 7702P PowerEdge R6515 server), outlining the scalability of our approach, QEMU in itself being fairly scalable [START_REF] Badaroux | To Pin or Not to Pin: Asserting the Scalability of QEMU Parallel Implementation[END_REF]. Table 2 summarizes all the different mean speedup/overhead ratios. According to the different types of experiments, our plugin is 7 to 10 times faster than the cache plugin. Regarding the overhead of using a plugin compared to QEMU vanilla, the cache plugin degrades much more the simulation time than our plugin (up to 60 times slower than QEMU vanilla with the PARSEC Suite compared to only 10 times slower with our model).

CONCLUSION

Although adding instrumentation in functional simulation will for sure degrade simulation time, taking benefit of the per block nature of the translation in DBT allows to define a strategy suited for instructions caches that minimizes the performance overhead. We raised two issues, one due to the DBT itself that we could mitigate, the other linked to time handling that we circumvent, but didn't solve per se. We implemented our approach in QEMU and validated against its existing instruction accurate cache model. For all the experiments we did, we obtained identical statistics, but at far better performances. Overall, fast instruction cache simulation isn't that easy, and even though doing it at scale for full-system is possible, asserting the validity of the metrics that are produced is still a work in progress.

Figure 1 :

 1 Figure 1: Outline of the dynamic binary translation process.

Figure 2 :

 2 Figure 2: Example of static hit/miss decision within a Tb.

Figure 3 :

 3 Figure 3: Stopped Tb execution due to a page-fault.

Figure 4 :

 4 Figure 4: Number of raised timer interrupts as a function of an arbitrary cache model simulation delay.

Figure 5 :

 5 Figure 5: Simplified representation of the QEMU TCG plugins mechanism.

Figure 6 :

 6 Figure 6: Number of instruction misses for lu_cb (log scale on 𝑥-axis).QEMU, Linux and the machine we run on are not time deterministic. We therefore run our experiments 20 times to measure performances, which, by the virtue of the central limit theorem and given the standard deviation we obtain, gives a high confidence in the computed execution time mean.In the following figures, cache refers to QEMU existing cache plugin, cacheTB refers to QEMU with a plugin that implements our solution and vanilla refers to QEMU without any plugin.

2 m m 3 m m a d i a ta x b ic g c h o le s k y c o rr e la ti o n c o v a ri a n c e d e ri c h e d o it g e n fd td -2 d g e m m g e m v e r g e s u m m v g ra m s c h m id t h e a t- 3 Figure 7 :

 37 Figure 7: Simulation time of the PolyBench/C programs (log scale on 𝑦-axis).

Figure 8 :

 8 Figure 8: Simulation time of lu_cb (log-log scale).the existing cache plugin. We obtained the exact same number of instructions and instruction misses with both plugins for the whole benchmark. Thus we have been able to confirm the good behavior of our model in monoprocessor.

Figure 9 :

 9 Figure 9: Simulation time of the PARSEC programs on 128 vCPUs (log scale on 𝑦-axis).Table2summarizes all the different mean speedup/overhead ratios. According to the different types of experiments, our plugin is 7 to 10 times faster than the cache plugin. Regarding the overhead of using a plugin compared to QEMU vanilla, the cache plugin degrades much more the simulation time than our plugin (up to 60 times slower than QEMU vanilla with the PARSEC Suite compared to only 10 times slower with our model).

2

  Bare metal programs do not suffer from this.

	# Insns in translation block
	0x7f1ffbe5692c : auipc	a5,237568
	0x7f1ffbe56930 : ld	a5,-612(a5)
	0x7f1ffbe56934 : sb	s0,0(a5)
	0x7f1ffbe56938 : ld	ra,8(sp)
	0x7f1ffbe5693a : auipc	a5,270336
	0x7f1ffbe5693e : sb	s0,1470(a5)	6 instructions
	0x7f1ffbe56942 : ld	s0,0(sp)
	0x7f1ffbe56944 : addi	sp,sp,16
	0x7f1ffbe56946 : ret
	# Executed insns until page fault
	0x7f1ffbe5692c : auipc	a5,237568
	0x7f1ffbe56930 : ld	a5,-612(a5)
	0x7f1ffbe56934 : sb	s0,0(a5)	←-first page-fault
	# New translation block after return from handler
	0x7f1ffbe56934 : sb	s0,0(a5)
	0x7f1ffbe56938 : ld	ra,8(sp)
	0x7f1ffbe5693a : auipc	a5,270336
	0x7f1ffbe5693e : sb	s0,1470(a5)
	0x7f1ffbe56942 : ld	s0,0(sp)
	0x7f1ffbe56944 : addi	sp,sp,16	3 instructions
	0x7f1ffbe56946 : ret
	# Executed insns until new page fault
	0x7f1ffbe56934 : sb	s0,0(a5)
	0x7f1ffbe56938 : ld	ra,8(sp)
	0x7f1ffbe5693a : auipc	a5,270336
	0x7f1ffbe5693e : sb	s0,1470(a5) ←-second page-fault
	# Again, new translation block after return from handler
	0x7f1ffbe5693e : sb	s0,1470(a5)
	0x7f1ffbe56942 : ld	s0,0(sp)
	0x7f1ffbe56944 : addi	sp,sp,16
	0x7f1ffbe56946 : ret
	# Executed insns until branch, no page-fault on ld
	0x7f1ffbe5693e : sb	s0,1470(a5)
	0x7f1ffbe56942 : ld	s0,0(sp)
	0x7f1ffbe56944 : addi	sp,sp,16
	0x7f1ffbe56946 : ret

Table 1 :

 1 Measure of the error due to early Tb exits.

		Nb of Tbs Nb of early Nb of executed	Wrongly
		executed	Tbs exits	insns counted insns
	Boot	35,704,017	254	196,831,388	669
	LU	989,522,360	10,447	16,439,546,310	6098

Table 2 :

 2 Mean simulation time ratios.

		PolyBench/C PARSEC
	Speedup cache to cacheTB	10.87	7.18
	Overhead cache to vanilla	23.67	59.85
	Overhead cacheTB to vanilla	2.07	10.16

To avoid ambiguity, we consistently use line for consecutive cache elements (a.k.a block) and block for what refers to translated target instructions.

ACKNOWLEDGMENTS

We would like to thank the partners of the ANR Rakes project and acknowledge the financial support of the French Agence Nationale de la Recherche (ANR-18-CE25-0017, https://anr.fr/Project-ANR-18-CE25-0017).


Maxence Denu
email: maxence.denu@grenoble-inp.fr

Pierre David

Aurélie Landry
email: aurelie.landry@univ-grenoble-alpes.fr

Fabien Mangione

Emerging digital technologies to support circular manufacturing systems implementation: a literature review

Keywords: circular economy, emerging technologies, industry 4.0, manufacturing, industrialization

Introduction

Faced with the repeated warnings issued by the IPCC (Intergovernmental Panel on Climate Change) about the emergency of the climate situation, Circular Economy is one of the solutions enlightened by diverse economic, political, and academic actors. In opposition to the linear economy (take-make-usedispose), Circular Economy could be defined as "an economic system that replaces the 'end-of-life' concept with reducing, alternatively reusing, recycling, and recovering materials in production/distribution and consumption processes. " [START_REF] Kirchherr | Conceptualizing the circular economy: An analysis of 114 definitions[END_REF].

To initiate the implementation of these circular economy strategies, it is necessary for new production systems to take place to allow the return of products, their disassembly, refurbishment or improvement. Numerous scholars predict that emerging technologies will be a key element to unlock the potential of a circular economy [START_REF] Bressanelli | Towards the Smart Circular Economy Paradigm: A Definition, Conceptualization, and Research Agenda[END_REF]- [START_REF] Rosa | Assessing relations between Circular Economy and Industry 4.0: a systematic literature review[END_REF]. Through the study of the case of SMEs, Prieto Sandoval [START_REF] Prieto-Sandoval | Key strategies, resources, and capabilities for implementing circular economy in industrial small and medium enterprises[END_REF] particularly indicates that the modernization of manufacturing systems is one of the major levers to achieving circularity in this type of company.

Considering the three pillars of sustainable development (economic, social and environmental), emerging technologies have important potential positive effects, but many risks must be anticipated from the adoption of these technologies to ensure social and environmental sustainability in particular [START_REF] Bonilla | Industry 4.0 and Sustainability Implications: A Scenario-Based Analysis of the Impacts and Challenges[END_REF]. The impacts of these technologies, when not implemented with a systemic and sustainable vision, could have important negative effects on several pillars. The consequences in terms of sustainability of Industry 4.0 also should not be taken for granted because it depends on the environment and conditions in which firms operate [START_REF] Ghobakhloo | Industry 4.0, innovation, and sustainable development: A systematic review and a roadmap to sustainable innovation[END_REF].

Nevertheless, the majority of works on the connection between CE and Industry 4.0 have primarily been conceptual or exploratory [START_REF] Jabbour | Better together": Evidence on the joint adoption of circular economy and industry 4.0 technologies[END_REF], and only few studies provide technical demonstrations of the observed contributions of these emerging technologies in the context of circularity [START_REF] Kerin | A review of emerging industry 4.0 technologies in remanufacturing[END_REF]. Indeed, the main published studies are surveys on the feelings of industrialists about the expected or obtained benefits of these future technologies [START_REF] Bai | Industry 4.0 technologies assessment: A sustainability perspective[END_REF], [START_REF] Kumar | Application of industry 4.0 technologies in SMEs for ethical and sustainable operations: Analysis of challenges[END_REF]. Kirchherr and van Santen [START_REF] Kirchherr | Research on the circular economy: A critique of the field[END_REF] particularly emphasize the need for industrialists to benefit from concrete, large-scale work to support the transition to these circular production systems.

2

Objectives of the communication

In this context, the cross-disciplinary CIRCULAR project aims to support the transition of industrial systems towards more circularity, through the transformation of end-of-life objects into high-valueadded products. Bringing together 12 research teams and 3 technological platforms in Grenoble, this project aims to participate in the deployment of methods of the industry of the future and emerging technologies to support companies in their transition to the production of objects that meet the needs of consumers, while reducing the use of resources necessary for their production. This communication will present the results of a literature review whose objective is to determine whether the technologies of the industry of the future will be adequate resources for the conduct of workshops in the context of circular economy, by adopting a multidisciplinary approach mixing engineering and human sciences. In this paper, we will present an analysis of circular economy practices, levers and barriers in the first part, before making a similar study on emerging technologies (notably with a focus on Industry 4.0 and related concepts). We will briefly point out the negative effects that could occur if the focus of a project is not based on a systemic and sustainable vision. Then, our objective will be to show the synergies and divergences between the two studied concepts, by comparing them with the one already highlighted in the literature. In particular, we will criticize the research methods adopted for the combination of circularity and emerging technologies, and show the obvious lack of empirical evidence on this topic. We will conclude by proposing a roadmap of topics to be studied around digital industrialization in the context of the circular economy. This literature review is part of a larger study, which aims to explore the concrete contributions of some of these emerging technologies (connected objects, autonomous logistics vehicles, information and communication technologies, collaborative robots, ...) through technical or experimental evidence to facilitate the implementation of the circular economy of the industry. In other words, we will aim to determine a methodology to study several technologies to draw up their benefits and potential risks associated in the context of their implementation of the circular economy, within a systemic analysis. Consideration of different perspectives including notably ethical and sustainable aspects will also be part of the study [START_REF] Kumar | Application of industry 4.0 technologies in SMEs for ethical and sustainable operations: Analysis of challenges[END_REF], [START_REF] Berrah | Performance industrielle : une évolution intégrant l'éthique dans le contexte de l'Industrie 4[END_REF].

Conclusion and perspectives

In a particularly difficult environmental and geopolitical context, the circular economy is perceived as a solution by many actors. To remove the many obstacles linked to circularity in a manufacturing context, the adoption of new technologies is seen as a real lever in this transition by academics. This communication will present a preliminary literature review with a multidisciplinary perspective, which will carefully focus on the environmental, social and economic pillars of sustainability, as well as on ethics. This exploratory research has raised a number of issues that require more study. Following objectives of this research project will be to promulgate empirical evidence, both technical and experimental, to determine the impact of emerging technologies in circular production systems.


Roberto Boccagna

Maurizio Bottini

Massimo Petracca

Matteo Di Giorgio

Alessia Amelio

Guido Camata

Prof Alan O'connor

Matteo Di

Unsupervised Methods for Railway Bridge Structural Monitoring

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

INTRODUCTION

Artificial Intelligence (AI) and especially Deep Learning (DL) algorithms are widely employed in the Civil and Structural Engineering context since large amounts of data can now be transmitted and processed in a relatively short time [START_REF] Tibaduiza Burgos | Damage identification in structural health monitoring: A brief review from its implementation to the use of data-driven applications[END_REF]]. Systems for early anomaly detection draw from a broad ensemble of Machine Learning (ML) techniques, the final purpose being the ability to automatically recognize possible damage that may affect a given structure so that appropriate measures can be taken in time.

Although most data analysis methods are well established, the scientific community has yet to agree on the sequence of processes to be performed, starting from the acquisition of raw data via sensors to get reliable predictions. This is especially a pivotal issue in Structural Health Monitoring (SHM), although general consensus exists regarding AI, with specific reference to the algorithms worth adopting in unsupervised or supervised cases. In the SHM framework, we rely on unsupervised algorithms whenever the training data come from structures in healthy conditions, as the dedicated program has no way to distinguish between undamaged and damaged scenarios. In such cases, using the data available, AI tries to infer the under-lying process that produces the observed output, hence learning to identify data that is statistically incompatible with that generation process, as the paradigm of SHM is that any structural change is reflected in the mechanism that produces the observed data. The dedicated literature recommends autoencoder neural networks as the type of algorithms best suited for accomplishing this task. Conversely, supervised algorithms are adopted whenever the training set also contains data corresponding to damage scenarios that may be generated by running dynamic numerical simulations using the Digital Twin (DT) of the physical structure under examination. In this case, the literature advises several multiclass DL algorithms for optimal performance [START_REF] Farrar | K. structural health monitoring: A machine learning perspective[END_REF]].

In the present work, we have restricted our focus to the unsupervised case, using the accelerometric time series we obtained by running railway bridge dynamics for a Finite Element Method (FEM) created using Asdea's Scientific ToolKit for OpenSees (STKO) advanced GUI for OpenSees [START_REF] Petracca | Stko user manual[END_REF]] for both the training and for the test phase. This was done with the intent to establish a standard procedure for managing strictly data-driven contexts, from data collection to output emission, with the idea of providing a robust solution for application in real monitoring cases.

The paper is organized as follows. Section 2 presents the bridge model and provides information on the simulation set-up. Section 3 explains how raw data are pre-processed and how the features from the extracted raw data are fed as inputs to the neural network. Section 4 is dedicated to the presentation of the results, and in particular, we show how the algorithm correctly classifies undamaged and damaged data series. Section 5 examines the conclusions drawn and suggests further developments.

THE MODEL

STKO modelization

The model used for testing our baseline anomaly detection solution was created using the STKO interface for OpenSees, and reproduces a steel truss railway bridge with riveted connections consisting of 3 spans (Figure 1). The structure is approxi-mately 93 m long, 5 m wide, and has two piers that are 4.8 m high. Each lateral span measures ∼ 28 m, while the central one is ∼ 35 m long. Further information about the model is provided below:

-the braces were modeled using T-profile truss elements -piers were modeled using 4-node shell elements (ASDShellQ4) -the ballast was modeled using springs -tracks were modeled using IPE beam elements -both top and bottom chords were modeled using double-T beam elements -verticals were modeled using IPE300 beam elements -diagonals were modeled using double-C beam elements Regarding boundary conditions, the piers were fixed to the base; one of the two abutments blocks linear displacements along the x, y and z directions and rotations about the x axis, while the other abutment blocks linear displacements along the y and z axis. The spans are attached to the piers by means of rigidLinks OpenSees elements.

The monitoring system of the bridge consists of a network of 12 triaxial accelerometric sensors, which sample at a rate of 1 kHz, located at highly representative nodes as depicted in Figure 1. Damage scenarios can be introduced by reducing the elastic moduli of 32 different structural elements by modifying them as a percentage from 0 to 100, giving the possibility to determine which elements to damage. Figure 2 shows the elements for which the user can reduce the elastic modulus, consisting of 8 diagonals in compression, 8 diagonals in tension, 8 chords in compression and 8 chords in tension.

Numerical Simulations

Numerical simulations were performed as follows. A single run refers to the passage of one train across the railway bridge which generates a single instance in terms of the AI algorithm. Partially following [START_REF] Parisi | Automated location of steel truss bridge damage using machine learning and raw strain sensor data[END_REF]'s strategy, trains were modeled using their mass and velocity, represented as random variables extracted from lognormal distributions with parameters respectively µ mass = 62 ton, σ mass = 5 ton and µ velocity = 8.33 m/s, σ velocity = 1 m/s, while the length of the trains is constant and fixed to 40 m. A single run lasts for T = 50 s, regardless of the exact moment the caboose of the train crosses the last bridge ele-ment, to obtain consistent time windows for processing. The output of the STKO program consists of accelerometric time series considering the components along x, y and z axis for each of the 12 control points, for a total of 36 time series for each run.

THE ALGORITHM

Feature Extraction

Univariate analysis was applied to each time series to aggregate the information embedded in the entire run through a minimal set of parameters. This is the first step in eliminating redundant data. A collection of temporal, spectral and statistical features was extracted using the Python TS-FEL library [Barandas et al. ( 2020)] as suggested in [START_REF] Buckley | A feature extraction & selection benchmark for structural health monitoring[END_REF]]1 . In particular, a set of 163 features was selected for each of the 36 accelerometric series obtained for each run, considering the whole set of default features the library computes after removing a large number of Fourier coefficients to discard frequencies higher than 30 Hz. The features obtained were then arranged in a matrix A of dimensions N runs × (N features × N sensors ), where N features = 163 and N sensors = 36 in this case.

Feature Reduction

The second step for database reduction consists of the simple application of a linear PCA algorithm at fixed retained variance.

Training Phase

Dataset standardization is first performed on A in such a way that

A i j → A i j -µ A j σ A j for any i, j (1)
where µ j and σ j are, respectively, the mean and the standard deviation computed along the jth column of A.

After covariance matrix C = 1 N runs -1 A ⊤ A is computed, the spectral decompo- sition of C = OΛO ⊤ is performed, being Λ = diag (λ 1 , .
. . , λ N features ×N sensors ) the matrix containing the eigenvalues of C in non-increasing order 2 . A is then projected along the first n components, which retain a given amount of the total variance (99% in our case):

X = AP n , P n = O [ : , : n] , (2)
where (we used Python notation) P n is the N runs × n matrix containing the first n columns of O.

Test Phase

Once the algorithm is trained, new data are standardized as in (1) using the same µ A j s and σ A j s computed over the matrix A. Then, the dataset is projected along the n eigendirections, multiplying by matrix P n as in (2).

The Neural Network

The rows {x i } N runs i=1 of matrix X represent single input instances for the neural network. The first algorithm tested was a vanilla autoencoder (AE) with input dimension n and just one hidden intermediate layer, which represents the latent space and has dimension d = ⌊n⌋ + 1. This choice was motivated by the purpose of testing one of the most simple neural network architecture before performing an optimal hyperparameters choice through dedicated tools. The graph is fully connected, as depicted in Figure 3, and a hyperbolic tangent is used as an activation function.

In the training phase, the normalized sum of Euclidean distances between the inputs and outputs was adopted as a loss function. Let f θ ,φ (x i ) be the reconstructed version of input x i , where θ and φ indicate the collection of parameters of the encoding and decoding portions respectively. For the reconstruction error, we take the following quantity:

e θ ,φ (x i ) = x i -f θ ,φ (x i ) 2 .
(3)

The loss function is obtained by averaging the whole training dataset:

L AE (X) = 1 N runs N runs ∑ i=1 e θ ,φ (x i ) , (4)
2 C is symmetric positive semi-definite, its eigenvalues are non-negative, and its eigenvectors form an orthogonal basis for R N features ×N sensors . so that biases and weights converge to the values that best allow the network to get outputs as similar as possible to the inputs belonging to the training set in the ℓ 2 norm. The second neural network adopted was a variational autoencoder (VAE), whose architecture is described in Figure 4. Unlike the simple vanilla autoencoder, the latent representation is not deterministic, as the latent vector is sampled from a multivariate, Gaussian distribution with mean vector µ and diagonal covariance matrix σ I. Let p θ (z) equal the prior probability of obtaining the latent vector, z, p θ (z | x i ) equal the posterior distribution of z given x i , and p θ (x i | z) be the conditional probability of x i given z, where θ again represents the collection of encoder parameters. The true posterior distribution is generally intractable and is then approximated by a distribution of q φ (z | x i ) (φ the collection of parameters for the decoder), which is chosen to be Gaussian, and represents the probability of observing the output x i given the latent variable z. The evidence lower bound (ELBO) is typically adopted as a loss function for VAEs, being a mixture of cross-entropy between the original and reconstructed dataset and Kullback-Leibler divergence that measures the functional distance between the true prior and the approximated posterior (see [START_REF] Kingma | Auto encoding variational bayes[END_REF] for derivation):

L VAE θ ,φ (X) = 1 N runs N runs ∑ i=1 ℓ θ ,φ (x i ) , (5)
ℓ θ ,φ (x i ) = E q φ (z|x i ) [log p θ (x i | z)] -KL q φ (z | x i ) ∥p θ (z) . ( 6
)
The first term at the right-hand side of ( 6) can be estimated through a reparametrization trick, while Kullback-Leibler divergence assumes a simple expression after forcing p θ (z) to be a standard Gaussian on R d . For details, we again recommend consulting [START_REF] Kingma | Auto encoding variational bayes[END_REF]].

For the sake of notational simplicity, hereafter, we will use e i ≡ e θ ,φ (x i ).

Statistics of the Reconstruction Error

From a probabilistic point of view, the collection of reconstruction errors {e i } N runs i=1 represents a set of independent identically distributed random variables on R + , as single train passages are mutually independent. The resulting statistics are a generalized chi-square distribution because components of x i and its reconstructed counterpart are, in principle, correlated. Various methods can be applied for deducing probability distribution p (e i ) from data, thus establishing a threshold for distinguishing undamaged from damaged data in the inference phase. Once threshold α ∈ (0, 1) is established, a given instance can be considered anomalous when its reconstruction error e is such that p (e) < α. We used kernel density estimation (KDE) method for inferring the pdf of the underlying process. In our application, the bandwidth h of the kernel, which is the main parameter of the method, was estimated via Silverman's rule of thumb h = 0.9 min σ e , IQR 1.34 N -1 5 runs , where σ e is the standard deviation of reconstruction errors, and IQR represents the interquartile range. The threshold for acceptability α was fixed to 0.005, a value that finds its validity a posteriori. Note that, in this initial phase of the construction of our framework, we expect a structure in a healthy state to produce false alarms with probability α. To further perfect the framework, a method for alarm validations is thus required and will be introduced in the future.

RESULTS

The training of the two algorithms was performed using N runs = 500 train passages, consequently inferring the reconstruction error pdf and obtaining values for e min and e max such that p (e < e min ) < α and p (e > e max ) > α. Figure 5 shows the shape of the pdfs deduced from the training data for AE and VAE, alongside the corresponding frequency histogram.

We initially focused on the predictive capability of the two algorithms as anomaly detectors by introducing damage on node DC 1 . Specifically, FEM dynamics were run after the elastic modulus of the mentioned beam element was reduced by a factor of 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, and 0.10, considering 10 runs for each of these percentages of damage. In addition, we ran 25 further simulations for the undamaged bridge in order to check if the trained algorithms would classify the corresponding inputs as non-anomalous. Figure 6 shows the collection of reconstruction errors for the training set (from 0 to 499, blue colored) and for the test set corresponding to the 25 undamaged data sets (500 to 524, green colored if regular, red colored if anomalous), using thresholds (the dashed, gray lines) corresponding to the values of e min and e max . For the AE, the mean reconstruction error of the training set is µ AE training ≃ 1.00 • 10 -2 , compared with a value of µ AE, undamaged test ≃ 9.53 • 10 -3 for the set of the 10 undamaged tests. For the VAE, we ob- Figure 7 shows the collection of reconstruction errors for the training set and for the aforementioned damage configurations (enumerated from 500 to 559, red colored, in increasing damage order), adopting the log-scale for the y-axis since re- construction errors spread on various scales. Notice that both algorithms recognize all the "damaged" instances as anomalous, although the reconstruction error appears to be non-monotone in the elastic modulus reduction. It is also worth noting that VAE is, in this case, more damage-sensitive compared to the AE-based solution since the errors corresponding to the test set cover a larger interval. This is also a consequence of VAE's pdf appearing more peaked around its maximum. The computational costs are comparable for the two neural network autoencoders employed.

We performed a second test running numerical dynamics simulations for each of the 32 damage locations as in Figure 2. In order to have strong data consistency, we fixed the mass and the velocity of the train crossing the bridge to the mean values of 62 tons and 8.33 m/s and reduced the elastic moduli of beam elements by a factor of 0.02. Figure 8 again illustrates reconstruction errors for the training sets and for the 32 different damage scenarios. The evidence supports the conclusion that our solutions correctly classify all the instances given as inputs as anomalous.

CONCLUSIONS

This work represents a simple yet reliable framework for SHM in the full data-driven case. The results obtained so far are promising as our detector succeeded in recognizing healthy states and classifying the various configurations of damage types/severity as anomalous. Despite this, we were unable to define a connection between elastic modulus reduction and the output obtained. In the future, more effort will be made towards selecting features and exploring the space of the neural network's hyperparameters to enable our solution to estimate the actual damage level from the value of the reconstruction error obtained. After that, we can move towards considering the supervised case and trying to identify possible damage locations from test data. The introduction of a digital twin of the monitored structure would also improve the overall accuracy of the solution.

Figure 1 :

 1 Figure 1: Schematic representation of the bridge and sensor positions.

Figure 2 :

 2 Figure 2: Locations of the 32 elements whose elastic modulus the user can reduce. Legend: DC (diagonal in compression), DT (diagonal in tension), CC (chord in compression), CT (chord in tension).

Figure 4 :

 4 Figure 3: Schematic representation of the vanilla autoencoder used.

Figure 5 :

 5 Figure 5: The pdf (the orange lines) inferred from training data in the case of AE (upper) and VAE (bottom). Histograms (in blue) were obtained using 20 bins.

Figure 6 :

 6 Figure 6: The reconstruction error for training and non-damaged test sets for AE (top) and VAE (bottom).

Figure 7 :

 7 Figure 7: The reconstruction error for training and damaged test sets as indicated in the text for AE (top) and VAE (bottom).

Figure 8 :

 8 Figure 8: The reconstruction error for training and the second damaged test sets as indicated in the text for AE (top) and VAE (bottom).

See https://tsfel.readthedocs.io/en/latest/ for a complete list of computed features.